Year: 2021

Abstract:
Finding a suitable density function is essential for density-based clustering algorithms such as DBSCAN and DPC. A naive density corresponding to the indicator function of a unit d-dimensional Euclidean ball is commonly used in these algorithms. Such a density suffers from an inability to capture local features in complex datasets. To tackle this issue, we propose a new kernel diffusion density function, which is adaptive to data of varying local distributional characteristics and smoothness. Furthermore, we develop a surrogate that can be efficiently computed in linear time and space and prove that it is asymptotically equivalent to the kernel diffusion density function. Extensive empirical experiments on benchmark and large-scale face image datasets show that the proposed approach not only achieves a significant improvement over classic density-based clustering algorithms but also outperforms the state-of-the-art face clustering methods by a large margin.

text: Finding a suitable density function is essential for density-based clustering algorithms such as DBSCAN and DPC. A naive density corresponding to the indicator function of a unit d-dimensional Euclidean ball is commonly used in these algorithms. Such a density suffers from an inability to capture local features in complex datasets. To tackle this issue, we propose a new kernel diffusion density function, which is adaptive to data of varying local distributional characteristics and smoothness. Furthermore, we develop a surrogate that can be efficiently computed in linear time and space and prove that it is asymptotically equivalent to the kernel diffusion density function. Extensive empirical experiments on benchmark and large-scale face image datasets show that the proposed approach not only achieves a significant improvement over classic density-based clustering algorithms but also outperforms the state-of-the-art face clustering methods by a large margin.

heading: 1 INTRODUCTION
text: Density-based clustering algorithms are now widely used in a variety of applications, ranging from high energy physics (Tramacere & Vecchio, 2012; Rovere et al., 2020), material sciences (Marquis et al., 2019; Reza et al., 2007), social network analysis (Shi et al., 2014; Khatoon & Banu, 2019) to molecular biology (Cao et al., 2017; Ziegler et al., 2020). In these algorithms, data points are partitioned into clusters that are considered to be sufficiently or locally high-density areas with respect to an underlying probability density or a similar reference function. We call them density functions throughout this paper. These techniques are attractive to practitioners, due to their nonparametric feature, which leads to flexibility in discovering clusters that have arbitrary shapes, whilst classic methods such as k-means and k-medoids (Friedman et al., 2001) can only detect convex (e.g., spherical) clusters. Seminal work in the context of density-based clustering includes DBSCAN (Ester et al., 1996) and DPC (Rodriguez & Laio, 2014), among many others (Ankerst et al., 1999; Cuevas et al., 2001; Comaniciu & Meer, 2002; Hinneburg & Gabriel, 2007; Stuetzle, 2003).
Most density-based clustering algorithms implicitly identify cluster centers and assign remaining points to the clusters by connecting with the higher density point nearby. To proceed with these methods it requires a density function, which is usually an estimate of the underlying true probability density or some variants of it. For example, a popular choice is the naive density function that is carried out by simply calculating the number of data points covered in the ε-neighborhood of each x. Note that such densities are not adaptive to different distribution regions. One of the challenging scenarios is when clusters in the data have varying local features, for example, size, height, spread, and smoothness. Therefore, the resulting density function has a tendency to flatten the peaks and valleys in the data distribution, which leads to underestimation of the number of clusters (see Figure 1). Many heuristics variations of DBSCAN and DPC have been proposed to magnify the local features, thus making the clustering task easier (Campello et al., 2013; Chen et al., 2018; Ertöz et al., 2003; Zhu et al., 2016). Most of these methods can be viewed as performing clustering on certain transformations of the naive density function. However, if the naive density function itself is quite problematic in the first place, these methods will become less effective.
Moreover, even if we apply adaptive alternatives to modify the classic density functions, there are other contentious issues of the generally used linear kernel density estimator (KDE). It often suffers from severe boundary bias (Marron & Ruppert, 1994) and is acknowledged as computationally
expensive. These phenomenon prevent the classic density functions being practically useful and reliable, especially for large-scale and complex clustering tasks.
To overcome these problems in density-based clustering algorithms, in this paper we propose a general approach to build the so-called kernel diffusion density function to replace classic density functions. The key idea is to construct the density from a user-specified bivariate kernel that has desired local adaptive properties. Instead of using the naive density function and its variants, we utilize the bivariate kernel to derive a transition probability. A diffusion process is induced by this transition probability, which admits a limiting and stationary distribution. This limiting distribution serves as a plausible density function for clustering with reduced error.
Under this framework, we provide examples of symmetric and asymmetric bivariate kernels to construct the kernel diffusion density function, which can tackle clustering complex and locally varying data. We apply the resulting adapted DBSCAN and DPC algorithms to widely different empirical datasets and show significant improvement in each of these analyses. The main contributions of this paper are summarized below:
• We introduce new bivariate kernel functions and construct the associated kernel diffusion processes. Based on the diffusion process, we propose a kernel diffusion density function to adapt density-based clustering algorithms such as DBSACN and DPC, which attains accuracy in the presence of varying local features.
• We derive a computationally much more efficient surrogate, and show analytically it is asymptotic equivalent to the proposed kernel diffusion density function.
• By extensive experiments, we demonstrate the superiority of kernel diffusion density function over naive density function and its variants when applying to DBSCAN and DPC, and show it outperforms state-of-the-art GCN-based methods on face clustering tasks.

