# Clustering for directed graphs using parametrized random walk diffusion
kernels

Harry Sevi
Centre Borelli, Université Paris-Saclay, CNRS, ENS Paris-Saclay
<EMAIL_ADDRESS>Matthieu Jonckheere
CNRS, LAAS
<EMAIL_ADDRESS>Argyris Kalogeratos
Centre Borelli, Université Paris-Saclay, CNRS, ENS Paris-Saclay
<EMAIL_ADDRESS>
###### Abstract

Clustering based on the random walk operator has been proven effective for
undirected graphs, but its generalization to directed graphs (digraphs) is
much more challenging. Although the random walk operator is well-defined for
digraphs, in most cases such graphs are not strongly connected, and hence the
associated random walks are not irreducible, which is a crucial property for
clustering that exists naturally in the undirected setting. To remedy this,
the usual workaround is to either naively symmetrize the adjacency matrix or
to replace the natural random walk operator by the teleporting random walk
operator, but this can lead to the loss of valuable information carried by
edge directionality. In this paper, we introduce a new clustering framework,
the _Parametrized Random Walk Diffusion Kernel Clustering_ (P-RWDKC), which is
suitable for handling both directed and undirected graphs. Our framework is
based on the diffusion geometry (Coifman & Lafon, 2006) and the generalized
spectral clustering framework (Sevi et al., 2022). Accordingly, we propose an
algorithm that automatically reveals the cluster structure at a given scale,
by considering the random walk dynamics associated with a parametrized kernel
operator, and by estimating its critical diffusion time. Experiments on $K$-NN
graphs constructed from real-world datasets and real-world graphs show that
our clustering approach performs well in all tested cases, and outperforms
existing approaches in most of them.

##  1 Introduction

Clustering is a fundamental unsupervised learning task, whose aim is to
analyze and reveal the cluster structure of unlabeled datasets, and has
widespread applications in machine learning, network analysis, biology, and
other fields (Kiselev et al., 2017; McFee & Ellis, 2014). Clustering for data
represented as a graph has been formulated in various ways. A well-established
one consists in minimizing a functional of the graph-cut (Von Luxburg, 2007;
Shi & Malik, 2000), leading to the spectral clustering (SC) framework and
various algorithms. SC is simple and effective, but has important limitations
(Nadler & Galun, 2006) that make it unreliable in a number of non-rare data
regimes. Overcoming these limitations is where the focus of the machine
learning community is (Tremblay et al., 2016; Zhang & Rohe, 2018; Dall’Amico
et al., 2021; Sevi et al., 2022).

A well-studied clustering approach for high-dimensional data, which is of
particular interest for this work, is using the operator associated with a
random walk (or, in other terms, with a Markov chain), called _random walk
operator_. Meilă & Shi (2001) viewed the pairwise similarities between
datapoints as edge flows of a Markov chain and proposed an ergodic random walk
interpretation of the spectral clustering. However, and beyond spectral
clustering, the first to conceive the idea of turning the distance matrix
between high-dimensional data into a Markov process were Tishby & Slonim
(2000). More specifically, they proposed to examine the decay of mutual
information during the relaxation of the Markov process. During the relaxation
procedure, the clusters emerge as quasi-stable structures, and then they get
extracted using the information bottleneck method (Tishby et al., 2000). Azran
& Ghahramani (2006) proposed to estimate the number of data clusters by
estimating the diffusion time of the random walk operator that reveals the
most significant cluster structure. Lin & Cohen (2010) proposed then to find a
low-dimensional embedding of the data using the truncated power iteration of a
random walk operator derived from the pairwise similarity matrix.

Clustering high-dimensional data using a random walk operator through the lens
of the diffusion geometry (Coifman & Lafon, 2006; Coifman et al., 2005) has
also been investigated. Nadler et al. (2006) proposed a unifying probabilistic
diffusion framework based on the probabilistic interpretation of spectral
clustering and dimensionality reduction algorithms using the eigenvectors of
the normalized graph Laplacian. Given the pairwise similarity matrix built
from high-dimensional data points, they defined a distance function between
any two points based on the random walk on the graph called the diffusion
distance (Coifman & Lafon, 2006; Pons & Latapy, 2005) and showed that the low-
dimensional representation of the data by the first few eigenvectors of the
corresponding random walk operator is optimal under a certain criterion.
Recently, an unsupervised clustering methodology has been proposed based on
the diffusion geometry framework (Maggioni & Murphy, 2019; Murphy & Polk,
2022), where it has been demonstrated how the diffusion time of the random
walk operator can be exploited to successfully cluster datasets for which
$k$-means, spectral clustering, or density-based clustering methods fail.

Although there has been a constant development of clustering approaches which
are based on the random walk operator, all such efforts (like those mentioned
above) have been proposed for undirected graphs. Moreover, the motivation of
most of them is to run forward the random walk to avoid the costly
eigendecomposition of the graph Laplacian. The extension of clustering
approaches based on the random walk operator to digraphs is much more subtle
and challenging. As presented earlier, random walk-based clustering approaches
rely either on the eigenvectors of the random walk operator or on its iterated
powers. In the directed setting, the random walk is well-defined, but it is
not reversible in general like in the undirected case (Levin & Peres, 2017).
Consequently, the associated eigenvectors are possibly complex, which makes
their interpretation and use difficult in the context of clustering. A
possible workaround would be to consider the approach based on the iterated
powers of the random walk operator which would allow one to obtain a real-
valued embedding and thus avoid the use of complex eigenvectors.

However, a second much more subtle and problematic bottleneck arises: the
random walk’s irreducibility. While in the undirected case, the random walk is
irreducible, i.e. any graph vertex can be reached from any other vertex, this
is not the case in general for random walks on digraphs.To overcome the
irreducibility issue, either a symmetrization procedure is usually employed in
the case where the digraph is derived by a $K$-NN graph construction or the
original random walk operator is replaced by the operator of the _teleporting
random walk_ (Page et al., 1999). For these two workarounds, valuable
information can potentially be discarded. To address this problem, we present
a new clustering algorithm based on the diffusion geometry framework, which is
suitable for any digraph, either obtained by $K$-NN constructions or by real
digraphs representing asymmetric relationships between vertices (e.g. citation
graphs). Our algorithm stems from a new type of graph Laplacians (Sevi et al.,
2022) that is parametrized by an arbitrary vertex measure capable of encoding
digraph information, and from which we derive a random walk operator. Besides,
we can exploit the necessary diffusion time of this parametrized random walk
to successfully reveal clusters at different scales.

The contribution of this work is manifold: i) We propose a new similarity
kernel operator that we term the _random walk diffusion kernel_ , which
derives directly from the diffusion distance. ii) We generalize this kernel
through the parametrized random walk operator that we introduced, which allows
us to extend the diffusion distance to digraphs. iii) From there, we propose
our clustering algorithm on digraphs, the _parametrized random walk diffusion
kernel clustering_ , based on the parametrized random walk diffusion kernel
operator considered as a data embedding. iv) We propose a method to estimate
the diffusion time that reveals at best a given number of clusters, regardless
of the graph origin. Finally, we show that our approach is efficient on both
synthetic and real-world data as it outperforms existing methods in various
cases.

##  2 Background

##### Graph theory.

Let $\mathcal{G}=(\mathcal{V},\mathcal{E},w)$ be a (weighted) directed graph
(digraph), where $\mathcal{V}$ is the finite set of $N=|\mathcal{V}|$
vertices, and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the finite
set of edges. Each edge $(i,j)$ is an ordered vertex pair indicating the
direction of a link from vertex $i$ to vertex $j$. Any function
$\nu:\mathcal{V}\rightarrow\mathbb{R}_{+}$, associating a nonnegative value to
each graph vertex, can be regarded as a positive vertex measure. The edge
weight function $w:\mathcal{V}\times\mathcal{V}\rightarrow\mathbb{R}_{+}$
associates a nonnegative real value to every vertex pair: $w(i,j)>0$, iff
$(i,j)\in\mathcal{E}$, otherwise $w(i,j)=0$. A digraph $\mathcal{G}$ is
represented by a weighted adjacency matrix
$\mathbf{W}=\\{w(i,j)\\}_{i,j=1}^{N}\in\mathbb{R}_{+}^{N\times N}$, where
$w(i,j)$ is the weight of the edge $(i,j)$. We define the out-degree and the
in-degree of the $i$-th vertex by $\textstyle
d_{i}^{\textnormal{out}}=\sum_{j=1}^{N}w(i,j)$ and
$d_{i}^{\textnormal{in}}=\sum_{j=1}^{N}w(j,i)$, respectively. Also, the
function $\mathbf{D}_{\nu}=\operatorname{diag}(\nu)$ is a square diagonal
matrix with the elements of the input vector $\nu$ in its diagonal. .

##### Random walk fundamentals.

What we call in short a random walk on a weighted graph $\mathcal{G}$, is
defined more formally as a natural random walk on the graph, which is a
homogeneous Markov chain $\mathcal{X}=(X_{t})_{t\geq 0}$ with a finite state
space $\mathcal{V}$, and with state transition probabilities proportional to
the edge weights. The entries of the transition matrix
$\mathbf{P}=[\,p(i,j)\,]_{i,j\in\mathcal{V}}$ are defined by:

$p(i,j)=\mathbb{P}(X_{t+1}=j\,|\,X_{t}=i)=\frac{w(i,j)}{\sum_{z\in\mathcal{V}}w(i,z)}.$

Algebraically, the transition matrix $\mathbf{P}\in\mathbb{R}^{N\times N}$ can
be expressed as $\mathbf{P}=\mathbf{D}_{{\textnormal{out}}}^{-1}\mathbf{W}$,
$\mathbf{D}_{\textnormal{out}}=\mathbf{D}_{d^{\textnormal{out}}}$
(respectively,
$\mathbf{D}_{\textnormal{in}}=\mathbf{D}_{d^{\textnormal{in}}}$), whose
spectrum $\textnormal{sp}(\mathbf{P})\in[-1,1]$. For a strongly connected
digraph $\mathcal{G}$, the random walk $\mathcal{X}$ is irreducible. If, in
addition, the irreducible random walk $\mathcal{X}$ admits the aperiodicity
condition, then $\mathcal{X}$ is ergodic, and therefore as
$t\rightarrow\infty$, the probability measures, $\forall i\in\mathcal{V}$ and
$\delta_{i}\in\\{0,1\\}^{N\times 1}$ being the vector output of the Kronecker
delta function, converges towards a unique stationary distribution denoted by
the row vector $\pi\in\mathbb{R}_{+}^{N}$ (Brémaud, 2013). Within the
undirected setting: $d_{i}^{\textnormal{out}}=d_{i}^{\textnormal{in}}=d_{i}$,
where $d\in\mathbb{R}_{+}^{N\times 1}$ is the vector of the vertex degrees;
moreover, the stationary distribution is proportional to the vertex degree
distribution, i.e. $\pi\propto d$.

##### Parametrized graph Laplacians.

A notable example of parametrized operators that we use later, is the
generalized graph Laplacians from Sevi et al. (2022), which have been proposed
as a new type of graph operators for (but not restricted to) digraphs.

######  Definition 2.1.

Generalized graph Laplacians (Sevi et al., 2022). Let $\mathbf{P}$ be the
transition matrix of a random walk on a digraph $\mathcal{G}$. Under an
arbitrary positive vertex measure $\nu$ on $\mathcal{G}$, consider the
positive vertex measure $\xi=\nu^{{\mkern-1.5mu\mathsf{T}}}\mathbf{P}$. Let
also the diagonal matrices $\mathbf{D}_{\nu}=\operatorname{diag}(\nu)$ and
$\mathbf{D}_{\xi}=\operatorname{diag}(\xi)$. Then, we can define three
generalized Laplacians of $\mathcal{G}$ as follows:

