Construction of the similarity matrix for the spectral clustering method: numerical experiments

TL;DR

This paper investigates how the construction of the similarity matrix, especially the choice of scale parameter and sparsity, affects the performance of spectral clustering through extensive numerical experiments.

Contribution

It introduces methods for selecting the scale parameter and constructing the similarity matrix based on graph properties like the MST, improving spectral clustering accuracy.

Findings

Optimal scale parameter selection enhances clustering performance.

Graph-based similarity matrix construction outperforms traditional methods.

Numerical experiments validate the proposed approaches on various datasets.

Abstract

Spectral clustering is a powerful method for finding structure in a dataset through the eigenvectors of a similarity matrix. It often outperforms traditional clustering algorithms such as $k$-means when the structure of the individual clusters is highly non-convex. Its accuracy depends on how the similarity between pairs of data points is defined. Two important items contribute to the construction of the similarity matrix: the sparsity of the underlying weighted graph, which depends mainly on the distances among data points, and the similarity function. When a Gaussian similarity function is used, the choice of the scale parameter $\sigma$ can be critical. In this paper we examine both items, the sparsity and the selection of suitable $\sigma$'s, based either directly on the graph associated to the dataset or on the minimal spanning tree (MST) of the graph. An extensive numerical experimentation on artificial and real-world datasets has been carried out to compare the performances of the methods.

Tables

Table 1. Table 1: Averaged sparsity levels of the methods.

	method	Set A	Set B	Set R	Set U
	$F_1$	52%	74%	67%	44%
	$F_2$	19%	24%	40%	24%
	$F_3$	16%	18%	39%	18%
	$E_i$	92%	95%	93%	94%
$K=1+\lfloor {\log }_{2}n\rfloor$	$N_i$	94%	95%	98%	95%
	$M_i$	96%	97%	99%	97%
	$E_i$	88%	92%	90%	91%
$K=1+\lfloor \sqrt{n}\rfloor$	$N_i$	90%	91%	94%	91%
	$M_i$	93%	94%	97%	95%

Table 2. Table 2: Averaged accuracies α ~ ~ 𝛼 \widetilde{\alpha} and ranking of the methods.

		Set A		Set B		Set R		Set U
	method	$\tilde{\alpha }$	$\rho$	$\tilde{\alpha }$	$\rho$	$\tilde{\alpha }$	$\rho$	$\tilde{\alpha }$	$\rho$
	$F_1$	0.87	11	0.59	27	0.84	11	0.52	27
	$F_2$	0.73	26	0.69	18	0.76	26	0.58	1
	$F_3$	0.72	27	0.72	11	0.72	27	0.58	1
	$E_1$	0.88	9	0.76	6	0.78	20	0.56	11
	$E_2$	0.84	13	0.76	6	0.78	20	0.56	11
	$E_3$	0.88	9	0.75	8	0.78	20	0.56	11
	$E_4$	0.87	11	0.75	8	0.83	12	0.56	11
	$N_1$	0.77	21	0.74	10	0.77	23	0.57	7
$K=1+\lfloor {\log }_{2}n\rfloor$	$N_2$	0.83	19	0.68	20	0.77	23	0.57	7
	$N_3$	0.80	20	0.69	18	0.77	23	0.57	7
	$N_4$	0.84	13	0.68	20	0.82	16	0.58	1
	$M_1$	0.92	7	0.68	20	0.83	12	0.55	21
	$M_2$	0.95	6	0.71	12	0.85	7	0.56	11
	$M_3$	0.92	7	0.68	20	0.83	12	0.56	11
	$M_4$	1.00	1	0.71	12	0.93	5	0.56	11
	$E_1$	0.84	13	0.85	1	0.85	7	0.54	26
	$E_2$	0.84	13	0.85	1	0.85	7	0.55	21
	$E_3$	0.84	13	0.85	1	0.85	7	0.55	21
	$E_4$	0.84	13	0.85	1	0.86	6	0.55	21
	$N_1$	0.76	23	0.71	12	0.82	16	0.58	1
$K=1+\lfloor \sqrt{n}\rfloor$	$N_2$	0.76	23	0.68	20	0.83	12	0.58	1
	$N_3$	0.76	23	0.66	26	0.82	16	0.58	1
	$N_4$	0.77	21	0.67	25	0.82	16	0.57	7
	$M_1$	0.98	3	0.70	17	0.94	1	0.56	11
	$M_2$	0.98	3	0.71	12	0.94	1	0.56	11
	$M_3$	0.98	3	0.71	12	0.94	1	0.56	11
	$M_4$	1.00	1	0.78	5	0.94	1	0.55	21

Table 3. Table 3: Accuracies of methods E 4 subscript 𝐸 4 E_{4} and M 4 subscript 𝑀 4 M_{4} on selected problems.

	method	$A_7$	$A_8$	$A_9$	$B_1$	$B_2$	$B_3$	${\left(I-\Lambda \right)}_1$	${\left(I-\Lambda \right)}_2$
$K=1+\lfloor {\log }_{2}n\rfloor$	$E_4$	1	0.54	0.35	0.38	0.9	0.96	0.52	0.24
	$M_4$	1	1	1	0.91	0.47	0.75	0.81	0.81