heading: 2 RELATED WORK
text: Density-Based Clustering There is vast literature on adapting density-based clustering algorithms to tackle large variations in different clusters in the data. DPC itself is such a refinement of DBSCAN, as it determines cluster centers not only by highest density values but also by taking into account their distances from each other, thus has a generally better performance in complex clustering tasks. Other attempts include rescaling the data to have relative reference measures instead of KDE (Zhu
et al., 2016; Chen et al., 2018), and using the number of shared-nearest-neighbors between two points to replace the geometric distance (Ertöz et al., 2003).
Diffusion Maps The technique of diffusion maps (Coifman et al., 2005; Coifman & Lafon, 2006) gives a multi-scale organization of the data according to their underlying geometric structure. It uses a local similarity measure to create a diffusion process on the data which integrates local geometry at different scales along time t. Generally speaking, the diffusion will segment the data into several smaller clusters in small t and group data into one cluster for large t. Applying eigenfunctions at a carefully selected time t leads to good macroscopic representations of the data, which is useful in dimension reduction and spectral clustering (Nadler et al., 2005).
Face Clustering Face clustering has been extensively studied as an important application in machine learning. Traditional algorithms include k-means, hierarchical clustering (Sibson, 1973) and ARO (Otto et al., 2017). Many recent works take advantage of supervised information and GCN models, achieving impressive improvement comparing to traditional algorithms. To name a few, CDP (Zhan et al., 2018) proposes to aggregate the features extracted by different models; L-GCN (Wang et al., 2019) predicts the linkage in an instance pivot subgraph; LTC (Yang et al., 2019) generates a series of subgraphs as proposals and detects face clusters thereon; and GCN(V+E) (Yang et al., 2020) learns both the confidence and connectivity by GCN. In this paper we demonstrate that the proposed density-based clustering algorithm with kernel diffusion, as a general clustering approach, even outperforms theses state-of-the-art methods that are especially designed for face clustering.

heading: 3 PRELIMINARIES

heading: 3.1 NOTATIONS
text: Let the dataset D = {x1, . . . , xn} ⊂ Rd be n i.i.d samples drawn from a distribution measure F with density f on Rd. Let Fn denote the corresponding empirical distribution measured with respect to D, i.e., Fn(A) = 1n ∑n i=1 1A(xi), where 1A(·) denotes the indicator function of set A. We write ||u|| as the Euclidean norm of vector u. Let B(x, ε) and Vd denote the d-dimensional ε-ball centered at x and the volume of the unit ball B(0, 1), respectively. Let Nk(x) denote the set of k-nearest neighbors of point x within the dataset D.