$\displaystyle\text{{generalized random walk Laplacian}:}\ \ \ \
\mathbf{L}_{\textnormal{RW},(\nu)}$
$\displaystyle=\mathbf{I-(I+\mathbf{D}_{\xi/\nu})^{-1}(P+\mathbf{D}_{\nu}^{-1}P^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu})}$
(1) $\displaystyle\text{{unnormalized generalized Laplacian}:}\ \ \ \ \ \ \ \
\ \mathbf{L}_{(\nu)}$
$\displaystyle=\mathbf{\mathbf{D}_{\nu+\xi}-(\mathbf{D}_{\nu}P+P^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu}})$
(2) $\displaystyle\text{{normalized generalized Laplacian}:}\ \ \ \ \ \ \ \ \
\mathcal{L}_{(\nu)}$
$\displaystyle={\mathbf{\mathbf{D}_{\nu+\xi}^{-1/2}}\mathbf{L}_{(\nu)}\mathbf{\mathbf{D}_{\nu+\xi}^{-1/2}}}.$
(3)

$\mathbf{L}_{\textnormal{RW},(\nu)}$, $\mathbf{L}_{(\nu)}$, and
$\mathcal{L}_{(\nu)}$ correspond to the directed graph Laplacians in Chung
(2005) when the transition matrix $\mathbf{P}$ is irreducible and the vertex
measure $\nu$ is the ergodic measure $\pi$. These graph Laplacians are
parametrized by an arbitrary vertex measure has the modeling capacity to
encode the graph directionality using the random walk dynamics of the original
digraph.

##  3 Parametrized random walk operator on graphs

In this section, we introduce the parametrized random walk operator that we
propose, and we show how it is derived from the generalized graph Laplacian
(presented at the end of Sec. 2).

######  Definition 3.1.

Parametrized random walk operator (P-RW). Let $\mathbf{P}$ be the transition
matrix of a random walk on a digraph $\mathcal{G}$. Under an arbitrary
positive vertex measure $\nu$ on $\mathcal{G}$, consider the positive vertex
measure $\xi=\nu^{{\mkern-1.5mu\mathsf{T}}}\mathbf{P}$. Let also the diagonal
matrices $\mathbf{D}_{\nu}=\operatorname{diag}(\nu)$ and
$\mathbf{D}_{\xi}=\operatorname{diag}(\xi)$. Finally, let $\mathcal{X}_{\nu}$
be the random walk on $\mathcal{G}$ with the associated random walk operator
(transition matrix) $\mathbf{P}_{\nu}$ defined by:

$\mathbf{P}_{(\nu)}=(\mathbf{I}+\mathbf{D}_{\xi/\nu})^{-1}(\mathbf{P}+\mathbf{D}_{\nu}^{-1}\mathbf{P}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu}).$
(4)

It is easy to check that $\mathbf{P}_{(\nu)}$ is a transition matrix (see
Appendix A.1.2). Besides, the parametrized random walk $\mathbf{P}_{(\nu)}$ is
directly related to the generalized random walk Laplacian
$\mathbf{L}_{\textnormal{RW},(\nu)}$ as shown in the next proposition.

######  Proposition 3.1.

Let $\mathbf{P}$ be the transition matrix of a random walk on a digraph
$\mathcal{G}$. Under an arbitrary positive vertex measure $\nu$ on
$\mathcal{G}$. The generalized random walk Laplacian and the parametrized
random walk on $\mathcal{G}$ are related as follows:

$\mathbf{L}_{\textnormal{RW},(\nu)}=\mathbf{I}-\mathbf{P}_{(\nu)}.$

###### Proof.

We use Eq. 1 of Def. 2.1, and we combine with Eq. 4 of Def. 3.1. ∎

As stated in Sevi et al. (2022), the generalized random walk Laplacian,
$\mathbf{L}_{\textnormal{RW},(\nu)}$, is self-adjoint in
$\ell^{2}(\mathcal{V},\nu+\xi)$. Consequently, the generalized random walk
$\mathbf{P}_{(\nu)}$ is also self-adjoint in $\ell^{2}(\mathcal{V},\nu+\xi)$,
and the associated random walk $\mathcal{X}_{\nu}$ is ergodic (under the
aperiodicity condition) and admits the ergodic distribution $\pi_{\nu}$.
Therefore, we have a random walk operator $\mathbf{P}_{(\nu)}$ parameterized
by a vertex measure $\nu$ that allows the encoding of the random walk dynamic
of the original digraph.

##  4 Diffusion geometry for digraphs

In this section, we review the concept of diffusion geometry (Coifman & Lafon,
2006). We show: i) how its core feature, _the diffusion distance_ , can be
expressed as a Mahalanobis distance involving a specific kernel matrix, which
we call _Random Walk Diffusion Kernel_ (RWDK); ii) we show how the SC’s
diagonalization step can be thought of as a function applied to the spectrum
of the Laplacian and the conceptual connection with the RWDK, and finally iii)
propose a clustering algorithm based on the generalization of this kernel.

###  4.1 The Diffusion distance as a Mahalanobis distance

The seminal work by (Coifman & Lafon, 2006) introduced the _diffusion geometry
framework_ , which uses diffusion processes as basic tool to find meaningful
geometric descriptions for dataset. The framework can provide different
geometric representations of the dataset by iterating the Markov transition
matrix, which is equivalent to running forward the random walk. The key
element of the diffusion geometry is the diffusion distance (Coifman & Lafon,
2006; Pons & Latapy, 2005) defined as follows.

######  Definition 4.1.

Diffusion distance. Let $\mathbf{P}$ be the transition matrix of a reversible
random walk on an undirected graph $\mathcal{G}$, with an ergodic distribution
$\pi$. The diffusion distance at time $t\in\mathbb{N}$ between the vertices
$i$ and $j$ is defined by:

$d_{t}^{2}(i,j)=\|p_{t}(i,*)-p_{t}(j,*)\|_{1/\pi}^{2}.$

As noted in (Coifman & Lafon, 2006), the diffusion distance emphasizes the
cluster structure, if present. Next, we show that the diffusion distance can
be seen as a Mahalanobis distance.

######  Proposition 4.1.

Diffusion distance as Mahalanobis distance. Let $\mathbf{P}$ be the transition
matrix of a reversible random walk on an undirected graph $\mathcal{G}$, with
an ergodic distribution $\pi$. $\mathbf{P}$ admits the following
eigendecomposition
$\mathbf{P}=\bm{\Phi}\mathbf{D}_{\lambda}\bm{\Psi}^{{\mkern-1.5mu\mathsf{T}}}$.
The diffusion distance between vertices $i$ and $j$, $d_{t}^{2}(i,j)$ at a
given time $t\in\mathbb{N}$ can be expressed as the following Mahalanobis
distance:

$d_{t}^{2}(i,j)=(\delta_{i}-\delta_{j})^{{\mkern-1.5mu\mathsf{T}}}\mathbf{K}_{t}(\delta_{i}-\delta_{j}),$

where the similarity positive definite kernel matrix $\mathbf{K}_{t}$ is
defined by:

$\mathbf{K}_{t}=\mathbf{P}^{2t}\mathbf{D}_{d}^{-1}.$ (5)

The diffusion distance as a dissimilarity measure reveals a similarity kernel
matrix $\mathbf{K}_{t}$ that we call _Random Walk Diffusion Kernel_ (RWDK)
matrix, which is simply a power of the transition matrix $\mathbf{P}$
normalized by the vertex degrees. Consequently, using the diffusion distance
$d_{t}^{2}(x,y)$ is equivalent to using the RWDK matrix $\mathbf{K}_{t}$.

###  4.2 Random Walk Diffusion Kernel, Normalized Graph Laplacian and Spectral
Clustering

Spectral clustering (SC) is one of the most widely used clustering methods due
to its simplicity, efficiency, and strong theoretical foundation (Shi & Malik,
2000; Ng et al., 2002; Von Luxburg, 2007; Peng et al., 2015; Boedihardjo et
al., 2021). Given a fixed number of clusters $k$, SC consists of three main
steps: i) construct the graph Laplacian matrix; ii) compute the eigenvectors
associated with the $k$ smallest eigenvectors of the Laplacian matrix and
store them as columns in a matrix; iii) apply $k$-means on the rows of the
latter matrix, which are regarded as embedded representations of the
datapoints. We aim to highlight how the SC’s second step and the RWDK do have
similar characteristics. Let us consider the normalized graph Laplacian matrix
$\mathcal{L}$ (Chung & Graham, 1997) with eigendecomposition
$\mathcal{L}=\sum_{j=1}^{N}\vartheta_{j}\bm{\phi}_{j}\bm{\phi}_{j}^{{\mkern-1.5mu\mathsf{T}}}$,
ordered eigenvalues $0\leq\vartheta_{1}\leq...\leq\vartheta_{N}\leq 2$ and
eigenvectors $\\{\bm{\phi}_{j}\\}_{j=1}^{N}$. Computing the eigenvectors
associated to the $k$ smallest eigenvalues of $\mathcal{L}$ amounts to
applying a function $f_{1}$ on the spectrum matrix of $\mathcal{L}$ such that
$f_{1}(x)=1$ if $x\leq\vartheta_{k}$, and $f_{1}(x)=0$ otherwise, namely
$\mathbf{H}=f_{1}(\mathcal{L})=\sum_{j=1}^{k}\vartheta_{j}\bm{\phi}_{j}\bm{\phi}_{j}^{{\mkern-1.5mu\mathsf{T}}}$.
Thanks to the relation
$\mathcal{L}=\mathbf{D}^{\frac{1}{2}}(\mathbf{I}-\mathbf{P})\mathbf{D}^{-\frac{1}{2}}$,
selecting the $k$ smallest eigenvalues of $\mathcal{L}$ is thus equivalent to
selecting the $k$ largest eigenvalues of $\mathbf{P}$, namely
$\mathbf{H}=f_{1}(\mathcal{L})=\mathbf{D}^{\frac{1}{2}}f_{1}(\mathbf{I}-\mathbf{P})\mathbf{D}^{-\frac{1}{2}}$.
On the other hand, the RWDK $\mathbf{K}_{t}$ corresponds to applying a
function $f_{2}(x,t)=x^{2t}$ on the spectrum matrix of $\mathbf{P}$. As
$t\rightarrow\infty$,
$f_{2}(\mathbf{P},t)\rightarrow\pi\bm{1}^{{\mkern-1.5mu\mathsf{T}}}$, and
consequently
$\mathbf{K}_{t}=f_{2}(\mathbf{P},t)\mathbf{D}^{-1}\rightarrow[\textnormal{tr}(\mathbf{D}_{d})]^{-1}\bm{1}\bm{1}^{{\mkern-1.5mu\mathsf{T}}}$.
As $t$ increases, $f_{2}(\mathbf{P})$ becomes smoother because $f_{2}$ acts as
a soft truncation of the spectrum matrix of $\mathbf{P}$, which preserves
those eigenvectors of $\mathbf{P}$ that are associated with the largest
eigenvalues. However, as $t$ increases, $f_{1}$ acts as a hard truncation of
the spectrum matrix of $\mathcal{L}$, which preserves the $k$ eigenvectors
associated with the lowest ones up to $\vartheta_{k}$. Consequently, the RWDK
$\mathbf{K}_{t}$ is an adaptive alternative to $\mathbf{H}$, and the aim is to
determine the iteration time $t$ that best reveals the $k$ clusters the user
looks for.

###  4.3 Parametrized random walk Kernel clustering

We have shown that the diffusion distance is directly related to the RWDK and
that for a given diffusion time $t$, the RWDK is an alternative to computing
the eigenvectors of the graph Laplacian. However, the diffusion distance is
only settled for undirected graphs (Coifman & Lafon, 2006; Boyd et al., 2021).
Nevertheless, thanks to the parametrized random walk operator built to handle
digraphs, one can easily extend the diffusion distance for digraphs without
any restriction. This is established by the following definition.

######  Definition 4.2.

Parametrized diffusion distance. Let $\mathcal{X}$ be a random walk on digraph
$\mathcal{G}$, with transition matrix $\mathbf{P}$. Let $\nu$ be an arbitrary
positive vertex measure on $\mathcal{G}$, and $\xi$ be the vertex measure
defined by $\xi=\nu^{{\mkern-1.5mu\mathsf{T}}}\mathbf{P}$. Define the diagonal
matrices $\mathbf{D}_{\nu}=\operatorname{diag}(\nu)$ and
$\mathbf{D}_{\xi}=\operatorname{diag}(\xi)$. On $\mathcal{G}$, let
$\mathcal{X}_{\nu}$ be a random walk associated with the random walk operator
$\mathbf{P}_{(\nu)}$ parametrized by an arbitrary measure $\nu$ defined in Eq.
4 and ergodic measure $\pi_{\nu}$. Let
$p_{t,\nu}(i,*)=\delta_{i}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{P}_{\nu}^{t}$ be
the conditional probability vector given the vertex $i$ at the diffusion time
$t$. The parametrized diffusion distance between the vertices $i$ and $j$, at
a given diffusion time $t$, is defined by:

$d_{t,\nu}^{2}(i,j)=\|p_{t,\nu}(i,*)-p_{t,\nu}(j,*)\|_{1/\pi_{\nu}}^{2}=(\delta_{i}-\delta_{j})^{{\mkern-1.5mu\mathsf{T}}}\mathbf{K}_{t,\nu}(\delta_{i}-\delta_{j}),$
(6)

with $\mathbf{K}_{t,\nu}$ is the parametrized random walk diffusion kernel
(P-RWDK) defined as:

$\mathbf{K}_{t,\nu}=\mathbf{P}_{\nu}^{t}\mathbf{D}_{\nu+\xi}^{-1}.$ (7)

Once we have defined the parameterized diffusion distance $d_{t,\nu}^{2}$, we
are in position to propose the novel _parametrized random walk diffusion
kernel clustering_ (P-RWDKC) for digraphs, which is based upon the
parametrized RWDK $\mathbf{K}_{t,\nu}$ described in Alg. 1. As noted, P-RWDKC
has an algorithmic scheme similar to the SC algorithm except for the choice of
the operator and the fact that, our method computes the parametrized RWDK
$\mathbf{K}_{t,\nu}$ instead of computing the eigenvectors of the graph
Laplacian.

Algorithm 1 Parametrized Random Walk Diffusion Kernel Clustering (P-RWDKC)

Input: $\mathbf{W}\in\mathbb{R}^{N\times N}$: adjacency matrix, $k$: number of
clusters, $\nu$: vertex measure, $t_{d}$: diffusion time
Output: $\bm{V}_{t_{d}}$: graph $k$-partition for diffusion time $t_{d}$

1: Compute the parametrized random walk operator $\mathbf{P}_{(\nu)}$, see Eq.
4

2: Compute the parametrized random walk diffusion kernel
$\mathbf{K}_{(t_{d},\nu)}$, see Eq. 7

3: Consider each $x_{i}\in\mathbb{R}^{N}$, $i=1,...,N$, to be the embedding of
the $i$-th vertex, represented by the $i$-th row of
$\mathbf{K}_{(t_{d},\nu)}$, and apply a clustering method (e.g. $k$-means) to
all these vectors asking for $k$ clusters

4: Obtain the $k$-partition
$\bm{V}_{t_{d}}=\big{\\{}V_{t_{d}}\big{\\}}_{j=1}^{k}$ of the graph vertices
based on the clustering result of Step 3

5: return $\bm{V}_{t_{d}}$

##  5 The parametrized random walk diffusion kernel clustering in practice

In the previous section, we described a novel and simple clustering algorithm
for digraphs based on an operator parametrized by an arbitrary vertex measure
$\nu$ and a diffusion time $t$. To render our algorithm more flexible for
practical use, two aspects need to be treated: i) designing the vertex measure
$\nu$, ii) determining the appropriate diffusion time.

###  5.1 Designing the vertex measure

Designing the vertex measure is one of the major aspects of P-RWDKC as we aim
to capture with it the random walk dynamics of the original digraph in our
parametrized random walk operator. To do so, we propose a vertex measure
derived from the iterated powers of a random walk consisting of the forward
and backward digraph’s flow information. Specifically, the proposed vertex
measure can be parametrized by three parameters
$(t\\!\in\\!\mathbb{N},\,\gamma\\!\in\\![0,1],\,\alpha\\!\in\\!\mathbb{R})$
and is formally given by:

$\nu_{(t,\gamma)}^{\alpha}(i)=\left({\textstyle\frac{1}{N}}\bm{1}_{N\times
1}^{{\mkern-1.5mu\mathsf{T}}}{\mathbf{P}}_{\\!\\!\gamma}^{t}{\delta}_{i}\right)^{\alpha},$
(8)

where $\bm{1}_{N\times 1}$ is the all-ones vector, recall that
$\delta_{i}\in\\{0,1\\}^{N\times 1}$ is the vector returned by the Kronecker
delta function at $i\in\mathcal{V}$, and

$\mathbf{P}_{\gamma}=\gamma\mathbf{P}_{\textnormal{out}}+(1-\gamma)\mathbf{P}_{\textnormal{in}},$
(9)

where
$\mathbf{P}_{\textnormal{in}}=\mathbf{D}_{{\textnormal{in}}}^{-1}\mathbf{W}^{{\mkern-1.5mu\mathsf{T}}}$
and
$\mathbf{P}_{\textnormal{out}}=\mathbf{D}_{{\textnormal{out}}}^{-1}\mathbf{W}$
(recall that
$\mathbf{D}_{{\textnormal{out}}}=\mathbf{D}_{d_{\textnormal{out}}}$ and
$\mathbf{D}_{{\textnormal{in}}}=\mathbf{D}_{d_{\textnormal{in}}}$).

The three parameters have an easy-to-see role, and their use is optional as
one can set them to values that neutralize their effect (i.e. $t=1$,
$\alpha=1$, $\gamma=0.5$). The random walk iteration parameter $t$ controls
the random walk diffusion time, $\gamma$ controls the mixing between
$\mathbf{P}_{\textnormal{out}}$ (forward information) and
$\mathbf{P}_{\textnormal{in}}$ (backward information), and $\alpha$ controls
the re-weighting of the vertex measure. Plugging $\nu_{(t,\gamma)}^{\alpha}$
to Eq. 4 yields the following expression for the parametrized random walk
$\mathbf{P}_{t,\gamma}^{(\alpha)}=\mathbf{P}_{(\nu_{(t,\gamma)}^{\alpha})}$
and hence the P-RWDK $\mathbf{K}_{t_{d},\nu_{(t,\gamma)}^{\alpha}}$.

###  5.2 Determining the appropriate diffusion time in the unsupervised
setting

For a known number of clusters $k$, we aim at determining the best diffusion
time $t_{d}$ for our P-RWDKC algorithm. As we work in an unsupervised setting,
and hence lacking ground truth, determining the right diffusion time is
challenging (Shan & Daubechies, 2022). Our main insights to deal with this
matter are derived from the concepts of diffusion geometry (Coifman & Lafon,
2006), the theory of nearly uncoupled Markov chains (Tifenbach, 2011; Sharpe &
Wales, 2021) and the metastability of Markov chains (Landim & Xu, 2015).

As stated in (Coifman & Lafon, 2006), assuming that the graph has clusters
and/or a multi-scale structure, the random walk diffusion reveals clusters at
key moments of the diffusion at the course of time. The emergence of clusters
can be understood through the prism of metastability theory: for a given
cluster structure, there is a critical time-scale when an irreducible random
walk becomes nearly reducible, diffusing only inside clusters. Typically, this
means to observe approximately that, at that diffusion time, our operator
enjoys high intra-cluster compactness and high inter-cluster separation.
Cluster validity indexes can be employed for measuring the ratio of the former
to the latter of the quantities (José-García & Gómez-Flores, 2021). Here, we
propose the use of the Calinski–Harabasz criterion (CH; also known as Variance
ratio criterion) (Caliński & Harabasz, 1974), which computes the ratio between
the distance of the cluster centroids to the global centroid, and the distance
of the datapoints of each cluster to its cluster centroid.

To compute this variance criterion, we propose two types of distances between
data points, either the Euclidean distance, or a distance between their
embedded representations, which in the context of this work means to consider
the rows of a graph operator. Given the set of $N$ datapoints
$\bm{X}=\\{\bm{x}_{1},...,\bm{x}_{N}\\}$ partitioned into $k$ clusters,
denoted by $\bm{V}=\\{V_{j}\\}_{j=1}^{k}$, we denote by
$\bm{\mu}_{j}=\frac{1}{|V_{j}|}\sum_{\bm{x}_{i}\in V_{j}}\bm{x}_{i}$ the
centroid of cluster $j$, by
$\bm{\mu}=\frac{1}{N}\sum_{x_{i}\in\bm{X}}\bm{x}_{i}$ the centroid of
$\bm{X}$, and by $d(\bm{x}_{i},\bm{x}_{j})$ the pairwise distance used for two
datapoints $\bm{x}_{i},\bm{x}_{j}\in\bm{X}$. Thus, the CH criterion endowed
with a given distance $d$ is defined by:

$\textnormal{CH}\big{(}\bm{X},\bm{V}\big{)}=\frac{N-k}{k-1}\frac{\sum_{j=1}^{k}|V_{j}|\,d(\bm{\mu}_{j},\bm{\mu})}{\sum_{j=1}^{k}\sum_{\bm{x}_{i}\in
V_{j}}d(\bm{x}_{i},\bm{\mu}_{j})}.$

In the standard case where we consider multidimensional data vectors,
$\bm{X}=\\{\bm{x}_{1},...,\bm{x}_{N}\\}\in\mathbb{R}^{d}$, we can directly use
the Euclidean distance
$d(\bm{x}_{i},\bm{x}_{j})=\|\bm{x}_{i}-\bm{x}_{j}\|^{2}$ . However, when we
process a graph which is only represented by its adjacency matrix, the CH
criterion is no longer directly applicable. For this reason, as part of our
framework, we propose to extend the CH criterion for this setting, based on
the Kullback–Leibler divergence as a distance Van Erven & Harremos (2014).

######  Definition 5.1.

Probability Density-based Calinski-Harabasz (DCH) criterion. Let
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a digraph with cardinality
$|\mathcal{V}|=N$. Let $\mathbf{P}$ be the transition matrix of a random walk
on $\mathcal{G}$. The digraph $\mathcal{G}$ is partitioned into $k$ clusters
denoted by $\bm{V}=\\{V_{j}\\}_{j=1}^{k}$. Let
$\bm{p}(i,*)=\delta_{i}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{P}$ be the
conditional probability vector given the vertex $i$ that we consider as the
representation of the vertex $i$. We define as
$\bm{X}=\big{\\{}\bm{p}(1,*),...,\bm{p}(N,*)\big{\\}}\in\mathbb{R}^{N}$ the
set of conditional probability vectors associated with the vertices of
$\mathcal{G}$. Let us denote by
$\bm{\mu}_{j}=\frac{1}{|V_{j}|}\sum_{\bm{p}(i,*)\in V_{j}}\bm{p}(i,*)$ the
centroid of cluster $j$, and
$\bm{\mu}=\frac{1}{N}\sum_{\bm{p}(i,*)\in\bm{X}}\bm{p}(i,*)$ the centroid of
$\bm{X}$. The Kullback–Leibler divergence between two discrete probability
distributions, $\bm{p}$ and $\bm{q}$, is defined as
$\mathcal{D}_{\textnormal{KL}}(\bm{p},\bm{q})=\sum_{y}p(y)\log\frac{p(y)}{q(y)},\textnormal{s.t.}\,q(y)\neq
0$. Given a set of datapoints $\bm{X}$ and a partition $\bm{V}$, the
probability density-based Calinski-Harabasz (DCH) criterion endowed with the
Kullback–Leibler divergence $\mathcal{D}_{\textnormal{KL}}(\bm{p},\bm{q})$ is
defined by:

$\textnormal{{DCH}}\big{(}\bm{X},\bm{V}\big{)}=\frac{N-k}{k-1}\frac{\sum_{j=1}^{k}|V_{j}|\,\mathcal{D}_{\textnormal{KL}}(\bm{\mu}_{j},\bm{\mu})}{\sum_{j=1}^{k}\sum_{\bm{p}(i,*)\in
V_{j}}\mathcal{D}_{\textnormal{KL}}(\bm{p}(i,*),\bm{\mu}_{j})}.$

Algorithm 2 Estimating the random walk diffusion time ${t^{*}}$

Input: Dataset $\bm{X}$; $\mathbf{W}$: adjacency matrix; $k$: number of
clusters; ${J}$: max number of iterations
Output: ${t^{*}}$: the estimated diffusion time that best reveals $k$ clusters

1: Set $\nu=\bm{1}$ for the vertex measure