$K=1+\lfloor \sqrt{n}\rfloor$	$E_4$	0.68	0.58	0.35	0.85	0.9	0.8	0.57	0.51
	$M_4$	1	1	1	0.95	0.47	0.93	0.81	0.81

Equations

$$L_N=D^{-1/2}W{D}^{-1/2}.$$

$$y_{i,j}=u_{i,j}\big{/}\big{(}\sum_{r=1,k}u_{i,r}^{\,2}\big{)}^{1/2},\quad\hbox{for}\quad i=1,\ldots,n,\ j=1,\ldots,k.$$

$${\rm NMI}=\frac{2\,\sum\limits_{i=1}^{nc}\sum\limits_{j=1}^{k}n_{i,j}\log\big{(}\dfrac{n\ n_{i,j}}{n_{i}\ \widehat{n}_{j}}\big{)}}{\sum\limits_{j=1}^{k}\widehat{n}_{j}\log\dfrac{\widehat{n}_{j}}{n}+\sum\limits_{i=1}^{nc}n_{i}\log\dfrac{n_{i}}{n}}.$$

$$w_{i,j}=1for(i,j)\in {\text{mit}\Delta }_e,$$

$$w_{i,j}=\exp\Big{(}-\dfrac{\delta_{i,j}^{2}}{2\sigma^{2}}\Big{)}\ \hbox{for}\ (i,j)\in{\mit\Delta}_{e},$$

$$w_{i,j}=\exp\Big{(}-\dfrac{\delta_{i,j}^{2}}{2\sigma_{i}\sigma_{j}}\Big{)}\ \hbox{for}\ (i,j)\in{\mit\Delta}_{e},$$

$$\begin{array}[]{l}t=\max\limits_{\textstyle{(i,j)\in\hbox{MST}{(\mit\Delta)}_{e}}}\ \delta_{i,j},\\[11.38092pt] s_{i}=\max\limits_{\textstyle{(i,j)\in{\mit\Delta}_{e}}}\ \delta_{i,j},\ \hbox{for}\ i=1,\ldots,n,\\[11.38092pt] s={\rm mean\ of\ the\ }s_{i}.\end{array}$$

$${\tilde{\theta }}_r=\frac{1}{\mu }\sum _{s=1}^{\mu }{\theta }_{r,s}and{\tilde{\alpha }}_r=\frac{1}{\mu }\sum _{s=1}^{\mu }{\alpha }_{r,s}$$

Full text

Keywords: Spectral clustering, Similarity matrix, Minimum spanning tree

1 Introduction

Clustering, a key step for many data mining problems, can be applied to a variety of different kinds of documents as long as a distance measure can be assigned to define the similarity among the data objects. A clustering technique classifies the data into groups, called clusters, in such a way that the objects belonging to a same cluster are more similar to each other than to the objects belonging to different clusters. Clearly, a specific formulation of the clustering problem depends considerably on the metric that is assigned.

We assume here that the data objects are the $n$ distinct points {\mbox{\boldmathx}}_{1},\ldots,{\mbox{\boldmathx}}_{n} of a subset ${\cal X}$ of ${\bf R}^{m}$. It is standard to measure the distances in ${\bf R}^{m}$ by using norms. Both the 1-norm, known as Manhattan, and the 2-norm, known as Euclidean, are suitable, while the $\infty$-norm is not (see [1], Sec. 3.2.1.4). We use the Euclidean distance d({\mbox{\boldmathx}}_{i},{\mbox{\boldmathx}}_{j})=\|{\mbox{\boldmathx}}_{i}-{\mbox{\boldmathx}}_{j}\|_{2}. The given data objects are represented in a natural way by an undirected weighted complete graph $\mit\Gamma$, having $n$ nodes and the distances d({\mbox{\boldmathx}}_{i},{\mbox{\boldmathx}}_{j}) as weights on the edges. We assume also that the required number $k$ of clusters is a priori fixed, with $1<k\ll n$. Thus, a clustering ${{\mit\Pi}=\{\pi_{1},\ldots,\pi_{k}\}}$ corresponds to a partitioning of the indices $\{1,\ldots,n\}$.

To solve the clustering problem many different methods which work directly on the space ${\cal X}$ have been devised, but they are not entirely satisfying. For example, the widely used Euclidean-based $k$-means is very sensitive to the initialization and, applied directly to the original data, cannot deal with clusters that are nonlinearly separable in ${\cal X}$, has difficulties in taking into account other features like the cardinality or the density of each cluster, and in some situations can be very slow [2]. For this reason another approach can be followed using a graph-based procedure, which performs the clustering by constructing the so called similarity graph $G$, having the same nodes of $\mit\Gamma$ and possibly a subset of the edges of $\mit\Gamma$. The weights on the edges of $G$ are defined by using a similarity function to model the local neighborhood relationships between the data points, i.e. a higher similarity is assigned to pairs of closely related objects than to objects which are only weakly related. The $n\times n$ nonnegative symmetric adjacency matrix $W$ associated to $G$ is the so-called similarity matrix. By using a spectral clustering algorithm the subspace ${\cal X}$, where the original data live, is mapped to a new subspace ${\cal Y}$. The transformed representation gives better results than the original one, because the corresponding procedures are more flexible, can adjust to varying local densities and discover arbitrarily shaped clusters ([1], Sec. 6.7.1). Then the clustering in the subspace ${\cal Y}$ is computed and mapped back to ${\cal X}$. Spectral clustering algorithms have become very popular (see [3] for a thorough survey). In the experiments we choose the algorithm described in [4].