heading: 3.2 DENSITY FUNCTION
text: Density-based algorithms perform clustering by specifying and segmenting high-value areas in a density function denoted by ρ. Usually, we calculate each of ρ(xi), and then identify cluster centers with (locally) highest values. Many popular algorithms such as DBSCAN and DPC employ the following naive density function:
ρnaive(x) = 1
nεd ∑ y∈D 1B(x,ε)(y) Vd . (1)
The naive density function ρnaive is actually an empirical estimation of f for carefully chosen ε. It is easy to observe, for clustering purpose we only care about ρnaive(x) up to a normalising constant, which makes it simply equivalent to counting the total number of data points in the ε-ball around x.
In practice, the data distribution may be very complex and contains varying local features that are difficult to be detected. The naive density in (1) with the same radius ε for all x usually suffers from unsatisfactory empirical performance, for example, failing to identify small clusters with fewer data points. One possible way to alleviate this problem is through a transformation into the following local contrast (LC) function (Chen et al., 2018):
ρLC(x) = 1
n ∑ y∈Nk(x) 1ρnaive(x)>ρnaive(y). (2)
In this way, ρLC compares the density of each data point with its k-nearest neighbors. To see the benefit of LC, let us consider x to be a cluster center. After local contrasting, ρLC(x) is likely to reach the value of k regardless of the size of this cluster.
However density functions like ρLC still highly depend on the underpinning performance of ρnaive. This restricts their applications in clustering data with challenging local features.

heading: 4 METHODOLOGY
text: In this section, we present a new type of density-based clustering algorithm, based on the notion of kernel diffusion density function. Towards this end, we will introduce a kernel diffusion density function, which takes account of local adaptability and is well-suited for clustering purpose. We provide details on how to derive this density function from a diffusion process induced by bivariate kernels. We also provide a surrogate density function that is computationally more efficient.