2: Compute the parametrized random walk operator $\mathbf{P}_{(\nu)}$, see Eq.
4

3: for $j=0$ to ${J}$ do

4: Apply $k$-means on $\mathbf{P}_{(\nu)}^{2^{j}}$

5: Obtain the $k$-partition $\bm{V}_{j}=\\{V_{q,j}\\}_{q=1}^{k}$ of the graph
vertices based on the clustering result of Step 4.

6: end for

7: Select $j^{\star}$=
$\underset{{j\in\\{0,...,J\\}}}{\textnormal{argmax}}\,\textnormal{CH}\big{(}\bm{X},\bm{V}_{j}\big{)}$.

8: return ${t^{*}}=2^{j^{\star}}$

We thus estimate the diffusion time in practice by evaluating the CH or DCH
criterion for the partitions associated with the dyadic powers of the
parametric random walk $\mathbf{P}_{(\nu)}^{2^{j}}$ with $\nu=\bm{1}$ (the
uniform measure) and $j\in\\{0,...,J\\}$ with $J=15$. Since the purpose is to
estimate the diffusion time, taking the vertex measure equal to a uniform
measure is sufficient for this situation. We summarize this procedure in Alg.
2.

##  6 Experiments

###  6.1 Setup and competitors

In this section, we demonstrate the effectiveness of our approach both on
digraphs obtained from high-dimensional data through graph construction
procedures, and real-world graphs. When dealing with high-dimensional data, we
use the $K$-nearest neighbor ($K$-NN) graph construction with
$K=\lfloor\log(N)\rfloor$, which produces (relatively) sparse and non-strongly
connected digraphs. A resulting $K$-NN graph is unweighted, directed, and
represented by its non-symmetric adjacency matrix
$\mathbf{W}=\\{w_{ij}\\}_{i,j=1}^{N}$, with entries
$w_{ij}=\mathds{1}\left\\{{\|x_{i}-x_{j}\|^{2}}\leq{\text{dist}_{K}(x_{i})}\right\\}$.
In the latter, $x_{i},x_{j}\in\mathbb{R}^{d}$ stand for the original
coordinates of the datapoints corresponding to the vertices $i$ and $j$,
$\text{dist}_{K}(x)$ is the Euclidean distance between $x$ and its $K$-th
nearest neighbor, and $\mathds{1}\\{\,\cdot\,\\}\in\\{0,1\\}$ is the indicator
function that evaluates the truth of the input condition.

Our clustering method, denoted by $\texttt{P-RWDKC}(\alpha,\gamma,t,t_{d})$,
is endowed with the parameters $\alpha\geq 0$, $\gamma\in[0,1]$, and
$t,t_{d}\geq 0$. The search grid used for each parameter is:
$\alpha\in\\{0,0.1,...,1\\}$, $t\in\\{0,1,...,100\\}$, and
$\gamma\in\\{0,0.1,...,1\\}$. Finally, the diffusion time parameter
$t_{d}\in\\{2^{0},...,2^{J}\\}$, with $J=15$, is estimated using Alg. 2. For
each method, we apply $k$-means clustering over the obtained embeddings (we
report the best score out of $100$ restarts). We select for each method the
optimal parameter values obtained through cross-validation over a grid search,
yielding the closest partition to the ground truth. The obtained partitions
are evaluated by the normalized mutual information (NMI) Strehl & Ghosh
(2002), which is a supervised cluster evaluation index that makes use of the
ground truth labels of the data.

In the experiments we compare with the following methods:
$\bullet$ $\texttt{DSC}\\!+\\!(\gamma)$ (Zhou et al., 2005; Palmer & Zheng,
2020) spectral clustering on strongly connected digraphs. To extend this
method to the graphs used in our experiments, we employ the teleporting random
walk (Page et al., 1999) endowed with the parameter $\gamma\in[0,1)$.
$\bullet$ $\texttt{DI-SIM}_{\texttt{L}}(\tau)$ and $\texttt{DI-
SIM}_{\texttt{R}}(\tau)$ (Rohe et al., 2016) are two variants that are based
on the left and the right singular vectors, respectively, of a given
regularized and normalized operator whose regularization is denoted by the
parameter $\tau\geq 0$. We use cross-validation to search the optimal
parameter with a grid search over $\tau\in\\{1,2,...,20\\}$.
$\bullet$ $\texttt{SC-SYM}_{1}$ and $\texttt{SC-SYM}_{2}$, two variants of the
classical spectral clustering (Von Luxburg, 2007) based on the unnormalized
and the normalized graph Laplacians obtained from the symmetrization of the
adjacency matrix $\mathbf{W}$.
$\bullet$ $\texttt{PIC}(t_{d})$ (Lin & Cohen, 2010) is the clustering approach
based on the power iteration of the random walk operator. We use the random
walk operator derived from the original adjacency matrix $\mathbf{W}$ of the
graph. The diffusion time $t_{d}$ is estimated using Alg. 2.
$\bullet$ RSC($\tau$) (Qin & Rohe, 2013; Zhang & Rohe, 2018) is the
regularized spectral clustering proposed to deal with sparse graphs
parametrized by $\tau\geq 0$. We use cross-validation with a grid search over
$\tau\in\\{1,2,...,20\\}$ to tune this parameter. The method applies only to
undirected graphs. To extend its use to digraphs, we need to symmetrize the
adjacency matrix of the original digraph.

(a) Ground truth I

(b) $6$ clusters result

(c) Ground truth II

(d) $2$ clusters result

Figure 1: Comparison between the ground truth classes at different scales and
the result of P-RWDKC on a synthetic toy-case. (a),(b) Small-scale: ground
truth with $6$ clusters and the P-RWDKC’s result. (c),(d) Large-scale: ground
truth with $2$ clusters and the finding of P-RWDKC.

###  6.2 Multi-scale synthetic Gaussians

In this section, we test the effectiveness of our approach at revealing
clusters at different scales thanks to the accurate estimation of the
diffusion time. We generate one instance of a point cloud
$X=\\{x_{i}\\}_{i=1}^{N}$, $\forall x_{i}\in\mathbb{R}^{2}$, of $N=300$ data
points drawn independently from the following mixture of six Gaussian
distributions
$\sum_{i=1}^{6}\alpha_{i}\mathcal{N}(\mu_{i},\sigma_{i}^{2}\bm{I})$ with
weights $\alpha_{i}\big{(}\sum_{i}\alpha_{i}=1\big{)}$. Specifically,
$\sigma_{i}=0.5,\alpha_{i}=1/6,\,\forall i$ and
$\mu_{1}=(-3,-2),\mu_{2}=(0,-2),\mu_{3}=(-1,1),\mu_{4}=(4,-2),\mu_{5}=(7,-2),\mu_{6}=(5,1)$.
The resulting data exhibit a multi-scale structure, as they can be seen as
either containing the $6$ original Gaussian clusters, or as having $2$
clusters where each of them is made up of $3$ smaller clusters. Figs. 1(a) and
1(c) show the ground truth data classes for the two scales. Figs. 1(b) and
1(d), show the respective clustering results obtained by P-RWDKC with a vertex
measure $\nu=\bm{1}$ and estimated diffusion times $t_{d_{1}}=64$ and
$t_{d_{2}}=128$. The clustering result is quite consistent compared to the
ground truth of each case, and hence providing evidence that our framework is
able at revealing clusters in a multi-scale fashion.

###  6.3 Real world data

In this section, we show that P-RWDKC is consistent with existing methods in
clustering high-dimensional data and real networks, and in nearly all tested
cases it outperforms them.

##### High-dimensional data benchmarks.

Here we report the results of experiments we conducted to evaluate the
performance of the proposed P-RWDKC method on $11$ benchmark datasets from the
UCI repository (Dheeru & Karra Taniskidou, 2017). We compare against
$\texttt{DSC}\\!+\\!$, $\texttt{SC-SYM}_{1}$ and $\texttt{SC-SYM}_{2}$, DI-
SIM, and PIC.

Tab. 1 summarizes the comparative results based on NMI. In nearly all cases,
the proposed P-RWDKC outperforms significantly the other methods on average.
Our approach performs better than SC-SYM1 and SC-SYM2. This allows us to state
that P-RWDKC, associated with the suitable vertex measure from Eq. 8 brings
indeed real added value to the clustering problem. Furthermore, P-RWDKC
outperforms PIC. Consequently, the RWDK operator produces better graph
embeddings than the original random walk operator defined the directed $K$-NN
graphs because the parametrized random walk is irreducible compared to the
random walk used in PIC that is not and thanks to the digraph information
encoded into the vertex measure.

Table 1: Clustering performance (NMI) on UCI datasets with optimal parameters
in brackets.

Dataset $N$ $d$ $k$ $\texttt{SC-SYM}_{1}$ $\texttt{SC-SYM}_{2}$ $\texttt{DI-
SIM}_{\texttt{L}}(\tau)$ $\texttt{DI-SIM}_{\texttt{R}}(\tau)$
$\texttt{DSC}\\!+\\!(\gamma)$ $\texttt{PIC}(t_{d})$
$\texttt{P-RWDKC}(\alpha,\gamma,t,t_{d})$ Iris 150 4 3 80.58 80.58 74.98 (1)
66.57 (1) 68.63 (0.80) 78.32 (32) 90.11 (0.4,1,49,32) Glass 214 9 6 38.59
38.92 38.95 (1) 36.41 (1) 39.72 (0.80) 42.79 (128) 44.39 (0.9,0.6,1,256) Wine
178 13 3 86.33 86.33 83.66 (1) 85.62 (1) 91.09 (0.80) 86.33 (4) 86.50
(0.7,0,21,32) WBDC 569 30 2 67.73 69.47 68.54 (2) 53.43 (1) 61.12 (0.10) 64.77
(8) 73.90 (0.9,0.3,100,2) Control Chart 600 60 6 81.17 81.17 82.94 (1) 77.72
(1) 79.45 (0.90) 82.79 (32) 87.49 (0.8,1,75,32) Parkinson 185 22 2 21.96 19.13
28.89 (1) 27.36 (13) 25.82 (0.95) 28.89 (2) 36.08 (1,0.3,52,2) Vertebral 310 6
3 39.26 39.26 52.06 (2) 41.76 (2) 56.63 (0.80) 49.13 (8) 62.40 (0.3,0.8,48,4)
Breast Tissue 106 9 6 54.03 54.43 54.04 (2) 49.33 (2) 51.64 (0.20) 54.18 (32)
60.43 (0.5,1,39,16) Seeds 210 7 3 73.90 73.90 76.29 (1) 73.06 (1) 74.80 (0.80)
70.79 (32) 78.95 (1.0,0.4,84,8) Image Seg. 2310 19 7 67.06 67.41 67.42 (1)
64.77 (1) 31.83 (0.99) 69.58 (214) 72.19 (0.8,1,80,256) Yeast 1484 8 10 30.58
31.11 31.37 (2) 28.89 (1) 27.50 (0.90) 32.62 (16) 34.05 (0.7,0.6,51,16)
Average – – – 58.29 58.34 59.92 54.77 56.37 60.01 66.04

Table 2: Clustering performance (NMI) on real-world datasets with optimal
parameters in brackets.

Dataset $N$ $k$ $\texttt{SC-SYM}_{1}$ RSC$(\tau)$ PIC$(t_{d})$
$\texttt{P-RWDKC}(\alpha,\gamma,t,t_{d})$ Political blogs (Adamic & Glance,
2005) 1222 2 1.74 73.25 (0.2) 57.37 (8) 75.53 (0.2,1,17,32) Cora (Bojchevski &
Günnemann, 2017) 2485 7 17.10 36.07 (2.5) 7.86 (4) 52.03 (0,0.9,72,8) College
Football (Girvan & Newman, 2002) 115 12 29.97 30.79 (2.9) 29.98 (2) 31.30
(0.7,0,17,2) Karate club (Zachary, 1977) 34 2 73.24 83.72 (1.7) 83.72 (1)
83.72 (0,0,1,1) Average – – 30.51 55.96 49.23 60.65

##### Real-world graph benchmarks.

To further validate the efficiency of P-RWDKC, we also evaluated our framework
on several real-world graphs. Tab.2 summarizes the comparative results
according to the NMI index. Due to space limitation, we include the
experiments and explanations in Appendix A.2.