The similarity graph $G$ changes depending whether only distances among points are accounted for or density considerations are applied (see the literature review in [5] for a list of the different approaches). In general the procedures which exploit specifically the density are more expensive.

A similarity function widely used in literature to construct the similarity matrix is the Gaussian one. It depends on a scale parameter $\sigma$ whose choice can be crucial. A possible way to select $\sigma$ suggests to run the spectral algorithm repeatedly for different values of $\sigma$ and select the one that provides the best clustering according to some criterium. Besides the obvious increase in the computational cost, the drawback of this procedure is the difficulty of choosing a reliable quality measure.

In this paper we take into consideration procedures for setting suitable values to $\sigma$ without performing multiple runs. These procedures use techniques based either directly on a subgraph $\mit\Delta$ of $\mit\Gamma$ associated to the dataset or on the minimal spanning tree (MST) of $\mit\Delta$. We focus on aspects like the sparsity level of the similarity matrix and the computational cost of its construction. The paper is organized as follows: after a very concise review of the spectral clustering method, of the performance indices and of the minimum spanning tree in Section 2, in Section 3 the different types of sparsity of the similarity matrices, the role of a scale parameter in their construction and the different resulting methods we take into consideration are outlined. Finally, in Section 4 the results of the experimentation on both artificial and real-world datasets are presented and discussed.

2 Preliminaries

The idea behind a graph-based procedure is to convert the representation of $n$ objects given in the original space ${\cal X}$ into a transformed representation in a new $k$-dimensional feature space ${\cal Y}$, $k$ being the required number of clusters, in such a way that the similarity structure of the data is preserved. To this aim, the clustering problem is reformulated making use of the similarity graph $G$. The weight on the edge connecting the $i$th and the $j$th nodes depends mainly on their distance d({\mbox{\boldmathx}}_{i},{\mbox{\boldmathx}}_{j}) and is large for close points, small for far-away points. Other aspects, like for example the nearby density, could be taken into consideration. For a correct application of the spectral algorithm, the graph $G$ should be undirected and connected [6]. Section 3 is dedicated to the description of the similarity graphs we consider.

2.1 The normalized spectral clustering algorithm

Once the $n\times n$ similarity matrix $W$ associated to graph $G$ has been constructed, the spectral clustering can be performed. Spectral clustering algorithms are based on the eigen-decomposition of Laplacian matrices. A detailed description of many such algorithms can be found in [3].

The unnormalized Laplacian matrix is $L=D-W$, where $D$ is the degree matrix of $W$, i.e. the diagonal matrix whose principal entry is given by $d_{i,i}=\sum_{j=1}^{n}w_{i,j}$. Matrix $D$ can be seen as a local average of the similarity and is used as a normalization factor for $W$. An exposition of the properties of matrix $L$ and of its normalized versions, on which spectral algorithms rely, can be found in the tutorial [6]. In our experimentation we use the algorithm proposed by [4] which exploits the normalized Laplacian

[TABLE]

This algorithm runs as follows:

construct matrix $L_{N}$, 2. 2.

compute the matrix U=[{\mbox{\boldmathu}}_{1}, {\mbox{\boldmathu}}_{2},…, {\mbox{\boldmathu}}_{k}], where the {\mbox{\boldmathu}}_{i} are the mutually orthogonal eigenvectors of $L_{N}$ which correspond to the largest $k$ eigenvalues (ordered downward). 3. 3.

perform a row normalization on $U$ to map the points {\mbox{\boldmathu}}_{i} to a unit hypersphere, by setting

[TABLE]

For $i=1,\ldots,n$, points {\mbox{\boldmathy}}_{i} of components $y_{i,j}$, $j=1,\ldots,k$, are obtained. The new feature space ${\cal Y}$, which is a subspace of ${\bf R}^{k}$, is so identified. 4. 4.

compute the clustering of the points of ${\cal Y}$, which is mapped back to the original points of ${\cal X}$.

The graph-based approach allows discovering clusters of arbitrary shape. This property is not specific of the particular clustering method used in the final phase 4. Since the transformed representation enhances the clustering properties of the data, even the simple $k$-means can be successfully applied. In our experimentation we apply the clustering algorithm described in [7], which is based on a sequence of rotations and results to be fast and efficient.

Since the spectral clustering algorithm performs the dimensional reduction through the eigenvectors of $L_{N}$, the computation of the $k$ eigenvectors (using Lanczos algorithm) can have substantial computational complexity for large $n$ and $k$. The sparsity of matrix $W$ is then an important issue from this point of view.

2.2 The performance index

To measure the clustering validity one can rely on both internal or external indices (see [8] and its references). Internal indices generally evaluate the clustering in terms of compactness and separability, measured through Euclidean distances. When the data belong to $\mathbf{R}^{2}$ and $\mathbf{R}^{3}$ and the clusters are not convex or are badly separated or have different densities, these indices might not comply with what the human eye would suggest.