heading: 4.1 DIFFUSION PROCESS AND KERNEL DIFFUSION DENSITY
text: Considering a bivariate kernel function k : D ×D → R+, such that:
• k(x, y) is positive semi-definite, i.e., k(x, y) ≥ 0. • k(x, y) is Fn-integrable with respect to both x and y.
We define d(x) = n ∫ D k(x, y)dFn(y) as a local measure of the volume at x and define
p(x, y) = k(x, y)
d(x) . (3)
It is easy to see that p(x, y) satisfies the conservation property n ∫ D p(x, y)dFn(y) = 1. As a result, p(x, y) can be viewed as a probability for a random walk on the dataset from point x to point y, which induces a Markov chain on D with n× n transition matrix P = [p(x, y)]. This technique is standard in various applications, known as the normalized graph Laplacian construction. For example, we can view D as a graph, L = I −P as the normalized graph Laplacian, and d(x) as a normalization factor. For t ≥ 0, the probability of transiting from x to y in t time steps is given by P t, the t-th power of P . Running the Markov chain forward in time, we observe the dataset at different scales, which is the diffusion process Xt on D. Let ρ(x, t): D × R+ → R+ be the associated probability density, which is governed by the following second-order differential equation with initial conditions:{
∂ ∂tρ(x, t) = −Lρ(x, t), ρ(x, 0) = ϕ0(x),
(4)
where ϕ0(x) is a probability density at time t = 0. In practice we can use any valid choice of ϕ0(x), e.g. the uniform density.
To give an explicit example of the diffusion process, consider a sub-class of k, i.e., isotropic kernels, where k(x, y) = K(||x− y||2/h) for some function K : R → R+. Here we can dual interpret h as a scale parameter to infer local information and as a time step h = ∆t at which the random walk jumps. Then we can define the forward Chapman-Kolmogorov operator TF as
TF (x) = n ∫ D p(x, y)ϕ0(y)dFn(y).
Note that TF is the data distribution at time t = h, thus can be viewed as continuous analogues of the left multiplication by the transition matrix P . Letting h → 0, the random walk converges to a continuous diffusion process with probability density evolves continuously in t. In this case, we can explicitly write the second-order differential equation in (4) as:
∂ ∂t ρ(x, t) = lim h→0 ρ̂(x, t+ h)− ρ(x, t) h = lim h→0 TF − I h ρ(x, t), (5)
where Lh = limh→0 (TF − I)/h is the conventional infinitesimal generator of the process. Now we are ready to introduce our kernel diffusion density function.
Definition 1. (Kernel diffusion density function) Suppose the Markov chain induced by P is ergodic, we define the kernel diffusion density function as the limiting probability density of the diffusion process Xt, i.e.,
ρKD(x) = lim t→∞ ρ(x, t). (6)
Intuitively, on the one hand, with increased values of t we expect the diffusion process Xt gradually reveals the geometric structure (such as high-density regions) of the data distribution F . To see this, note that the transition probability P reflects connectivity between data. We can interpret a cluster as an underlying geometric structure in which the probability of staying in this region is high during a transition. In the diffusion process, the probability of following a path along a structure usually increases with t, as the involved data points are dense and highly connected. Therefore, the path consists of short and high probability jumps. Whilst paths that do not follow any structure consists of long and low probability jumps, which lowers their overall probability. As a result, geometry structures of F is magnified during the diffusion.
On the other hand, by talking certain sophisticated forms of k(x, y) that take into account of local adaptivity, we also slow down the diffusion to avoid trivial geometry structures such as one big cluster for all the data points. In this way, we can eventually identify the correct geometry structure at the right scale.

heading: 4.2 LOCALLY ADAPTIVE KERNELS
text: To address the local adaptability in kernel diffusion density function, we propose the following two bivariate kernels. Both of them are very simple variations of the most commonly used classic kernels.
Symmetric-Gaussian kernel:
k(x, y) = exp ( − ∥x− y∥ 2
h
) 1B(x,ε)(y). (7)
Here h and ϵ are both hyper-parameters. We call this kernel symmetric since k(x, y) = k(y, x) .
Asymmetric-Gaussian kernel:
k(x, y) = exp ( − ∥x− y∥ 2
h
) 1Nk(x)(y). (8)
Here h and k are hyper-parameters. Note that in this case k(x, y) is asymmetric as y ∈ Nk(x) does not imply x ∈ Nk(y). Bivariate kernels defined in (7) and (8) are just combinations of classic Gaussian kernel and εneighbourhood or k-nearest neighbours kernels, respectively. With these simple combinations, we truncate Gaussian kernel at local areas, and the contribution of each point y to the construction of the density function ρKD(x) depends not only on the distance ||y − x|| but also on the local geometry structure around x. Hence, the new kernels are adaptive at different x, which is expected to lead to better clustering performance against local features. We remark that the Asymmetric-Gaussian kernel takes into account a varying neighborhood around each x, thus is more adaptive comparing to the Symmetric-Gaussian kernel.
Although here we only provide two examples of locally adaptive kernels, other options can be easily created in a similar spirit under this framework, e.g., changing the Gaussian kernels to other kernels or changing the ε-neighbourhood (k-nearest neighbours) kernels to other locally truncated functions. Once k(x, y) is determined, we can derive the corresponding density function ρKD. Next, we just need to simply apply any density clustering procedure like DPC or DBSCAN based on ρKD instead of the naive density function ρnaive.
In Section 5, we assess the empirical performance of the proposed kernel diffusion density function with the above two locally adaptive kernels. They outperform existing density-based algorithms and other state-of-the-art methods.

heading: 4.3 FAST KERNEL DIFFUSION DENSITY
text: The kernel diffusion density function ρKD can be calculated as the stationary distribution of a Markov chain induced by the transition matrix P . Numerically, we can solve it by iteratively right multiplying P with ρ(x, t) until convergence, or applying a QR decomposition on P . These methods are expensive in terms of computational cost, especially when the sample size n is large.
To tackle this problem, we propose the following surrogate of ρKD(x) which is computationally more efficient. Definition 2. (Fast kernel diffusion density function) Let p(y, x) be the transition probability from point y to point x, as defined in equation (3). We define the fast kernel diffusion density function as
ρFKD(x) = ∫ D p(y, x)dFn(y), (9)
It is straightforward that ρFKD can be obtained in linear time and memory space, as we only need to compute the column averages of matrix P .
Here we show that ρFKD is not only computationally efficient but also suitable for detecting local features. This is illustrated through the following Theorem 1. Consider a special case that k(x, y) = 1B(x,ε)(y). Then it is easy to verify that
ρFKD(x) = 1
Cd ∑ y∈B(x,ε)
1
ρnaive(y) ,
where Cd = nεdVd is a normalising constant. In this way, we build a connection between ρFKD and the naive density function ρnaive in this special example. Theorem 1. Consider the above special case that k(x, y) = 1B(x,ε)(y). In addition, assume the dataset D = {x1, ..., xn} can be split into m disjoint clusters: i.e., X = D1 ⋃ ... ⋃ Dm and for each x ∈ D, B(x, ε) only contain data points that belong to the same cluster as x. Denote ρ̄j =
1 |Dj | ∑ x∈Dj ρFKD(x) as the average density in cluster j. We have
ρ̄1 = · · · = ρ̄m = 1.
Theorem 1 demonstrates that the averaged ρFKD in each cluster are the same regardless of cluster sizes and other local features. This shows that ρFKD elevates the density of small clusters, which is essential for finding the density peaks of small clusters.
Previously we claim that ρFKD is a surrogate of the kernel diffusion density ρKD. Next, we want to reveal the relationship between these two density functions from an asymptotic viewpoint. To proceed, we will need the following assumption. Assumption 1. There exists some positive constant c < 1 that is independent of n, such that ρFKD(x) ≤ c holds for every x ∈ D.
This is a very mild assumption, since it always holds that ρFKD(x) < 1, and the average of ρFKD(x) over the dataset is ∫ D ρFKD(x)dFn(x) = 1/n, which vanishes as n → ∞. Now we are ready to present the following theorem that characterise the closeness between ρFKD and ρKD. Theorem 2. Suppose that Assumption 1 holds and the Markov chain induced by the kernel k(x, y) is ergodic. We have
ρKD(x)
ρFKD(x)
a.s.−→ 1
As shown in the Appendix, the almost sure convergence in Theorem 2 is of a fast rate at n−1. Thus it is safe for us to use it to replace ρKD in finite sample experiments. This result is also verified by our numerical experiments in Section 5.

heading: 5 EXPERIMENTS
text: In this section, we empirically evaluate the proposed kernel diffusion density functions against ρnaive and ρLC in density-based clustering algorithms, and also compare them with other state-of-the-art
methods. We denote ρsymKD and ρ sym FKD as the kernel diffusion density functions and its fast surrogate, with symmetric-Gaussian kernel, respectively. Similarly, we denote ρasymKD and ρ asym FKD as the proposed two density functions with asymmetric-Gaussian kernel, respectively. We examine their performance on a wide range of datasets. The clustering results are measured Pairwise F-score (Banerjee et al., 2005), BCubed F-score (Amigó et al., 2009) or NMI (Cover, 1999).

heading: 5.1 PERFORMANCE ON BENCHMARK DATASETS
text: We now discuss the performance on 13 benchmark datasets (∼100 to ∼5, 000 data points) from UCI repository. The metadata is summarised in the Appendix.
As summarised in Table 1, both ρsymKD and ρ asym KD uniformly outperform ρnaive and ρLC in terms of clustering accuracy in terms of F-scores. results based on NMI are deferred to the Appendix. The proposed kernel diffusion density functions with asymmetric Gaussian kernel, ρasymKD , which enjoys better local adaptivity analytically, achieves the best results on most datasets and outperforms ρnaive and ρLC by a large margin. It is worth noticing that the two fast surrogates, ρ sym FKD and ρ asym FKD , achieve comparable results with their original counterparts, ρsymKD and ρ asym KD . In the Appendix similar results are observed for the same set of density functions applied to DBSCAN.

heading: 5.2 PERFORMANCE ON FACE IMAGE DATASETS
text: Clustering face images according to their latent identity becomes an important application in recent years. It is challenging in the sense that face image datasets usually contain thousands of identities, corresponding to thousands of clusters. Meanwhile, the number of images for each identity (cluster) is quite different, corresponding to the variety of cluster sizes. We assess the performance of the proposed approach on two popular face image datasets: emore_200k (Zhan et al., 2018) and MS1M (Guo et al., 2016).
emore_200k. The dataset contains 2,577 identities with 200,000 images following the protocol in Zhan et al. (2018). We set ε to 0.8, k to 200, and h to 0.5 for density-based methods. Results are summarized in Table 2. The proposed diffusion density functions are applied to DPC, and compared with k-means, HAC (Sibson, 1973), ARO (Otto et al., 2017), and CDP (Zhan et al., 2018). Again, we observe significant improvement in the proposed density functions over ρnaive and ρLC. It is also worth pointing out that, density-based clustering with proposed kernel diffusion density functions also outperform the
state-of-the-arts approaches such as CDP by a large margin.
MS1M. The dataset contains 8,573 identities with around 584,000 images following the protocols in Yang et al. (2020). We set ε to 0.8, k to 200, and h to 0.5 for density-based methods. We reported the results of clustering performance in Table 3. Precision versus Recall curves for different density functions (applied to DPC) are plotted in Figure 2. In Table 3, the proposed kernel diffusion density functions outperform ρnaive and ρLC. Note that GCN-based methods such as L-GCN (Wang et al., 2019), LTC (Yang et al., 2019) and GCN (V+E) (Yang et al., 2020) achieve generally better clustering performance than unsupervised methods due to their supervised nature. However, it is quite encouraging to see that the proposed kernel diffusion approaches, although are also unsupervised clustering methods, considerably outperform the GCN-based methods.

heading: 5.3 SENSITIVITY ANALYSIS
text: Next, we examine the sensitivity of the proposed kernel diffusion density functions to hyperparameters and compare it with ρnaive and ρLC. The results are obtained via extensive experiments on emore_200k and MS1M, which are shown in Figure 3. We can see that the clustering performance of ρsymKD is much more stable than ρnaive and ρLC when we vary the value of ε. Whilst ρ asym KD is robust to the parameter k, and both ρsymKD and ρ asym KD are quite robust to the parameter h.
5.4 COMPUTATIONAL COST
We carried out a series of experiments on MS1M to demonstrate the computational efficiency of the fast surrogate ρFKD in terms of time and space. With a collection of subsampled data from MS1M at different percentile levels, we run both the kernel diffusion density ρKD and the fast surrogate ρFKD. As we can observe from Figure 4, the running time and memory usage of ρKD increase dramatically with the sample size. Whilst ρFKD retains a very low level of computational cost. This suggests that ρFKD, which achieves an excellent computational efficiency, should be favored in practice.

heading: 6 CONCLUSION
text: Density-based clustering has a profound impact on machine learning and data mining. However, the underpinning naive density function suffers from detecting varying local features, causing extra errors in the clustering. We propose a new set of density functions based on the kernel diffusion process to resolve this problem, which is adaptive to density regions of varying local distributional features. We demonstrate that DBSCAN and DPC adapted by the proposed approach have improved clustering performance comparing to their classic versions and other state-of-the-art methods.

heading: A APPENDIX
text: In this supplementary file, we provide technical proofs of the theoretical results in Section 4.3, and present extra empirical experiments regarding our kernel diffusion approach with symmetric and asymmetric Gaussian kernels applied to DBSCAN. All the numerical experiments are carried out on a standard work station with a Intel 64-cores CPU and two Nvidia P100 GPUs.
A.1 PSEUDO CODE FOR DBSCAN AND DPC.
Algorithm 1 DBSCAN 1: Input: SetOfPoints X , Eps ε, MinPts k 2: H:={x ∈ X : |B(x, ε) ∩X| ≥ k}; 3: G:=undirected graph with vertices H and edge between x, x′ ∈ H if |x− x′| ≤ ε; 4: Output: connected components of G
The connected compoenents of the graph are determined a clusters,a dn the remaining points are unclustered and considered as noise-points.
Algorithm 2 DPC 1: Input: SetOfPoints, TruncDis dc 2: Compute di,j for ∀i, j ∈ SetOfPoints; 3: For i from 1 to SetOfPoints.size: 4: ρi := ∑ i ̸=j 1{di,j−dc}
5: δi := minj:ρj>ρi(di,j) 6: Plot decision map M with ρ as the horizontal axis and δ as the vertical axis; 7: Mark Point i with relatively higher ρi and δi as a cluster center; 8: Mark Point i with relatively lower ρi but relatively higher δi as a noise point; 9: Assign the rest point with the label the same as the nearest cluster center;
10: Output: SetOfPoints
A.2 PROOFS OF THEORETICAL RESULT.
Proof of Theorem 1. Since {D1, · · · , Dm} are disjoint, we have p(x, y) = 0 if x and y belong to different clusters. By the definition of matrix P , for each x ∈ Dj , we have∫
D
p(x, y)dFn(y) = 1,
which implies that ∫ x∈Dj ∫ y∈Dj p(x, y)dFn(x)dFn(y) = |Dj |.
Therefore,
ρ̄j |Dj | = ∫ x∈Dj ∫ y∈Dj p(x, y)dFn(x)dFn(y) = |Dj |,
which implies that ρ̄j = 1 for any j = 1, , . . . ,m.
Before proceeding to the proof of Theorem 2, we need following auxiliary lemma that relates the stationary distribution of a Markov chain to an arbitrary vector g. Lemma A.1. Let P be transition probability matrix of a finite inreducible discrete time Markov chain with n states, which admits a stationary distribution, denoted by vector π. We write e = (1, . . . , 1)T ∈ Rn as the a column vector of ones. The following holds for any vector g such that gT e ̸= 0:
(1) (I − P + egT ) is non-singular.
(2) Let H = (I − P + egT )−1, then πT = gTH.
Proof. Since π is the stationary distribution, we have πT e = 1. Applying Theorem 3.3 in (Hunter, 1982) yields that matrix (I − P + egT ) is non-singular. Next recall that πTP = πT , therefore we have
πT (I − P + egT ) = πT − πTP + πT egT
= πT egT
= gT ,
which implies πT = gTH .
Proof of Theorem 2. Note that for for each x ∈ D, the linear reference function ρFKD(x) =∫ D p(y, x)dFn(y) is the corresponding column average of the transition matrix P . We write the i-th column vector of P as pi = ( p(x1, xi), . . . , p(xn, xi) )T .
Therefore ρFKD(xi) = eT pi/n
Since the Markov chain induced by the kernel k(x, y) is ergodic, the density ρ(x, t) of the diffusion process Xt converges to the limiting stationary distribution of the Markov chain, denoted by π.
We can write the n-vectors of g and π in the following form:
g = (g1, . . . , gn) T = n ( ρFKD(x1), . . . , ρFKD(xn) )T and π = ( ρKD(x1), . . . , ρKD(xn) )T ,
where gi = eT pi is the i-th column sums of matrix P . As a result, we have∫ D ρFKD(x)dFn(x) = 1 n eT g = 1, and n ∫ D ρKDdFn(x) = e Tπ = 1.
By the definition of g, we know
(egT )2 = negT and eTP = gT .
It follows from Lemma A.1 that (I − P + egT ) is non-singular and πT = gTH, where H = (I − P + egT )−1. We define B = I + egT . By simple algebra calculation, we can find B is non-singular with
B−1 = I − eg T
n+ 1 .
As a result, it is easy to see that that gTB−1 = gT
n+ 1 and
H−1 = B − P = (I − PB−1)B.
Use the Neumann series, we have
H = B−1(I − PB−1)−1 = B−1 ∞∑ i=0 (PB−1)i.
Thus
πT − gT /n = gT ( H − I
n
) = gT [ B−1 ∞∑ i=0 (PB−1)i − I n ] .
Since we assume for any x ∈ D, ĝ(x) < c for some 0 < c < 1. This leads to
gT pj ≤ nceT pj = ncgj .
Therefore, let κj be the j-th compoenent of gTPB−1, it is straightforward
κj ≤ nc
n+ 1 gj ≤ cgj .
This implies for every x ∈ D, |ρKD(x)− ρFKD(x)| ≤ ρFKD(x) ∣∣∣∣ 1n+ 1 ∞∑ i=0 ci − 1 n ∣∣∣∣ ≤ ρFKD(x) ∣∣∣∣ 1(n+ 1)(1− c) − 1n ∣∣∣∣.
Hence we have
P (
lim n→∞
ρKD(x) ρFKD(x) = 1
) = 1,
which completes the proof.
A.3 ADDITIONAL NUMERICAL EXPERIMENTS
Naive density with different bandwidths. To illustrate the fail of naive density function in scenario as in Figure 1, we also plot it with a range of different values of hyperparameters ϵ below. We can observe that it is difficult to detect the three true underlying clusters in all the cases.
Hyperparameters The parameter ε (radius of the ball, used in ρnaive, ρLC, ρsymKD and ρ sym FKD) is tuned by searching within the range between 0.1 and 1 with am increment of 0.1, parameter k (number of nearest neighbors, used in ρLC, ρ asym KD and ρ asym FKD ) is tuned by searching within the range between 10% and 50% number of samples, with an increment of 10%.
Metadata of benchmark datasets. The number of samples n, the number of clusters c, and feature dimension d for each benchmark dataset are listed in Table 4 below.
Benckmark datasets with DBSCAN. We provide the performance of the conventional density functions, ρnaive and ρLC, and the proposed kernel diffusion density functions with symmetric and asymmetric Gaussian kernels, ρ∗KD and ρ ∗ FKD (∗ ∈ {sym, asym}), applied to DBSCAN on 13 benchmark datasets. The results are summarised in Table 6. Similar to DPC, we see that both ρsymKD and ρasymKD uniformly outperform ρnaive and ρLC in terms of clustering quality. ρ asym KD , which has better local adaptivity analytically, achieves the best results on most datasets and outperforms others by a significant margin in Breast-o, Control, Haberma and Seeds.
NMI for benchmark datasets. Below we present in Table 6 the clustering results for benchmark datasets based on NMI metric.
Number of clusters. In Table 7, we present the number of clusters returned by the density-based methods for the benchmark datasets. It can be observed that clustering with the proposed diffusion density functions returned a significantly better estimate of the number of clusters, comparing to that with classic density functions such as ρnaive and ρLC .

Review summaries:
The authors propose a new density function for density clustering models such as DBSCAN and Density Peaks Clustering (DPC). Given a kernel $K$ and the corresponding normalized random walk transition matrix $P = diag(K\mathbf{1})^{-1}K$, the authors propose to use the density function which corresponds to the stationary distribution attained in the limit. Since this function is expensive to compute, the authors propose a surrogate density function which is easier to compute.  The main experiments compare various density functions for DPC on UCI classification datasets and the DPC on density functions to face detection methods on two face identification datasets. 

This paper proposes a different density function for popular density-based clustering algorithms, i.e. truncated symmetric and asymmetric Gaussian kernels for DBSCAN and DPC. It offers some diffusion-based argument for why that is a better density function for clustering than naive function, and provides a computationally more efficient surrogate function as well.  Performance of proposed method on many clustering experiments are provided.

The paper introduces a new approach to performing clustering, based on constructing a Markov process on the input space, in which transitions are determined according to normalised kernel similarity scores, as in Spectral Clustering (SC). Unlike SC, however, the proposal utilises the stationary distribution of the process, or an approximation thereof, evaluated on the data set, as a density function to be used within existing density based clustering algorithms, such as DBSCAN and Density Peaks Clustering (DPC). Experimental results are given to illustrate the methods potential improvement of standard kernel density estimates, and also its potential improvement over purpose-built models for face clustering.

The paper presents a methodology for defining the density reference function required in density-based clustering methods such as DBSCAN and DPC. The density function is obtained as the limiting probability density of a diffusion process and it takes into account the local characteristics of the dataset. A surrogate density is also proposed that is fast to compute. The method is compared to typical density functions used in DBSCAN and DPC.