##  7 Conclusion

We have proposed the _parametrized random walk diffusion kernel clustering_
(P-RWDKC) that applies to both directed and undirected graphs. First, we
introduced the parametrized random walk (P-RW) operator. We then show that the
diffusion distance, is a Mahalanobis distance involving a special kernel
matrix, called random walk diffusion kernel (RWDK) which is simply a power of
the transition matrix (normalized by the vertex degrees). From this, we show
that RWDK is an alternative to the eigendecomposition step of the spectral
clustering pipeline. We propose to extend the diffusion geometry framework to
digraphs by combining RWDK and P-RW. Our proposal comes with a clustering
algorithm stemming from our framework. Finally, we demonstrated empirically
with extensive experiments on several datasets that P-RWDKC outperforms
existing approaches for digraphs.

##  8 Acknowledgment

Harry Sevi and Argyris Kalogeratos were funded by the Industrial Analytics and
Machine Learning (IdAML) Chair hosted at ENS Paris-Saclay, University Paris-
Saclay. Matthieu Jonckheere was funded by the International Centre for
Mathematics and Computer Science (CIMI) in Toulouse.

## References

  * Adamic & Glance (2005) Lada A Adamic and Natalie Glance.  The political blogosphere and the 2004 us election: divided they blog.  In _Proceedings of the International Workshop on Link discovery_ , pp. 36–43, 2005.
  * Azran & Ghahramani (2006) Arik Azran and Zoubin Ghahramani.  A new approach to data driven clustering.  In _Proceedings of the International Conference on Machine Learning_ , pp. 57–64, 2006.
  * Boedihardjo et al. (2021) March Boedihardjo, Shaofeng Deng, and Thomas Strohmer.  A performance guarantee for spectral clustering.  _SIAM Journal on Mathematics of Data Science_ , 3(1):369–387, 2021.
  * Bojchevski & Günnemann (2017) Aleksandar Bojchevski and Stephan Günnemann.  Deep gaussian embedding of graphs: Unsupervised inductive learning via ranking.  _arXiv preprint arXiv:1707.03815_ , 2017.
  * Boyd et al. (2021) Zachary M Boyd, Nicolas Fraiman, Jeremy Marzuola, Peter J Mucha, Braxton Osting, and Jonathan Weare.  A metric on directed graphs and markov chains based on hitting probabilities.  _SIAM Journal on Mathematics of Data Science_ , 3(2):467–493, 2021.
  * Brémaud (2013) Pierre Brémaud.  _Markov chains: Gibbs fields, Monte Carlo simulation, and queues_ , volume 31.  Springer Science & Business Media, 2013.
  * Caliński & Harabasz (1974) Tadeusz Caliński and Jerzy Harabasz.  A dendrite method for cluster analysis.  _Communications in Statistics-theory and Methods_ , 3(1):1–27, 1974.
  * Chung (2005) Fan Chung.  Laplacians and the Cheeger inequality for directed graphs.  _Annals of Combinatorics_ , 9(1):1–19, 2005.
  * Chung & Graham (1997) Fan RK Chung and Fan Chung Graham.  _Spectral graph theory_.  Number 92. American Mathematical Soc., 1997.
  * Coifman & Lafon (2006) Ronald R Coifman and Stéphane Lafon.  Diffusion maps.  _Applied and computational harmonic analysis_ , 21(1):5–30, 2006.
  * Coifman et al. (2005) Ronald R Coifman, Stephane Lafon, Ann B Lee, Mauro Maggioni, Boaz Nadler, Frederick Warner, and Steven W Zucker.  Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps.  _Proceedings of the National Academy of Sciences_ , 102(21):7426–7431, 2005.
  * Dall’Amico et al. (2021) Lorenzo Dall’Amico, Romain Couillet, and Nicolas Tremblay.  A unified framework for spectral clustering in sparse graphs.  _J. Mach. Learn. Res._ , 22:217–1, 2021.
  * Dall’Amico et al. (2020) Lorenzo Dall’Amico, Romain Couillet, and Nicolas Tremblay.  Optimal laplacian regularization for sparse spectral community detection.  In _IEEE International Conference on Acoustics, Speech and Signal Processing_ , pp. 3237–3241, 2020.
  * Dheeru & Karra Taniskidou (2017) D Dheeru and E Karra Taniskidou.  UCI repository of learnin learning databases.  _University of California, Irvine, School of Information and Computer Sciences_ , 2017.
  * Girvan & Newman (2002) Michelle Girvan and Mark EJ Newman.  Community structure in social and biological networks.  _Proceedings of the National Academy of Sciences_ , 99(12):7821–7826, 2002.
  * José-García & Gómez-Flores (2021) Adán José-García and Wilfrido Gómez-Flores.  A survey of cluster validity indices for automatic data clustering using differential evolution.  In _Proceedings of the Genetic and Evolutionary Computation Conference_ , pp. 314–322, 2021.
  * Kiselev et al. (2017) Vladimir Yu Kiselev, Kristina Kirschner, Michael T Schaub, Tallulah Andrews, Andrew Yiu, Tamir Chandra, Kedar N Natarajan, Wolf Reik, Mauricio Barahona, Anthony R Green, et al.  Sc3: consensus clustering of single-cell rna-seq data.  _Nature methods_ , 14(5):483–486, 2017.
  * Landim & Xu (2015) Claudio Landim and Tiecheng Xu.  Metastability of finite state markov chains: a recursive procedure to identify slow variables for model reduction.  _arXiv preprint arXiv:1512.06597_ , 2015.
  * Levin & Peres (2017) David A Levin and Yuval Peres.  _Markov chains and mixing times_ , volume 107.  American Mathematical Society, 2017.
  * Lin & Cohen (2010) Frank Lin and William W Cohen.  Power iteration clustering.  In _Proceedings of the International Conference on Machine Learning_ , 2010.
  * Maggioni & Murphy (2019) Mauro Maggioni and James M Murphy.  Learning by unsupervised nonlinear diffusion.  _J. Mach. Learn. Res._ , 20(160):1–56, 2019.
  * McFee & Ellis (2014) Brian McFee and Dan Ellis.  Analyzing song structure with spectral clustering.  In _ISMIR_ , pp. 405–410. Citeseer, 2014.
  * Meilă & Shi (2001) Marina Meilă and Jianbo Shi.  A random walks view of spectral segmentation.  In _International Workshop on Artificial Intelligence and Statistics_ , pp. 203–208, 2001.
  * Murphy & Polk (2022) James M Murphy and Sam L Polk.  A multiscale environment for learning by diffusion.  _Applied and Computational Harmonic Analysis_ , 57:58–100, 2022.
  * Nadler & Galun (2006) Boaz Nadler and Meirav Galun.  Fundamental limitations of spectral clustering.  _Advances in Neural Information Processing Systems_ , 2006.
  * Nadler et al. (2006) Boaz Nadler, Stephane Lafon, Ronald R. Coifman, and Ioannis G Kevrekidis.  Diffusion maps, spectral clustering and eigenfunctions of fokker-planck operators.  In _Advances in Neural Information Processing Systems_ , 2006.
  * Ng et al. (2002) Andrew Y Ng, Michael I Jordan, and Yair Weiss.  On spectral clustering: Analysis and an algorithm.  In _Advances in Neural Information Processing Systems_ , 2002.
  * Page et al. (1999) Lawrence Page, Sergey Brin, Rajeev Motwani, and Terry Winograd.  The Pagerank citation ranking: Bringing order to the web.  Technical report, Stanford InfoLab, 1999.
  * Palmer & Zheng (2020) William R Palmer and Tian Zheng.  Spectral clustering for directed networks.  In _International Conference on Complex Networks and Their Applications_ , pp. 87–99. Springer, 2020.
  * Peng et al. (2015) Richard Peng, He Sun, and Luca Zanetti.  Partitioning well-clustered graphs: Spectral clustering works!  In _Conference on Learning Theory_ , pp. 1423–1455, 2015.
  * Pons & Latapy (2005) Pascal Pons and Matthieu Latapy.  Computing communities in large networks using random walks.  In _International symposium on computer and information sciences_ , pp. 284–293. Springer, 2005.
  * Qin & Rohe (2013) Tai Qin and Karl Rohe.  Regularized spectral clustering under the degree-corrected stochastic blockmodel.  _preprint arXiv:1309.4111_ , 2013.
  * Rohe et al. (2011) Karl Rohe, Sourav Chatterjee, and Bin Yu.  Spectral clustering and the high-dimensional stochastic blockmodel.  _The Annals of Statistics_ , 39(4):1878–1915, 2011.
  * Rohe et al. (2016) Karl Rohe, Tai Qin, and Bin Yu.  Co-clustering directed graphs to discover asymmetries and directional communities.  _Proceedings of the National Academy of Sciences_ , 113(45):12679–12684, 2016.
  * Sevi et al. (2022) Harry Sevi, Matthieu Jonckheere, and Argyris Kalogeratos.  Generalized spectral clustering for directed and undirected graphs.  _arXiv preprint arXiv:2203.03221_ , 2022.
  * Shan & Daubechies (2022) Shan Shan and Ingrid Daubechies.  Diffusion maps: Using the semigroup property for parameter tuning.  _arXiv preprint arXiv:2203.02867_ , 2022.
  * Sharpe & Wales (2021) Daniel J Sharpe and David J Wales.  Nearly reducible finite markov chains: Theory and algorithms.  _The Journal of Chemical Physics_ , 155(14):140901, 2021.
  * Shi & Malik (2000) Jianbo Shi and Jitendra Malik.  Normalized cuts and image segmentation.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 22(8):888–905, 2000.
  * Strehl & Ghosh (2002) Alexander Strehl and Joydeep Ghosh.  Cluster ensembles—a knowledge reuse framework for combining multiple partitions.  _Journal of Machine Learning Research_ , 3(Dec):583–617, 2002.
  * Tifenbach (2011) Ryan M Tifenbach.  _A combinatorial approach to nearly uncoupled Markov chains_.  PhD thesis, National University of Ireland Maynooth, 2011.
  * Tishby & Slonim (2000) Naftali Tishby and Noam Slonim.  Data clustering by Markovian relaxation and the information bottleneck method.  _Advances in Neural Information Processing Systems_ , 2000.
  * Tishby et al. (2000) Naftali Tishby, Fernando C Pereira, and William Bialek.  The information bottleneck method.  _arXiv preprint physics/0004057_ , 2000.
  * Tremblay et al. (2016) Nicolas Tremblay, Gilles Puy, Rémi Gribonval, and Pierre Vandergheynst.  Compressive spectral clustering.  In _Proceedings of the International Conference on Machine Learning_ , pp. 1002–1011, 2016.
  * Van Erven & Harremos (2014) Tim Van Erven and Peter Harremos.  Rényi divergence and kullback-leibler divergence.  _IEEE Transactions on Information Theory_ , 60(7):3797–3820, 2014.
  * Von Luxburg (2007) Ulrike Von Luxburg.  A tutorial on spectral clustering.  _Statistics and Computing_ , 17(4):395–416, 2007\.
  * Zachary (1977) Wayne W Zachary.  An information flow model for conflict and fission in small groups.  _Journal of anthropological research_ , 33(4):452–473, 1977.
  * Zhang & Rohe (2018) Yilin Zhang and Karl Rohe.  Understanding regularized spectral clustering via graph conductance.  _preprint arXiv:1806.01468_ , 2018.
  * Zhou et al. (2005) Dengyong Zhou, Jiayuan Huang, and Bernhard Schölkopf.  Learning from labeled and unlabeled data on a directed graph.  In _Proceedings of the International Conference on Machine Learning_ , pp. 1036–1043, 2005.

##  Appendix A Appendix

###  A.1 Proofs

####  A.1.1 Proof of Proposition 4.1

######  Proposition A.1.

Let $\mathcal{X}$ be a random walk on an undirected graph $\mathcal{G}$ with
transition matrix $\mathbf{P}$ and ergodic distribution $\pi$. The transition
matrix $\mathbf{P}$ admits the following eigendcomposition
$\mathbf{P}=\bm{\Phi}\mathbf{D}_{\lambda}\bm{\Psi}^{{\mkern-1.5mu\mathsf{T}}}$.
The diffusion distance between vertices $i$ and $j$, $d_{t}^{2}(i,j)$ at a
given time $t\in\mathbb{N}$ can be written as the following Mahalanobis
distance