External indices evaluate how well the obtained clustering ${\mit\Pi}_{obt}$ matches the target, i.e. the assumed clustering ${\mit\Pi}_{tar}$ of the data. The results obtained in the experimentation of Section 4 have been evaluated by means of the external Normalized Mutual Information index (NMI) [9]. Three other external indices have also been used for comparison purpose, namely the purity index [10], the Rand index [11] and the Clustering Error [12]. Since their grading on the set of all the methods that will be introduced in Section 3 was fully matching, in Section 4 we report only the clustering evaluations based on the NMI index.

Let ${{\mit\Pi}_{tar}=\{\widehat{\pi}_{1},\ldots,\widehat{\pi}_{k}\}}$ be the target clustering, ${{\mit\Pi}_{obt}=\{\pi_{1},\ldots,\pi_{nc}\}}$ the computed clustering of $nc$ components and define $n_{i}=\#\;\pi_{i}$, $\widehat{n}_{j}=\#\;\widehat{\pi}_{j}$ and $n_{i,j}=\#\,(\pi_{i}\bigcap\widehat{\pi}_{j})$, for $i=1,\ldots,nc$ and $j=1,\ldots,k$. The number $n_{i,j}$ counts how many points in the $i$th cluster of ${\mit\Pi}_{obt}$ belong also to the $j$th cluster of ${\mit\Pi}_{tar}$. The NMI is given by

[TABLE]

In the literature there is also a version where the normalization is performed with the geometrical mean of the entropies instead of the arithmetical one. In the case $nc<k$, the algorithm is considered failed.

The larger the value NMI, the better the clustering performance according to the assumed target. Perfectly matching clusterings have NMI=1, but lower indices can still be considered acceptable when the difference is due to few far-away points, possibly outliers. The treatment of outliers is not banal, as can be seen from the last three datasets in Figure 1, where it is not easy to decide whether the points of the two bridges should be considered outliers or belonging to the clusters of the other points (and which one) or even should form two separated clusters. For this reason, no specific treatment for the outliers detection is performed in the experiments.

2.3 Minimum spanning tree

Given a connected undirected weighted graph with $n$ nodes, its minimum spanning tree (MST) is a subgraph which has the same nodes, but connected by no more than $n-1$ edges totalling the minimal weighting. Classical Prim’s algorithm [13] computes the MST with a computational cost of order $O(n^{2})$ if the graph is dense and of order $O(e+n\log n)$, $e$ being the number of edges, if the graph is sparse, provided that adjacency lists are used for the graph representation and a Fibonacci heap is used as a working data structure [13].

Prim’s algorithm starts from an arbitrary root node, constructs a tree which grows until it spans all the nodes and returns the $n$-vector $p$ such that $p_{i}$ is the parent of the $i$th node. MST is unique if there are no ties in the pairwise distances. Since the graph is connected, also its MST is connected. From this property it follows that the graph obtained by removing all edges longer than the longest edge of its MST is still connected [6].

Clustering procedures which exploit the MST of a graph have been proposed [14, 15]. They act by removing from the graph the edges with the largest distances, until only $k$ disconnected components are obtained, avoiding to drop edges which lead to subsets of a single element considered outliers. Anyway, MST based clustering algorithms are known to be highly instable in presence of noise and/or outliers [16].

We are interested in MST not as a clustering procedure, but as a technique for extracting suitable values of the parameters that will be used in the construction of the similarity graphs, as described in the next section.

3 The methods

Sparsity is an important issue from the computational point of view: a sparse graph can be represented by adjacency lists which allow both a reduction of the required storage and the implementation of more efficient matrix-vector product algorithms. Hence, when $n$ is large, a sparse representation of the similarity among the data conveying sufficient information to get an acceptable clustering should be used.

The more immediate sparsity structures are obtained from distance-based graphs which take into account the pairwise distances \delta_{i,j}=d({\mbox{\boldmathx}}_{i},{\mbox{\boldmathx}}_{j}) of the $n$ points {\mbox{\boldmathx}}_{i} in ${\bf R}^{m}$. The sparse graph obtained by dropping from the complete graph ${\mit\Gamma}$ some selected edges is denoted by ${\mit\Delta}$ and the similarity graph, which has the same nodes and edges of ${\mit\Delta}$ and weights defined through a similarity function, is denoted by $G$. The two cases will be considered:

$\bullet$  the full case when ${\mit\Delta}={\mit\Gamma}$, i.e. no edge is dropped;

$\bullet$  the sparse case where ${\mit\Delta}$ is a proper subgraph of ${\mit\Gamma}$ obtained according to a chosen sparsity model.

The ratio $\nu/n^{2}$, where $\nu$ is the number of dropped edges from ${\mit\Gamma}$, measures the sparsity level. The greater $\nu$, the higher the sparsity.

3.1 The sparsity model

In this paper we consider the distance-based sparsity obtained through techniques like the ones outlined in [6] where the sparsity level is imposed through a sparsity parameter. The two following techniques are here considered.

(a) $\epsilon$-neighbor sparsity, where $\epsilon$ is a selected threshold. In the sparse graph ${\mit\Delta}$ the nodes $i$ and $j$ are connected only if their pairwise distance $\delta_{i,j}$ is not greater than $\epsilon$. As noted in [6], if $\epsilon$ is set equal to the largest weight of the MST of ${\mit\Gamma}$, then graph $G$ results connected, but this choice could give a too large $\epsilon$ if some very far-away outliers are present. For this reason smaller values should be chosen. We denote these methods as $E$-methods.

(b) $K$-nearest neighbor sparsity, where $K$ is the sparsity parameter. In the sparse directed graph node $i$ is connected to node $j$ if {\mbox{\boldmathx}}_{j} is among the $K$-nearest neighbors of {\mbox{\boldmathx}}_{i}. Since the nearest neighborhood relationship is not symmetric, a symmetrization step must be provided to get an undirected graph. The symmetrization can be performed by connecting nodes $i$ and $j$ with an undirected edge

(b1) if both points {\mbox{\boldmathx}}_{i} and {\mbox{\boldmathx}}_{j} are among the $K$-nearest neighbors of the other one. We denote these methods as $M$-methods (the $M$ stands for mutual).

(b2) if either point {\mbox{\boldmathx}}_{i} or point {\mbox{\boldmathx}}_{j} is among the $K$-nearest neighbors of the other one. We denote these methods as $N$-methods (the $N$ stands for non mutual).

Contrary to $\epsilon$-neighbor sparse graphs, $K$-nearest neighbor sparse graphs allow an a-priori control of the achievable sparsity level. Moreover, $K$-nearest neighbor sparse graphs allow an easier detection of outliers. In fact, if before symmetrization no directed edge exists from the node $j\neq i$ to the node $i$ for all $j$, the point {\mbox{\boldmathx}}_{i} could be recognized as an outlier.

If the similarity graph obtained applying techniques (a) and (b) results to be unconnected, because too many connections have been lost, an aggregation step must be performed to guarantee sufficient connection.

For what concerns the $E$-methods, in the experimentation we set $\epsilon$ equal to the mean of the distances of each point from its $K$th closest neighbor. In any case the chosen integer $K$ depends on $n$ but is much smaller than $n$. In [6] the value $K_{\ell}=1+\lfloor\log_{2}n\rfloor$ is suggested, on the basis of asymptotic connectivity considerations. In our experiments we consider also the value $K_{s}=1+\lfloor\sqrt{n}\rfloor$ which leads to a lower sparsity of the similarity graph without a significant increase of the computational cost for the similarity matrix construction (see Subsection 3.4).

3.2 The similarity function

Let ${\mit\Delta}_{e}$ denote the set of edges of the sparse graph ${\mit\Delta}$ obtained according one of the techniques described above. Obviously the edge $(i,i)\not\in{\mit\Delta}_{e}$ and $\delta_{i,j}\neq 0$ for $(i,j)\in{\mit\Delta}_{e}$. A simple definition of the similarity function would be $w_{i,j}=1/\delta_{i,j}$ for $(i,j)\in{\mit\Delta}_{e}$, but this function is not satisfactory in presence of very close points. Similarity functions considered in literature are:

the unit similarity function

[TABLE]

and the Gaussian similarity function

[TABLE]

where the scale parameter $\sigma$ controls the decay rate of the distances.

The values $w_{i,j}$ are the weights of the similarity graph $G$. The similarity matrix $W$ is the adjacency matrix of $G$, hence its entries are $w_{i,j}$ for $(i,j)\in{\mit\Delta}_{e}$ and zero otherwise.

When the Gaussian similarity function (2) is used, a correct choice of $\sigma$ may be critical for the efficiency of the spectral algorithm. A too small value of $\sigma$ would give weights $w_{i,j}$ very close to 0, so that all the points would appear equally far-away. On the contrary, a too large value of $\sigma$ would give weights $w_{i,j}$ very close to 1, so that all the points would appear equally close. In both cases it would be difficult to discriminate between close and distant points. Hence an intermediate value between the smallest and the largest $\delta_{i,j}$, with $i\neq j$, must be chosen.

Instead of a global $\sigma$, different local scales $\sigma_{i}$ can be used for the different objects {\mbox{\boldmathx}}_{i}. Then (2) is replaced by

[TABLE]

in order to tune pairwise distances according to the local statistics of the surrounding neighbors [17]. These local scales are suggested especially for high-dimensional problems with large behavior variations, where the choice (3) is trusted in better results than (2).

A first simple idea for determining a reasonable value of $\sigma$ takes into account the decay rate of the exponential function $f(x)=\exp(-x^{2}/(2\sigma^{2}))$ whose inflection point $x=\sigma$ discriminates between a first part of rapid decay and a second flatter part. This suggests to look for gaps in the curve of the distances $\delta_{i,j}$ sorted in descending order. If we can individuate a large gap, we can exploit it as a suitable value for $\sigma$.

Another idea takes into account the histogram of the distances $\delta_{i,j}$. If the data form clusters, the histogram is multi-modal and, if $\sigma$ is chosen around the first mode, the affinity values of the points which form a cluster can be expected to be significantly larger than others. This suggests to choose for $\sigma$ a value somewhat smaller than the first mode of the histogram. In [18] this technique is said to work well for spherical-like clusters.

Although simple, these ideas are difficult to implement, because of the difficulty of detecting the right gap or mode. We have verified that their applicability is in effect restricted to a small number of cases, and for this reason we have discarded them.