$d_{t}^{2}(i,j)=(\delta_{i}-\delta_{j})^{{\mkern-1.5mu\mathsf{T}}}\mathbf{K}_{t}(\delta_{i}-\delta_{j}),$

where the similarity positive definite kernel matrix $\mathbf{K}_{t}$ is
defined by

$\mathbf{K}_{t}=\mathbf{P}^{2t}\mathbf{D}_{d}^{-1}.$ (10)

###### Proof.

$\displaystyle d_{t}^{2}(i,j)$
$\displaystyle=\|p_{t}(i,*)-p_{t}(j,*)\|_{1/\pi}^{2},$
$\displaystyle=\|(\mathbf{P}^{t})^{{\mkern-1.5mu\mathsf{T}}}(\delta_{i}-\delta_{j})\|^{2}_{1/\pi},$
$\displaystyle d_{t}^{2}(i,j)$
$\displaystyle=(\delta_{i}-\delta_{j})^{{\mkern-1.5mu\mathsf{T}}}\mathbf{P}^{t}\mathbf{D}_{\pi}^{-1}(\mathbf{P}^{t})^{{\mkern-1.5mu\mathsf{T}}}(\delta_{i}-\delta_{j}).$

By setting
$\mathbf{K}_{t}=\mathbf{P}^{t}\mathbf{D}_{\pi}^{-1}(\mathbf{P}^{t})^{{\mkern-1.5mu\mathsf{T}}}$,
using the eigendecomposition of $\mathbf{P}$ and the fact that $\mathbf{P}$ is
self-adjoint in $\ell^{2}(\mathcal{V},\pi)$, we have

$\displaystyle\mathbf{K}_{t}$
$\displaystyle=\mathbf{P}^{t}\mathbf{D}_{\pi}^{-1}(\mathbf{P}^{t})^{{\mkern-1.5mu\mathsf{T}}}$
$\displaystyle=\bm{\Phi}\mathbf{D}_{\lambda}^{t}\bm{\Psi}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\pi}^{-1}\bm{\Psi}\mathbf{D}_{\lambda}^{t}\bm{\Phi}^{{\mkern-1.5mu\mathsf{T}}},$
$\displaystyle=\bm{\Phi}\mathbf{D}_{\lambda}^{2t}\bm{\Phi}^{{\mkern-1.5mu\mathsf{T}}}\>(\bm{\Psi}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\pi}^{-1}\bm{\Psi}=\mathbf{I}),$
$\displaystyle=\bm{\Phi}\mathbf{D}_{\lambda}^{2t}\bm{\Psi}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{d}^{-1}\>(\bm{\Phi}=\mathbf{D}_{d}^{-1}\bm{\Psi}),$
$\displaystyle\mathbf{K}_{t}$
$\displaystyle=\mathbf{P}^{2t}\mathbf{D}_{d}^{-1}.$

∎

####  A.1.2 Supplementary proofs

######  Proposition A.2.

$\mathbf{P}_{(\nu)}$ is a transition matrix and reversible.

###### Proof.

We have the following equality

$\mathbf{P}_{(\nu)}=\mathbf{(I+\mathbf{D}_{\nu}^{-1}\mathbf{D}_{\xi})^{-1}(P+\mathbf{D}_{\nu}^{-1}P^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu})}=\mathbf{(\mathbf{D}_{\nu}+\mathbf{D}_{\xi})^{-1}(\mathbf{D}_{\nu}P+P^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu})}.$

As a result, we need to show that
$\mathbf{(\mathbf{D}_{\nu}+\mathbf{D}_{\xi})^{-1}(\mathbf{D}_{\nu}P+P^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu})}$
is a transition matrix, i.e. we show that

$\sum_{j=1}^{|\mathcal{V}|}\mathbf{P}_{(\nu),ij}=1,\quad\forall
i\in\mathcal{V}$
$\displaystyle\sum_{j=1}^{|\mathcal{V}|}\mathbf{P}_{(\nu),ij}$
$\displaystyle=\sum_{j=1}^{|\mathcal{V}|}\bigg{(}\mathbf{(\mathbf{D}_{\nu}+\mathbf{D}_{\xi})^{-1}(\mathbf{D}_{\nu}P+P^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu})}\bigg{)}_{ij}$
$\displaystyle=\sum_{j=1}^{|\mathcal{V}|}\bigg{(}\sum_{k=1}^{|\mathcal{V}|}(\mathbf{D}_{\nu}+\mathbf{D}_{\xi})_{ik}^{-1}(\mathbf{D}_{\nu}\mathbf{P}+\mathbf{P}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu})_{kj}\bigg{)}$
$\displaystyle\sum_{j=1}^{|\mathcal{V}|}\mathbf{P}_{(\nu),ij}$
$\displaystyle=(\mathbf{D}_{\nu}+\mathbf{D}_{\xi})_{ii}^{-1}\sum_{j=1}^{|\mathcal{V}|}(\mathbf{D}_{\nu}\mathbf{P}+\mathbf{P}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu})_{ij}$
(11)

$\displaystyle\sum_{j=1}^{|\mathcal{V}|}(\mathbf{D}_{\nu}\mathbf{P}+\mathbf{P}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu})_{ij}$
$\displaystyle=\sum_{j=1}^{|\mathcal{V}|}\sum_{k=1}^{|\mathcal{V}|}(\mathbf{D}_{\nu})_{ik}\mathbf{P}_{kj}+\sum_{j=1}^{|\mathcal{V}|}\sum_{k=1}^{|\mathcal{V}|}(\mathbf{P}^{{\mkern-1.5mu\mathsf{T}}})_{ik}(\mathbf{D}_{\nu})_{kj}$
$\displaystyle=\sum_{j=1}^{|\mathcal{V}|}(\mathbf{D}_{\nu})_{ii}\mathbf{P}_{ij}+\sum_{j=1}^{|\mathcal{V}|}(\mathbf{P}^{{\mkern-1.5mu\mathsf{T}}})_{ij}(\mathbf{D}_{\nu})_{jj}$
$\displaystyle=\sum_{j=1}^{|\mathcal{V}|}\nu(i){p}(i,j)+\sum_{j=1}^{|\mathcal{V}|}\nu(j)p(j,i)$
$\displaystyle\sum_{j=1}^{|\mathcal{V}|}(\mathbf{D}_{\nu}\mathbf{P}+\mathbf{P}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu})_{ij}$
$\displaystyle=\nu(i)+\xi(i),\quad\forall\,i\in\mathcal{V}$

Using Eq. 11, we finally obtain
$\sum_{j=1}^{|\mathcal{V}|}\mathbf{P}_{(\nu),ij}=1$. ∎

###  A.2 Experiments on real-world graph benchmarks

Most of the real-world graphs are sparse with heterogeneous degrees with a
high variance of the degree. As a result, spectral clustering mostly fails for
these graphs Zhang & Rohe (2018); Rohe et al. (2011); Dall’Amico et al.
(2020). To avoid numerical issues due to the high variance degree, we proposed
a slight modification of the parametrized random walk operator (transition
matrix).

$\mathbf{P}_{(\nu)}=(\mathbf{I}+\mathbf{D}_{\xi/\tilde{\nu}})^{-1}\mathbf{D}_{\textnormal{out}}^{-1}(\mathbf{W}+\mathbf{D}_{\nu}^{-1}\mathbf{W}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{D}_{\nu}).$
(12)

with $\xi=\nu^{{\mkern-1.5mu\mathsf{T}}}\mathbf{W}$ and the vertex measure
$\tilde{\nu}=[\tilde{\nu}_{1},...,\tilde{\nu}_{N}]^{{\mkern-1.5mu\mathsf{T}}}$
defined by $\tilde{\nu}_{i}=\nu_{i}d_{i}^{\textnormal{out}},\forall
i\in\mathcal{V}$.

We report the results of experiments that evaluate the performance of the
proposed P-RWDKC method on $4$ real-world networks. Among these, $2$ are
directed (Political blogs, Cora) and $2$ are undirected (Karate Club, College
Football). For all networks, the number of clusters is considered known.
Moreover, when necessary, clustering is performed on the largest connected
component of the graph. We compare against $\texttt{SC-SYM}_{1}$, RSC, and
PIC.

Tab. 2 summarizes the comparative results according to the NMI index. In all
cases, the proposed P-RWDKC outperforms the other methods, and when looking at
the overall average performance the difference is significant. Our approach
performs significantly better than $\texttt{SC-SYM}_{1}$. This is caused by a
mixed effect of the symmetrization of the digraph as well as the hard
truncation phenomenon discussed in Sec. 4.2. RSC is competitive against our
method. Nevertheless, the adjacency matrix involved in RSC is dense and
clearly modified, which can have a deleterious impact compared to our method
(e.g. see the results on Cora).

###  A.3 Alternative for the design of the vertex measure

In Sec. 5.1 we presented a design for the vertex measure to be used by
P-RWDKC. Here we discuss an alternative design. As before, the vertex measure
can be parametrized by three parameters
$(t\\!\in\\!\mathbb{N},\,\gamma\\!\in\\![0,1],\,\alpha\\!\in\\!\mathbb{R})$
and is formally given by:

$\nu_{(t,\gamma)}^{\alpha}(i)=\left({\textstyle\frac{1}{N}}\bm{1}_{N\times
1}^{{\mkern-1.5mu\mathsf{T}}}{\mathbf{P}}_{\\!\\!\gamma}^{t}{\delta}_{i}\right)^{\alpha},$
(13)

where $\bm{1}_{N\times 1}$ is the all-ones vector ,
$\delta_{i}\in\\{0,1\\}^{N\times 1}$ is the vector output of the Kronecker
delta function at $i\in\mathcal{V}$, and

$\mathbf{P}_{\gamma}=\mathbf{D}_{\gamma}^{-1}\mathbf{W}_{\gamma},\quad\mathbf{D}_{\gamma}=\textnormal{diag}(\mathbf{W}_{\gamma}\bm{1}),$
(14)

where $\mathbf{W}_{\gamma}$ is defined by

$\mathbf{W}_{\gamma}=\gamma\mathbf{W}+(1-\gamma)\mathbf{W}^{{\mkern-1.5mu\mathsf{T}}}\,,\quad\gamma\in[0,1].$
(15)

The random walk iteration parameter $t$ controls the random walk diffusion,
$\gamma$ controls the influence between the original adjacency $\mathbf{W}$
(forward information) and its transpose
$\mathbf{W}^{{\mkern-1.5mu\mathsf{T}}}$ (backward information), and $\alpha$
controls the re-weighting of the vertex measure.

####  A.3.1 Vertex measure when $\gamma=1/2$ and $t\rightarrow\infty$

In the setting where, $\gamma=\frac{1}{2}$ and $t\rightarrow\infty$, we are
able to characterize explicitly the vertex measure. Indeed,
$\mathbf{W}_{\gamma}$ becomes symmetric, the associated random walk is thus
ergodic,
$\lim_{t\rightarrow\infty}\delta_{i}^{{\mkern-1.5mu\mathsf{T}}}\mathbf{P}_{1/2}^{t}=\pi_{1/2},\forall
i$ and the parametrized vertex measure defined in Eq. 13 becomes
$\nu_{(t,\gamma)}^{\alpha}(i)=\big{(}\pi_{{1}/{2}}(i)\big{)}^{\alpha}$ for any
vertex $i$. Consequently, in this setting, the vertex measure is only
parametrized by $\alpha$.

####  A.3.2 Additional Experiments

In this section, we demonstrate the effectiveness of our approach based on the
same setting described in Sec. 6, using the vertex measure in Sec. A.3.1.
Tabs. 3 and 5 summarize the comparative results based on NMI. Similarly to the
results from Secs. 6.3 and 6.3, we observe that the proposed P-RWDKC
outperforms significantly the other methods in nearly all cases. The P-RWDKC,
associated with the vertex measure defined in Sec. A.3.1 stays competitive
against P-RWDKC from Sec. 6.3 with less degree of freedom. The main takeway to
consider is that the construction of this kernel and the design of the vertex
measures are the key elements to creating relevant embedding for clustering
graphs.