To find a suitable value for the scale parameter, it is suggested in literature to run the spectral clustering algorithm repeatedly for different values of $\sigma$ and select the one that provides the best clustering according to a chosen quality measure. For example, in [4] it is suggested to choose for $\sigma$ the value which gives the tightest clusters in ${\cal Y}$. To implement this technique, a decreasing sequence could be used, starting with a value of $\sigma$ smaller than $\max_{i,j}\delta_{i,j}$. However, the drawback of this procedure is the choice of the quality measure to use, which might be a non-monotone function of $\sigma$. In this case a reliable determination of an acceptable $\sigma$ could not be obtained using only a small number of tries. Moreover, a test based on an internal performance index, like for example the DB index [19], does not guarantee that the final clustering is consistent with what a human would set as a target.

For a single run of the spectral clustering algorithm, a careful a-priori selection of $\sigma$ should be considered and we suggest in the next subsection to exploit either the minimum spanning tree of $\mit\Delta$, denoted by MST$(\mit\Delta)$, or directly $\mit\Delta$. In the first case MST$(\mit\Delta)$ is assumed as a reliable representation of the data and $\sigma$ is set equal to its largest weight. In the second case, we suggest for local $\sigma_{i}$’s the distance from {\mbox{\boldmathx}}_{i} of its $\ell$th nearest neighbor, $\ell$ being a chosen index much smaller than $n$. In [17] the index $\ell=7$ is suggested, but this choice, independent from $n$, appears somewhat arbitrary. Our choice of $\sigma_{i}$’s will be described in the next Section. From local scales a single global $\sigma$ can be obtained by averaging the $\sigma_{i}$’s.

If the similarity functions (2) or (3) are used, the graph $G$ might result to be sparse in practice also when no particular sparsity structure is explicitly imposed on ${\mit\Gamma}$, depending on the magnitude of the scale parameters, due to a possible fast decay of the exponential function. So, when edges with negligible weights are dropped and the adjacency lists are used to represent $G$, for all computational purposes this graph can be considered sparse.

3.3 Definition of the methods

By the term “method” we intend the combination of a sparsity model and a similarity function (either unit, global Gaussian or local Gaussian). The different similarity matrices $W$ so constructed are then processed by the spectral algorithm.

$\bullet$ In the full case, no sparsity is imposed on ${\mit\Gamma}$, i.e. ${\mit\Delta}={\mit\Gamma}$. In the following we refer to these methods as $F$-methods. The unit similarity function is not used.

—  Method $F_{1}$ implements the similarity function (2). In [6] the largest weight of MST$({\mit\Gamma})$ is suggested as $\sigma$. In our experimentation this choice appeared often too large, matrix $W$ resulted nearly full and gave really poor clusterings. Hence if $\sigma$ results larger than the mean $\delta_{\rm mean}$ of all the pairwise distances ${\delta}_{i,j}$, we set $\sigma$ equal to $\delta_{\rm mean}$.

—  Method $F_{2}$ implements the similarity function (3). The $i$th local scale $\sigma_{i}$ is set equal to the distance from {\mbox{\boldmathx}}_{i} of its $K$th nearest neighbor, with $K=K_{\ell}$. The sparsity level obtained with this value of $K$ is too low in most cases of the datasets considered in the experimentation, so the value $K=K_{s}$, which would give lower sparsity levels, has not been used.

—  Method $F_{3}$ implements the similarity function (2) with $\sigma$ equal to the mean of the previous $\sigma_{i}$’s.

$\bullet$ In the sparse case, let ${\mit\Delta}_{e}$ and MST$(\mit\Delta)_{e}$ be the edge sets of $\mit\Delta$ and of MST$(\mit\Delta)$, respectively. We set

[TABLE]

The name of a method indicates the sparsity model with a subscript which indicates the similarity function: specifically, subscript 1 indicates similarity (1), subscript 2 indicates similarity (2) with $\sigma=t$, subscript 3 indicates similarity (3) with $\sigma_{i}=s_{i}$, subscript 4 indicates similarity (2) with $\sigma=s$. For example $E_{1}$ denotes the $E$-method applied with the unit similarity function and $M_{3}$ denotes the $M$-method applied with local Gaussian function (3). Each method is applied with both $K=K_{\ell}$ and $K=K_{s}$.

On the whole, the methods considered for the experiments are $q=27$.

3.4 Computational costs

We now consider the computational cost of the described methods, under the assumption that all the distances are available. Hence the cost for the construction of ${\mit\Gamma}$ is not accounted for. When the graph is sparse, adjacency lists are used for the graph representation. The cost to construct the similarity matrix $W$ is given by $\gamma=\gamma_{1}+\gamma_{2}+\gamma_{3}$, where

– $\gamma_{1}$ is the cost of sparsifying ${\mit\Gamma}$,

– $\gamma_{2}$ is the cost of computing the scale parameters,

– $\gamma_{3}$ is the cost of evaluating the similarity function.

For $F$-methods no sparsity is performed, hence $\gamma_{1}=0$ and $\gamma_{3}=O(n^{2})$ since the similarity function is evaluated on a nearly full matrix. For method $F_{1}$, $\gamma_{2}$ is of order $O(n^{2})$ as it follows from Subsection 2.3. Also for methods $F_{2}$ and $F_{3}$, $\gamma_{2}=O(n^{2})$, indeed, given an integer $K<n$, $O(n)$ computational steps suffice to select the $K$th nearest neighbor to a point and to identify the set of the $h$th nearest neighbors, for $h=1,\ldots,K$ (see [13]).

For the same reason, the detection of the edges of ${\mit\Delta}$ for $E$-methods has a cost $\gamma_{1}=O(n^{2})$. For method $E_{1}$ the costs $\gamma_{2}$ and $\gamma_{3}$ are zero. The number $\eta$ of edges of ${\mit\Delta}$ is not a-priori quantifiable, and $\eta=O(n^{2})$. For method $E_{2}$ $\gamma_{2}=O(\eta+n\log n)$ and $\gamma_{3}=O(\eta)$. For methods $E_{3}$ and $E_{4}$ both $\gamma_{2}$ and $\gamma_{3}$ are of order $O(\eta)$. In any case all the costs result to be of order $O(n^{2})$.

For $N$-methods and $M$-methods, the number of edges of ${\mit\Delta}$ is of order $O(K\,n)$ and their selection has a cost $\gamma_{1}$ of order $O(n^{2})$. When the unit similarity function is used, the costs $\gamma_{2}$ and $\gamma_{3}$ are zero. Otherwise, both $\gamma_{2}$ and $\gamma_{3}$ are of order $O(K\,n)$.

In brief, all the above described methods have costs $\gamma$ of order $O(n^{2})$. The methods which exploit the sparsity of the graphs have computational costs dominated by the cost $\gamma_{1}$ of the sparsifying step.

4 The experimentation

The experiments have been conducted on an Intel(R)Core(TM)i7-4770CPU @3.40GHz using double precision arithmetic.

4.1 The datasets

Datasets belonging to four sets $A$, $B$, $R$, $U$ have been considered. The target clustering ${\mit\Pi}_{tar}$ is known for every dataset.

A, B. Sets $A$ and $B$ contain 12 artificial datasets of ${\bf R}^{2}$ which can be found in the Spatial Data Mining Project of [20]. Their clusters have circular or noncircular shape, uniform or nonuniform density. In Figure 1 they are shown according to growing difficulty. Different gray levels evidence the assumed target clusters. The first 9 datasets belong to set $A$ and have outliers and well separated clusters; the last 3 datasets belong to set $B$, have bridges which connect different clusters and lend themselves to possible different targets. The datasets have from 76 to 289 points, normalized in such a way that the maximum distance between points of each figure is equal to 1.

R. Set $R$ contains 6 artificial datasets generated in ${\bf R}^{3}$ which consist of two interlacing rings with increasing data dispersion and have 900 points each (see three of them in Figure 2), and the 2 datasets generated in ${\bf R}^{8}$ with 200 points each which can be found in [5] under the name I-$\Lambda$.

U. Set $U$ contains 4 real-world datasets iris, wine, vote and seeds which have been taken from the UCI database [21].

Dataset iris is a typical test case for many statistical classification techniques in machine learning. It consists of samples from three species of iris, with four features measured for each sample: the length and the width of sepals and petals. Iris contains $n=150$ points in ${\bf R}^{4}$ and 3 clusters are expected. Actually, the dataset contains only two well separated clusters: one cluster contains the measurements of only one species, while the other cluster contains the measurements of the other two species and is not evidently separable without further supervised information.

Dataset wine consists of samples of chemical analysis of wines derived from three different cultivars. The analysis determined 13 attributes for each sample. Wine contains $n=178$ points in ${\bf R}^{13}$ and 3 clusters are expected.

Dataset vote consists of 435 samples. Each sample lists the votes (a vote can be ”yes”, ”no” or can be missing) given by one U.S. congressman (out of 168 republican and 267 democrat congressmen) on 16 different questions. Vote contains $n=435$ points in ${\bf R}^{16}$ and 2 clusters are expected.

Dataset seeds consists of samples measuring the geometrical properties of the kernels of three different varieties of wheat: Kama, Rosa and Canadian, randomly selected for the experiment. A soft X-ray technique was used to construct seven real-valued attributes. Seeds contains $n=210$ points in ${\bf R}^{7}$ and 3 clusters are expected.

4.2 Results

The $q$ methods listed in Subsection 3.3 have been applied to each dataset obtaining almost always clusterings ${\mit\Pi}_{obt}$ with the expected number of clusters. For $r=1,\ldots,q$ and $s=1,\ldots,\mu$, the $r$th method has been applied to the $s$th dataset, with $\mu=9$ for the set $A$, $\mu=3$ for the set $B$, $\mu=8$ for the set $R$, $\mu=4$ for the set $U$. The following quantities have been computed:

– the sparsity level $\theta_{r,s}$ of the similarity graph $G$,

– the accuracy $\alpha_{r,s}$, measured by the NMI index of the clustering ${\mit\Pi}_{obt}$.

For each set $A$, $B$, $R$ and $U$, the averaged sparsity level and the averaged accuracy of the $r$th method, $r=1,\ldots,q,$

[TABLE]

are computed varying the $\mu$ datasets in the set.