Table 3: Clustering performance (NMI) on UCI datasets with optimal parameters
in brackets.

Dataset $N$ $d$ $k$ $\texttt{SC-SYM}_{1}$ $\texttt{SC-SYM}_{2}$ $\texttt{DI-
SIM}_{\texttt{L}}(\tau)$ $\texttt{DI-SIM}_{\texttt{R}}(\tau)$
$\texttt{DSC}\\!+\\!(\gamma)$ $\texttt{PIC}(t_{d})$
$\texttt{P-RWDK}(\alpha,t_{d})$ Iris 150 4 3 80.58 80.58 74.98 (1) 66.57 (1)
68.63 (0.80) 78.32 (32) 90.11 (0.6,4) Glass 214 9 6 38.59 38.92 38.95 (1)
36.41 (1) 39.72 (0.80) 42.79 (128) 42.89 (0.5,256) Wine 178 13 3 86.33 86.33
83.66 (1) 85.62 (1) 91.09 (0.80) 86.33 (4) 84.73 (0,32) WBDC 569 30 2 67.73
69.47 68.54 (2) 53.43 (1) 61.12 (0.10) 64.77 (8) 71.32 (1,2) Control Chart 600
60 6 81.17 81.17 82.94 (1) 77.72 (1) 79.45 (0.90) 82.79 (32) 84.22 (0.6,32)
Parkinson 185 22 2 21.96 19.13 28.89 (1) 27.36 (13) 25.82 (0.95) 28.89 (2)
36.08 (0.5,4) Vertebral 310 6 3 39.26 39.26 52.06 (2) 41.76 (2) 56.63 (0.80)
49.13 (8) 59.57 (0.1,4) Breast Tissue 106 9 6 54.03 54.43 54.04 (2) 49.33 (2)
51.64 (0.20) 54.18 (32) 55.29 (0.6,16) Seeds 210 7 3 73.90 73.90 74.89 (1)
73.06 (1) 74.80 (0.80) 70.79 (32) 75.52 (0,8) Image Seg. 2310 19 7 67.06 67.41
67.42 (1) 64.77 (1) 31.83 (0.99) 69.58 (214) 70.70 (0.5,64) Yeast 1484 8 10
30.58 31.11 31.37 (2) 28.89 (1) 27.50 (0.90) 32.62 (16) 33.83 (0.1,16) Average
– – – 58.29 58.34 59.92 54.77 56.37 60.01 63.92