The sparsity level is measured with respect to the machine zero $eps=2^{-52}$, i.e. $\theta_{r,s}=\nu/n^{2}$, where $\nu$ is the number of elements in $W$ which are smaller than $eps$. When the Gaussian similarity function is used, the reference to $eps$ instead of the real zero makes significant sparsity levels even in the case of a (theoretical) dense graph, as can be seen in Table 1. Except for the $F$-methods, $\widetilde{\theta}_{r}$ vary negligibly when the methods belong to the same class, pointing out that once a sparsity model has been chosen, the number of edges in ${\mit\Delta}$ is almost unaffected by the choice of $\sigma$. The index $i$ stands for any $i=1,\dots,4$.

Let $\rho_{r}$ denote the rank of the $r$th method, i.e. the position of $\widetilde{\alpha}_{r}$, in the not increasing sequence \widetilde{\mbox{\boldmath\alpha}}=\{\widetilde{\alpha}_{1},\ldots,\widetilde{\alpha}_{q}\}. Methods with equal accuracy receive the same ranking number according to the standard competition “1224” policy (i.e. the ranking of each method is equal to 1 plus the number of the methods which have a larger accuracy). Table 2 shows the averaged accuracies $\widetilde{\alpha}_{r}$ and the corresponding ranking $\rho_{r}$. From Tables 1 and 2, it appears that there is not a strict relation between the sparsity level of a method and its accuracy, in the sense that the accuracy seems to depend on “which” more than on “how many” edges are retained in the graph ${\mit\Delta}$. For example on class $A$, the $M$-methods with both values of $K$ get good results while the $N$-methods with $K=K_{\ell}$ which have an intermediate sparsity level, obtain poor results.

Figures 3 and 4 plot the averaged accuracies shown in Table 2, for $E$, $N$ and $M$-methods, grouped according to the value of $K$. From these figures it appears evident that the choice of the similarity function influences more the averaged accuracy, when $K=K_{\ell}$.

It is worth noting that the averaged accuracy of all the different methods for problems of class $U$ is very low and the range of its values is very small. This class of problems is not worthwhile to analyze and compare the considered methods. Restricting ourselves to the classes of problems $A$, $B$ and $R$, from Table 2 it appears that the $F$-methods and the $N$-methods are outperformed by $M$-methods on problems $A$ and $R$ and by $E$-methods on problems $B$. For this reason we limit our analysis to $E$ and $M$-methods.

Since the averaged values of the accuracy are not sufficient to show the different features of these methods, we analyze the behavior of selected methods on selected problems. We choose a representative method for both $E$ and $M$-methods. More precisely, due to their good performance on average, the methods $E_{4}$ and $M_{4}$, for the two values $K_{\ell}$ and $K_{s}$ are selected.

In Table 3 we list the accuracies of the four methods under investigation on three problems of class $A$, namely $A_{7}$, $A_{8}$, $A_{9}$, on all problems of class $B$ and on the two $I-\Lambda$ problems of class $R$. On all the remaining problems of class $A$ and $R$ these four methods find the target clusters. By inspection of Table 3, it appears that the $M$-methods always outperform $E$-methods except for problems with bridges. For $E$-methods the choice of the radius of the neighbors can be crucial especially when there are clusters with different densities. In order to better understand the meaning of Table 3, see Figures 5, 6, 7 and 8 where critical clustering outcomes are plotted and the corresponding values of accuracy $\alpha$ are given.

Figure 5 shows the poor performance of method $E_{4}$ with $K=K_{s}$ on problems $A_{7}$, $A_{8}$ and $A_{9}$. While for problem $A_{7}$ the choice of a smaller radius allows finding the target clustering, for problems $A_{8}$ and $A_{9}$ also the choice $K=K_{\ell}$ produces analogous results, may be due to nonuniform densities. Figure 6 shows the clustering obtained by method $E_{4}$ (on the left) having irrelevant differences from the target one and the wrong clustering obtained by method $M_{4}$ (on the right) for problem $B_{2}$ which has short bridges. Similar results are obtained for $K=K_{s}$. Figures 7 and 8 show that the choice of $K$ is critical for both methods $E_{4}$ and $M_{4}$ on problem $B_{3}$ which has bridges and clusters of nonuniform densities. To further investigate the negative influence of the bridges, we have removed bridges from the problems of class $B$, obtaining the target clusterings with all the four methods under consideration. In this case the nonuniform densities of problems $B_{1}$ and $B_{3}$ do not affect the behavior of $E$-methods, due to the sufficiently large distances among the clusters.

4.3 Conclusions

In this paper several techniques to construct the similarity matrix have been considered. The corresponding methods have been tested on different datasets by applying a normalized spectral algorithm. The methods based on the dense graph and those obtained by exploiting the non mutual $K$-nearest neighbor sparsity are clearly outperformed by those based on the $\epsilon$-neighbor sparsity and mutual $K$-nearest neighbor sparsity. The performance of the more effective $E$ and $M$-methods depends on the characteristics of the problems. Among the different $M$-methods, the one with $K=K_{s}$ and Gaussian similarity function given in (2) with $\sigma=s$ given in (4) appears to be the most reliable, even if in few cases it is outperformed by the corresponding $E$-method.