To further evaluate the P-RWDKC framework, instead of using the ground truth
and cross-validation, here we choose the parameter values that optimize the
Calinski-Harabasz (CH). We first compute the set of candidate partitions, one
for each ($\alpha$, $t_{d}$) combination in the considered parameter grid.
Then we select as best the model with the highest CH and compute its NMI. We
operate in the same way for the methods that have parameters (i.e. all but
$\texttt{SC-SYM}_{1}$, $\texttt{SC-SYM}_{2}$, and PIC whose results are the
same as in Tabs. 3 and 5. The comparative results are shown in Tabs. 4 and 6.
Notice that P-RWDKC outperforms significantly the other methods in nearly all
cases (see also the average performance in the last row of the table.).
Compared to Tab. 3 here the NMI of P-RWDKC stays just a little lower. This
indicates that the unsupervised tuning of the model parameters offers
comparable graph partition quality to the previous case where we applied
cross-validation using the ground truth.

Table 4: Clustering performance (NMI) on UCI datasets with estimated
parameters (shown in brackets) according to the CH (or DCH) index.

Dataset $N$ $d$ $k$ $\texttt{SC-SYM}_{1}$ $\texttt{SC-SYM}_{2}$ $\texttt{DI-
SIM}_{\texttt{L}}(\tau)$ $\texttt{DI-SIM}_{\texttt{R}}(\tau)$
$\texttt{DSC}\\!+\\!(\gamma)$ $\texttt{PIC}(t_{d})$
$\texttt{P-RWDK}(\alpha,t_{d})$ Iris 150 4 3 80.58 80.58 74.98 (1) 66.57 (1)
68.63 (0.80) 78.32 (32) 90.11 (0.6,4) Glass 214 9 6 38.59 38.92 37.39 (2)
35.87 (1) 36.58 (0.80) 42.79 (128) 42.13 (0.3,256) Wine 178 13 3 86.33 86.33
83.66 (1) 85.62 (1) 91.09 (0.80) 86.33 (4) 84.73 (0.9,32) WBDC 569 30 2 67.73
69.47 68.54 (2) 53.43 (1) 61.12 (0.10) 64.77 (8) 70.24 (0.5,2) Control Chart
600 60 6 81.17 81.17 82.94 (2) 77.44 (1) 79.45 (0.90) 82.79 (32) 83.19
(0.2,32) Parkinson 185 22 2 21.96 19.13 28.89 (1) 27.36 (13) 22.97 (0.95)
28.89 (2) 36.08 (0.7,4) Vertebral 310 6 3 39.26 39.26 45.89 (2) 39.62 (1)
54.24 (0.80) 49.13 (8) 59.34 (0.3,4) Breast Tissue 106 9 6 54.03 54.43 54.04
(2) 49.27 (2) 51.64 (0.20) 54.18 (32) 54.90 (0.2,16) Seeds 210 7 3 73.90 73.90
74.89 (1) 73.06 (1) 74.80 (0.80) 70.79 (32) 74.89 (0.1,8) Image Seg. 2310 19 7
67.06 67.41 67.42 (1) 64.77 (1) 31.46 (0.99) 69.58 (2${}^{14})$ 68.79 (0.1,64)
Yeast 1484 8 10 30.58 31.11 31.22 (1) 28.89 (1) 27.47 (0.90) 32.62 (16) 33.60
(0,16) Average – – – 58.29 58.34 58.82 54.57 51.95 60.01 63.45

Table 5: Clustering performance (NMI) on real-world datasets with optimal
parameters in brackets.

Dataset $N$ $k$ $\texttt{SC-SYM}_{1}$ RSC$(\tau)$ PIC$(t_{d})$
$\texttt{P-RWDKC}(\alpha,t_{d})$ Political blogs (Adamic & Glance, 2005) 1222
2 1.74 73.25 (0.2) 57.37 (8) 74.02 (0.1,32) Cora (Bojchevski & Günnemann,
2017) 2485 7 17.10 36.07 (2.5) 7.86 (4) 49.71 (0.2,16) College Football
(Girvan & Newman, 2002) 115 12 29.97 30.79 (2.9) 29.98 (2) 30.79 (0.9,2)
Karate club (Zachary, 1977) 34 2 73.24 83.72 (1.7) 83.72 (1) 83.72 (0,1)

Table 6: Clustering performance (NMI) on real-world datasets with estimated
parameters (shown in brackets) according to the CH (or DCH) index.

Dataset $N$ $k$ $\texttt{SC-SYM}_{1}$ RSC$(\tau)$ PIC$(t_{d})$
$\texttt{P-RWDKC}(\alpha,t_{d})$ Political blogs (Adamic & Glance, 2005) 1222
2 1.74 73.25 (0.2) 57.37 (8) 73.93 (0,32) Cora (Bojchevski & Günnemann, 2017)
2485 7 17.10 30.65 (0.5) 7.86 (4) 49.71 (0.2,16) College Football (Girvan &
Newman, 2002) 115 12 29.97 29.84 (2.9) 29.98 (2) 30.30 (0.9,2) Karate club
(Zachary, 1977) 34 2 73.24 83.72 (1.7) 83.72 (1) 83.72 (0,1)


# Kick-motion Training with DQN in AI Soccer Environment
††thanks: This work was supported by the Institute for Information
communications Technology Promotion (IITP) grant funded by the Korean
government (MSIT) (No.2020-0-00440, Development of Artificial Intelligence
Technology that continuously improves itself as the situation changes in the
real world).

Bumgeun Park  Cho Chun Shik Graduate School of Mobility
Korea Advanced Institute of Science and Technology (KAIST)
Daejeon, South Korea
<EMAIL_ADDRESS>Jihui Lee  Division of Future Vehicle
Korea Advanced Institute of Science and Technology (KAIST)
Daejeon, South Korea
<EMAIL_ADDRESS>Taeyoung Kim  Cho Chun Shik Graduate School of Mobility
Korea Advanced Institute of Science and Technology (KAIST)
Daejeon, South Korea
<EMAIL_ADDRESS>Dongsoo Har  Cho Chun Shik Graduate School of Mobility
Korea Advanced Institute of Science and Technology (KAIST)
Daejeon, South Korea
<EMAIL_ADDRESS>
###### Abstract

This paper presents a technique to train a robot to perform kick-motion in AI
soccer by using reinforcement learning (RL). In RL, an agent interacts with an
environment and learns to choose an action in a state at each step. When
training RL algorithms, a problem called the curse of dimensionality (COD) can
occur if the dimension of the state is high and the number of training data is
low. The COD often causes degraded performance of RL models. In the situation
of the robot kicking the ball, as the ball approaches the robot, the robot
chooses the action based on the information obtained from the soccer field. In
order not to suffer COD, the training data, which are experiences in the case
of RL, should be collected evenly from all areas of the soccer field over
(theoretically infinite) time. In this paper, we attempt to use the relative
coordinate system (RCS) as the state for training kick-motion of robot agent,
instead of using the absolute coordinate system (ACS). Using the RCS
eliminates the necessity for the agent to know all the (state) information of
entire soccer field and reduces the dimension of the state that the agent
needs to know to perform kick-motion, and consequently alleviates COD. The
training based on the RCS is performed with the widely used Deep Q-network
(DQN) and tested in the AI Soccer environment implemented with Webots
simulation software.

###### Index Terms:

Reinforcement learning (RL), Deep Q-Network (DQN), AI soccer, Curse of
dimensionality (COD), Coordinate transformation matrix (CTM)

##  I Introduction

Reinforcement Learning (RL), which is one of the widely used machine learning
techniques, makes an agent maximize the accumulated rewards in given
environment [1]. To do so requires the use of experiences, e.g., experience
replay, during training. The agent selects an action based on the current
state and receives the reward from the environment. Through this interaction,
the agent can learn which action is more useful to achieve the goal. Recently,
a deep neural network (DNN) for nonlinear approximation enables the RL to
handle more complicated tasks such as playing a range of Atari games [2],
playing the game of GO [3, 4], controlling the robot arm [5, 6, 7], and
planning path for mobile robots [8, 9].

AI soccer environment, which is one of the environments for training multi-
agent RL algorithms, is introduced in the previous work [10]. The main target
of AI Soccer is winning the soccer game by controlling each robot. Each
robot’s movement consists of many small motions such as ”going to a target
spot”, ”blocking robot or ball”, ”positioning”, ”heading”, and ”kicking”.
Therefore, the entire AI soccer training can be divided into various tasks of
training motions. The state, action, and reward, which are three basic
elements for RL, should be designed depending on the motion to be trained.

The curse of dimensionality(COD) occurs when the dimension of the state
increases and the amount of training data is inappropriately small. The
inappropriately small data for training causes the RL models to be overfitted
and consequently leads to deteriorated performance.

In this paper, we consider training kick-motion with a deep Q-network (DQN).
To train kick-motion successfully, the state should include various
information about the robot and ball on the soccer field, which is described
by their current position, speed, and the direction of the movement of ball or
robot. However, designing the state to include large information can cause the
COD, unless experiences are gained evenly over all areas of the soccer field.
To eliminate the need for the agent to know all the game information over
entire soccer field and to reduce the dimension of the state, we attempt to
design the state for efficient kick-motion training by using a coordinate
transformation matrix (CTM) used to transform coordinates in robotics or
computer vision [11]. The main contributions of this paper are as follows:

  1. 1.

The attempted designing of the state avoids deterioration arising from the COD
and improve the performance of the RL algorithm.

  2. 2.

The RCS taken for AI Soccer eliminates the need for excessive exploration and
reduces the dimension of the state.

##  II Background

###  II-A AI soccer environment

In this paper, experiments are conducted on the AI soccer environment employed
for AI World Cup event [12]. The environment consists of a ball, two teams of
robots and a soccer field. Each team consists of five robots each driven by
two wheels. Like the human soccer game, each team tries to beat the other team
by controlling their robots to take the best action for each game state. The
robots are controlled at each timestep by setting wheel speeds. The two
separately controllable wheels extend the range of actions and increase the
degree of freedom. Value of the best action can be measured by corresponding
reward. Each robot, as an agent, takes an action based on the situation of the
field to get a higher accumulated reward. This is a very challenging task,
because various strategies are needed for the agents to adapt to the vast
state space which is observed in dynamic game situation.

###  II-B Curse of dimensionality (COD)

Traditional RL algorithms suffer from the COD when the dimension of the space
in the environment is excessively high [13]. As the dimension of the state
increases, the number of experiences needed to train an agent grows
exponentially. The lack of experiences makes the RL model overfitted. Although
using function approximation such as DNN alleviates this problem, the
performance of the algorithm is still deteriorated. In the AI soccer
environment, large information describing the (game) state and two-wheeled
actions with a high degree of freedom increase the dimension of the state
space and action space, which readily causes the COD.

###  II-C Deep Q-network (DQN)

The DQN algorithm, which is a widely used RL algorithm, is a deep RL method
that approximates an action-value function, called Q-function by using a DNN.
The Q-function with parameters $\theta$ is defined as follows:

$Q^{\theta}(s,a)=\mathbb{E}_{s,a,\theta}[\sum_{t=1}^{\infty}\gamma^{t-1}r_{t}]$
(1)

where a, s, r, and $\gamma$ are action, state, reward, and discount factor
respectively. DQN is trained to minimize the loss function, which can be
presented as follows:

$L(\theta)=\mathbb{E}[(y_{t}-Q^{\theta}(s_{t},a_{t}))^{2}]$ (2)

The target Q function $y_{t}$ is calculated by reward and next state
$s^{\prime}$ as follows:

$y_{t}=r(s_{t},a_{t})+\gamma
max_{a^{\prime}}Q^{\theta^{-}}(s^{\prime},a^{\prime})$ (3)

where $\theta^{-}$ represents the parameters of a fixed and separate target
network.

Figure 1: Illustration of the absolute coordinate system (ACS). The origin of
the ACS is placed in the center of the field. The orange dotted line presents
the direction of ball movement. The yellow arrow presents the direction the
robot movement.

Figure 2: Illustration of the relative coordinate system (RCS). The origin of
the RCS is placed in the center of the robot, and the direction the robot is
aligned with the $X^{\prime}$-axis.

##  III Relative coordinate system for description of game state

In the task of training the robot to perform kick-motion, the robot chooses at
each timestep the best action out of a set of actions, taking into account the
information about the ball and the robot such as the current coordinates, the
velocity of the ball and the robot, and the approaching angle of the ball to
the robot. To avoid the COD and train the robot to kick the ball at any
location, it is required to use evenly distributed experiences over entire
field. For this purpose, we adopt the relative coordinate system (RCS) to
represent the state, instead of the absolute coordinate system (ACS). The RCS
can alleviate the COD and make training process straightforward. Most
important benefit obtained from the use of the RCS is that the robot can
achieve the task at any location after training at only one location.

In the AI soccer environment, all raw information given to the robot from the
environment is presented in the ACS format. The origin of the ACS is defined
as the center of the field with $X$-axis and $Y$-axis as shown in Fig. 1. In
this case, the state given to the robot is defined as follows:

$S_{ACS}=[R_{x},R_{y},R_{\theta},B_{x},B_{y},V_{R_{x}},V_{R_{y}},V_{B_{x}},V_{B_{y}},\theta]$
(4)

Each element in (4) respectively represents $X$ and $Y$ coordinates of the
robot, the heading angle of the robot, $X$ and $Y$ coordinates of the ball,
$X$ and $Y$ direction velocity of the robot, $X$ and $Y$ direction velocity of
the ball, and the angle of the direction in which the ball moves. The high
dimensional state space formulated in (4) can deteriorate the performance of
RL algorithms.

With the RCS, the dimensional space is reduced and training becomes efficient.
In the RCS, the origin is defined at the center of the robot with
$X^{\prime}$-axis and $Y^{\prime}$-axis and the direction the robot movement
is the same as the direction of $X^{\prime}$-axis shown in Fig. 2. With the
RCS, the values of $R_{x}$, $R_{y}$, and $R_{\theta}$ are always zero since
the origin of the RCS is the center of the robot. The value of $V_{R}{}_{y}$
is also always zero because the robot can take only one of forward, backward,
or stop, during the task of the kick-motion. As a result, the state can be
defined as follows:

$S_{RCS}=[B_{x}^{RCS},B_{y}^{RCS},V_{R_{x}}^{RCS},V_{B_{x}}^{RCS},V_{B_{y}}^{RCS},\theta^{RCS}]$
(5)

where $(\cdot)^{RCS}$ means $(\cdot)$ in the RCS. The relationship between the
RCS and the ACS can be presented using the CTM as shown in (6).

$M=\begin{bmatrix}\cos(R_{\theta})&-\sin(R_{\theta})&R_{x}\\\
\sin(R_{\theta})&\cos(R_{\theta})&R_{y}\\\ 0&0&1\end{bmatrix}$ (6)

The position, velocity, and angle in the RCS can be computed with the CTM and
those in the ACS, and can be shown as (7), (8), and (9).

$\begin{bmatrix}B_{x}^{RCS}&V_{B_{x}}^{RCS}\\\ B_{y}^{RCS}&V_{B_{y}}^{RCS}\\\
1&1\end{bmatrix}=M^{-1}\begin{bmatrix}B_{x}&V_{B_{x}}\\\ B_{y}&V_{B_{y}}\\\
1&1\end{bmatrix}$ (7)

$V_{R_{x}}^{RCS}={V_{R_{x}}}^{2}+{V_{R_{y}}}^{2}$ (8)

$\theta^{RCS}=\theta-R_{\theta}$ (9)

The change of the coordinate system from the ACS to the RCS reduces the
dimension of the state by 60%, which means significantly easier training
process.

Figure 3: Illustration of a situation in which the ball and the robot are in
contact. The blue line connects from the center of the robot to the center of
the ball and the yellow arrow presents the direction toward the target point.

##  IV Experiments

This section describes AI Soccer environment and presents the result of
experiments conducted to show the effect of the COD.

###  IV-A Experimental Environment

Experiments are conducted for the kick-motion task in the AI soccer
environment developed with Webot physics simulator [14]. The goal of the kick-
motion task is that robot kicks the ball to a target point. The robot is
facing the target point and the ball is located at a random point with a
distance of 2m from the robot. With the start of an episode of training, the
ball starts to move at a random speed to a random point near the robot. At
each timestep, the robot chooses an action based on the state. Each episode
consist of 40 experiences over 40 timesteps. The RL framework for this task
can be described as follows:

  * States:

State consists of 6 values shown as (5) in Section III.

  * Rewards:

Reward is given as a sum of 3 types of reward,
$r_{contact}+r_{theta}+r_{velocity}$. The $r_{contact}$ is a reward of
contact, which is -0.1 before the robot contact the ball and 0 after the
contact. The $r_{Theta}$ is a reward according to distance of the contact
point from the center of the robot, which is represented as
$3\cos{\theta_{contact}}$, where $\theta_{contact}$ is the angle between the
blue and the yellow lines in Fig. 3. The $r_{velocity}$ is a reward according
to the speed of the contact with the ball, which is represented as
$3.92V_{contact}$, where $V_{contact}$ is the velocity of the robot at the
contact.

  * Actions:

List of available actions: going straight forward with steady-state speeds of
2.55 (m/s), stop, going straight backward with steady-state speeds of -2.55
(m/s).

###  IV-B Experimental Results

The performances with the kick-motion task are evaluated in terms of the sum
of the reward for each training episode. To reduce granularity of the reward
variation, the moving average of the sum of the reward across 250 episodes is
taken.

Figure 4: Total reward according to the number of dummy variables. It shows
that increasing the dimension of state space slows down the training speed and
as a result degrades performance. The blue line represents the model with the
RCS and the others represent comparative models. Especially, the purple line
is equivalent to the model with the ACS.

To present the intuitive understanding of the COD, models used to compare with
the proposed method are trained with states including ’dummy variables’, which
are some meaningless random variables. The states for the comparison models
are the state in the RCS with ’dummy variables’. In Fig. 4, the name of ’$+n$’
means that $n$ dummy variables are added to the state.

Figure 4 shows the results of the experiments. The (total accumulated) reward
of the action for RCS based state exceeds -4 from 1500-th episode, which is
about 500 episodes faster than by ACS based method. The RCS based method
achieves the highest reward value of 6.8, although the reward value of the RCS
based method is similar to the others in the latter part of the training,
e.g., 3600-4700 episodes.

##  V Conclusion

Efficient training of the RL algorithms with high dimensional state space is
hard to achieve, due to the COD. In this paper, a method to reduce the
dimension of the state space for training the kick-motion task in the AI
soccer environment is attempted. The attempted method uses the CTM to convert
the state in the ACS to the state in the RCS. The dimension of the state space
required for the agent to achieve the kick-motion is reduced from ten to six
since 4 values in the state are always zero by the RCS based method. The
experimental result shows that the kick-motion can be achieved with the
reduced state by the RCS based method. The result obtained along with the
comparative models with different number of variables, which are trained with
the states including dummy variables, shows that the lower dimension of the
state space can avoid deterioration due to the COD.

## Acknowledgment

This work was supported by the Institute for Information communications
Technology Promotion (IITP) grant funded by the Korean government (MSIT)
(No.2020-0-00440, Development of Artificial Intelligence Technology that
continuously improves itself as the situation changes in the real world).

## References

  * [1] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * [2] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin Riedmiller.  Playing atari with deep reinforcement learning.  arXiv preprint arXiv:1312.5602, 2013.
  * [3] David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al.  Mastering the game of go with deep neural networks and tree search.  nature, 529(7587):484–489, 2016.
  * [4] David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, et al.  Mastering the game of go without human knowledge.  nature, 550(7676):354–359, 2017.
  * [5] Minah Seo, Luiz Felipe Vecchietti, Sangkeum Lee, and Dongsoo Har.  Rewards prediction-based credit assignment for reinforcement learning with sparse binary rewards.  IEEE Access, 7:118776–118791, 2019.
  * [6] Luiz Felipe Vecchietti, Taeyoung Kim, Kyujin Choi, Junhee Hong, and Dongsoo Har.  Batch prioritization in multigoal reinforcement learning.  IEEE Access, 8:137449–137461, 2020.
  * [7] Luiz Felipe Vecchietti, Minah Seo, and Dongsoo Har.  Sampling rate decay in hindsight experience replay for robot control.  IEEE Transactions on Cybernetics, 2020.
  * [8] Inhwan Kim, Sarvar Hussain Nengroo, and Dongsoo Har.  Reinforcement learning for navigation of mobile robot with lidar.  In 2021 5th International Conference on Electronics, Communication and Aerospace Technology (ICECA), pages 148–154. IEEE, 2021.
  * [9] Woohyeon Moon, Bumgeun Park, Sarvar Hussain Nengroo, Taeyoung Kim, and Dongsoo Har.  Path planning of cleaning robot with reinforcement learning.  arXiv preprint arXiv:2208.08211, 2022.
  * [10] Taeyoung Kim, Luiz Felipe Vecchietti, Kyujin Choi, Sanem Sariel, and Dongsoo Har.  Two-stage training algorithm for ai robot soccer.  PeerJ Computer Science, 7:e718, 2021.
  * [11] Richard Hartley and Andrew Zisserman.  Multiple view geometry in computer vision.  Cambridge university press, 2003.
  * [12] Chansol Hong, Inbae Jeong, Luiz Felipe Vecchietti, Dongsoo Har, and Jong-Hwan Kim.  Ai world cup: Robot-soccer-based competitions.  IEEE Transactions on Games, 13(4):330–341, 2021.
  * [13] Abhijit Gosavi.  Reinforcement learning: A tutorial survey and recent advances.  INFORMS Journal on Computing, 21(2):178–192, 2009.
  * [14] Olivier Michel.  Cyberbotics ltd. webots™: professional mobile robot simulation.  International Journal of Advanced Robotic Systems, 1(1):5, 2004\.