2.3.Clustering#

Clustering ofunlabeled data can be performed with the modulesklearn.cluster.

Each clustering algorithm comes in two variants: a class, that implementsthefit method to learn the clusters on train data, and a function,that, given train data, returns an array of integer labels correspondingto the different clusters. For the class, the labels over the trainingdata can be found in thelabels_ attribute.

Input data

One important thing to note is that the algorithms implemented inthis module can take different kinds of matrix as input. All themethods accept standard data matrices of shape(n_samples,n_features).These can be obtained from the classes in thesklearn.feature_extractionmodule. ForAffinityPropagation,SpectralClusteringandDBSCAN one can also input similarity matrices of shape(n_samples,n_samples). These can be obtained from the functionsin thesklearn.metrics.pairwise module.

2.3.1.Overview of clustering methods#

../_images/sphx_glr_plot_cluster_comparison_001.png

A comparison of the clustering algorithms in scikit-learn#

Method name

Parameters

Scalability

Usecase

Geometry (metric used)

K-Means

number of clusters

Very largen_samples, mediumn_clusters withMiniBatch code

General-purpose, even cluster size, flat geometry,not too many clusters, inductive

Distances between points

Affinity propagation

damping, sample preference

Not scalable with n_samples

Many clusters, uneven cluster size, non-flat geometry, inductive

Graph distance (e.g. nearest-neighbor graph)

Mean-shift

bandwidth

Not scalable withn_samples

Many clusters, uneven cluster size, non-flat geometry, inductive

Distances between points

Spectral clustering

number of clusters

Mediumn_samples, smalln_clusters

Few clusters, even cluster size, non-flat geometry, transductive

Graph distance (e.g. nearest-neighbor graph)

Ward hierarchical clustering

number of clusters or distance threshold

Largen_samples andn_clusters

Many clusters, possibly connectivity constraints, transductive

Distances between points

Agglomerative clustering

number of clusters or distance threshold, linkage type, distance

Largen_samples andn_clusters

Many clusters, possibly connectivity constraints, non Euclideandistances, transductive

Any pairwise distance

DBSCAN

neighborhood size

Very largen_samples, mediumn_clusters

Non-flat geometry, uneven cluster sizes, outlier removal,transductive

Distances between nearest points

HDBSCAN

minimum cluster membership, minimum point neighbors

largen_samples, mediumn_clusters

Non-flat geometry, uneven cluster sizes, outlier removal,transductive, hierarchical, variable cluster density

Distances between nearest points

OPTICS

minimum cluster membership

Very largen_samples, largen_clusters

Non-flat geometry, uneven cluster sizes, variable cluster density,outlier removal, transductive

Distances between points

Gaussian mixtures

many

Not scalable

Flat geometry, good for density estimation, inductive

Mahalanobis distances to centers

BIRCH

branching factor, threshold, optional global clusterer.

Largen_clusters andn_samples

Large dataset, outlier removal, data reduction, inductive

Euclidean distance between points

Bisecting K-Means

number of clusters

Very largen_samples, mediumn_clusters

General-purpose, even cluster size, flat geometry,no empty clusters, inductive, hierarchical

Distances between points

Non-flat geometry clustering is useful when the clusters have a specificshape, i.e. a non-flat manifold, and the standard euclidean distance isnot the right metric. This case arises in the two top rows of the figureabove.

Gaussian mixture models, useful for clustering, are described inanother chapter of the documentation dedicated tomixture models. KMeans can be seen as a special case of Gaussian mixturemodel with equal covariance per component.

Transductive clustering methods (in contrast toinductive clustering methods) are not designed to be applied to new,unseen data.

Examples

2.3.2.K-means#

TheKMeans algorithm clusters data by trying to separate samples in ngroups of equal variance, minimizing a criterion known as theinertia orwithin-cluster sum-of-squares (see below). This algorithm requires the numberof clusters to be specified. It scales well to large numbers of samples and hasbeen used across a large range of application areas in many different fields.

The k-means algorithm divides a set of\(N\) samples\(X\) into\(K\) disjoint clusters\(C\), each described by the mean\(\mu_j\)of the samples in the cluster. The means are commonly called the cluster“centroids”; note that they are not, in general, points from\(X\),although they live in the same space.

The K-means algorithm aims to choose centroids that minimise theinertia,orwithin-cluster sum-of-squares criterion:

\[\sum_{i=0}^{n}\min_{\mu_j \in C}(||x_i - \mu_j||^2)\]

Inertia can be recognized as a measure of how internally coherent clusters are.It suffers from various drawbacks:

  • Inertia makes the assumption that clusters are convex and isotropic,which is not always the case. It responds poorly to elongated clusters,or manifolds with irregular shapes.

  • Inertia is not a normalized metric: we just know that lower values arebetter and zero is optimal. But in very high-dimensional spaces, Euclideandistances tend to become inflated(this is an instance of the so-called “curse of dimensionality”).Running a dimensionality reduction algorithm such asPrincipal component analysis (PCA) prior tok-means clustering can alleviate this problem and speed up thecomputations.

../_images/sphx_glr_plot_kmeans_assumptions_002.png

For more detailed descriptions of the issues shown above and how to address them,refer to the examplesDemonstration of k-means assumptionsandSelecting the number of clusters with silhouette analysis on KMeans clustering.

K-means is often referred to as Lloyd’s algorithm. In basic terms, thealgorithm has three steps. The first step chooses the initial centroids, withthe most basic method being to choose\(k\) samples from the dataset\(X\). After initialization, K-means consists of looping between thetwo other steps. The first step assigns each sample to its nearest centroid.The second step creates new centroids by taking the mean value of all of thesamples assigned to each previous centroid. The difference between the oldand the new centroids are computed and the algorithm repeats these last twosteps until this value is less than a threshold. In other words, it repeatsuntil the centroids do not move significantly.

../_images/sphx_glr_plot_kmeans_digits_001.png

K-means is equivalent to the expectation-maximization algorithmwith a small, all-equal, diagonal covariance matrix.

The algorithm can also be understood through the concept ofVoronoi diagrams. First the Voronoi diagram ofthe points is calculated using the current centroids. Each segment in theVoronoi diagram becomes a separate cluster. Secondly, the centroids are updatedto the mean of each segment. The algorithm then repeats this until a stoppingcriterion is fulfilled. Usually, the algorithm stops when the relative decreasein the objective function between iterations is less than the given tolerancevalue. This is not the case in this implementation: iteration stops whencentroids move less than the tolerance.

Given enough time, K-means will always converge, however this may be to a localminimum. This is highly dependent on the initialization of the centroids.As a result, the computation is often done several times, with differentinitializations of the centroids. One method to help address this issue is thek-means++ initialization scheme, which has been implemented in scikit-learn(use theinit='k-means++' parameter). This initializes the centroids to be(generally) distant from each other, leading to probably better results thanrandom initialization, as shown in the reference. For detailed examples ofcomparing different initialization schemes, refer toA demo of K-Means clustering on the handwritten digits data andEmpirical evaluation of the impact of k-means initialization.

K-means++ can also be called independently to select seeds for otherclustering algorithms, seesklearn.cluster.kmeans_plusplus for detailsand example usage.

The algorithm supports sample weights, which can be given by a parametersample_weight. This allows to assign more weight to some samples whencomputing cluster centers and values of inertia. For example, assigning aweight of 2 to a sample is equivalent to adding a duplicate of that sampleto the dataset\(X\).

Examples

2.3.2.1.Low-level parallelism#

KMeans benefits from OpenMP based parallelism through Cython. Smallchunks of data (256 samples) are processed in parallel, which in additionyields a low memory footprint. For more details on how to control the number ofthreads, please refer to ourParallelism notes.

Examples

References#

2.3.2.2.Mini Batch K-Means#

TheMiniBatchKMeans is a variant of theKMeans algorithmwhich uses mini-batches to reduce the computation time, while still attemptingto optimise the same objective function. Mini-batches are subsets of the inputdata, randomly sampled in each training iteration. These mini-batchesdrastically reduce the amount of computation required to converge to a localsolution. In contrast to other algorithms that reduce the convergence time ofk-means, mini-batch k-means produces results that are generally only slightlyworse than the standard algorithm.

The algorithm iterates between two major steps, similar to vanilla k-means.In the first step,\(b\) samples are drawn randomly from the dataset, to forma mini-batch. These are then assigned to the nearest centroid. In the secondstep, the centroids are updated. In contrast to k-means, this is done on aper-sample basis. For each sample in the mini-batch, the assigned centroidis updated by taking the streaming average of the sample and all previoussamples assigned to that centroid. This has the effect of decreasing therate of change for a centroid over time. These steps are performed untilconvergence or a predetermined number of iterations is reached.

MiniBatchKMeans converges faster thanKMeans, but the qualityof the results is reduced. In practice this difference in quality can be quitesmall, as shown in the example and cited reference.

../_images/sphx_glr_plot_mini_batch_kmeans_001.png

Examples

References#

2.3.3.Affinity Propagation#

AffinityPropagation creates clusters by sending messages betweenpairs of samples until convergence. A dataset is then described using a smallnumber of exemplars, which are identified as those most representative of othersamples. The messages sent between pairs represent the suitability for onesample to be the exemplar of the other, which is updated in response to thevalues from other pairs. This updating happens iteratively until convergence,at which point the final exemplars are chosen, and hence the final clusteringis given.

../_images/sphx_glr_plot_affinity_propagation_001.png

Affinity Propagation can be interesting as it chooses the number ofclusters based on the data provided. For this purpose, the two importantparameters are thepreference, which controls how many exemplars areused, and thedamping factor which damps the responsibility andavailability messages to avoid numerical oscillations when updating thesemessages.

The main drawback of Affinity Propagation is its complexity. Thealgorithm has a time complexity of the order\(O(N^2 T)\), where\(N\)is the number of samples and\(T\) is the number of iterations untilconvergence. Further, the memory complexity is of the order\(O(N^2)\) if a dense similarity matrix is used, but reducible if asparse similarity matrix is used. This makes Affinity Propagation mostappropriate for small to medium sized datasets.

Algorithm description#

The messages sent between points belong to one of two categories. The first isthe responsibility\(r(i, k)\), which is the accumulated evidence thatsample\(k\) should be the exemplar for sample\(i\). The second is theavailability\(a(i, k)\) which is the accumulated evidence that sample\(i\) should choose sample\(k\) to be its exemplar, and considers thevalues for all other samples that\(k\) should be an exemplar. In this way,exemplars are chosen by samples if they are (1) similar enough to many samplesand (2) chosen by many samples to be representative of themselves.

More formally, the responsibility of a sample\(k\) to be the exemplar ofsample\(i\) is given by:

\[r(i, k) \leftarrow s(i, k) - max [ a(i, k') + s(i, k') \forall k' \neq k ]\]

Where\(s(i, k)\) is the similarity between samples\(i\) and\(k\).The availability of sample\(k\) to be the exemplar of sample\(i\) isgiven by:

\[a(i, k) \leftarrow min [0, r(k, k) + \sum_{i'~s.t.~i' \notin \{i, k\}}{r(i',k)}]\]

To begin with, all values for\(r\) and\(a\) are set to zero, and thecalculation of each iterates until convergence. As discussed above, in order toavoid numerical oscillations when updating the messages, the damping factor\(\lambda\) is introduced to iteration process:

\[r_{t+1}(i, k) = \lambda\cdot r_{t}(i, k) + (1-\lambda)\cdot r_{t+1}(i, k)\]
\[a_{t+1}(i, k) = \lambda\cdot a_{t}(i, k) + (1-\lambda)\cdot a_{t+1}(i, k)\]

where\(t\) indicates the iteration times.

Examples

2.3.4.Mean Shift#

MeanShift clustering aims to discoverblobs in a smooth density ofsamples. It is a centroid based algorithm, which works by updating candidatesfor centroids to be the mean of the points within a given region. Thesecandidates are then filtered in a post-processing stage to eliminatenear-duplicates to form the final set of centroids.

Mathematical details#

The position of centroid candidates is iteratively adjusted using a techniquecalled hill climbing, which finds local maxima of the estimated probabilitydensity. Given a candidate centroid\(x\) for iteration\(t\), thecandidate is updated according to the following equation:

\[x^{t+1} = x^t + m(x^t)\]

Where\(m\) is themean shift vector that is computed for each centroidthat points towards a region of the maximum increase in the density of points.To compute\(m\) we define\(N(x)\) as the neighborhood of sampleswithin a given distance around\(x\). Then\(m\) is computed using thefollowing equation, effectively updating a centroid to be the mean of thesamples within its neighborhood:

\[m(x) = \frac{1}{|N(x)|} \sum_{x_j \in N(x)}x_j - x\]

In general, the equation for\(m\) depends on a kernel used for densityestimation. The generic formula is:

\[m(x) = \frac{\sum_{x_j \in N(x)}K(x_j - x)x_j}{\sum_{x_j \in N(x)}K(x_j -x)} - x\]

In our implementation,\(K(x)\) is equal to 1 if\(x\) is small enoughand is equal to 0 otherwise. Effectively\(K(y - x)\) indicates whether\(y\) is in the neighborhood of\(x\).

The algorithm automatically sets the number of clusters, instead of relying on aparameterbandwidth, which dictates the size of the region to search through.This parameter can be set manually, but can be estimated using the providedestimate_bandwidth function, which is called if the bandwidth is not set.

The algorithm is not highly scalable, as it requires multiple nearest neighborsearches during the execution of the algorithm. The algorithm is guaranteed toconverge, however the algorithm will stop iterating when the change in centroidsis small.

Labelling a new sample is performed by finding the nearest centroid for agiven sample.

../_images/sphx_glr_plot_mean_shift_001.png

Examples

References#

2.3.5.Spectral clustering#

SpectralClustering performs a low-dimension embedding of theaffinity matrix between samples, followed by clustering, e.g., by KMeans,of the components of the eigenvectors in the low dimensional space.It is especially computationally efficient if the affinity matrix is sparseand theamg solver is used for the eigenvalue problem (Note, theamg solverrequires that thepyamg module is installed.)

The present version of SpectralClustering requires the number of clustersto be specified in advance. It works well for a small number of clusters,but is not advised for many clusters.

For two clusters, SpectralClustering solves a convex relaxation of thenormalized cutsproblem on the similarity graph: cutting the graph in two so that the weight ofthe edges cut is small compared to the weights of the edges inside eachcluster. This criteria is especially interesting when working on images, wheregraph vertices are pixels, and weights of the edges of the similarity graph arecomputed using a function of a gradient of the image.

noisy_imgsegmented_img

Warning

Transforming distance to well-behaved similarities

Note that if the values of your similarity matrix are not welldistributed, e.g. with negative values or with a distance matrixrather than a similarity, the spectral problem will be singular andthe problem not solvable. In which case it is advised to apply atransformation to the entries of the matrix. For instance, in thecase of a signed distance matrix, is common to apply a heat kernel:

similarity=np.exp(-beta*distance/distance.std())

See the examples for such an application.

Examples

2.3.5.1.Different label assignment strategies#

Different label assignment strategies can be used, corresponding to theassign_labels parameter ofSpectralClustering."kmeans" strategy can match finer details, but can be unstable.In particular, unless you control therandom_state, it may not bereproducible from run-to-run, as it depends on random initialization.The alternative"discretize" strategy is 100% reproducible, but tendsto create parcels of fairly even and geometrical shape.The recently added"cluster_qr" option is a deterministic alternative thattends to create the visually best partitioning on the example applicationbelow.

assign_labels="kmeans"

assign_labels="discretize"

assign_labels="cluster_qr"

coin_kmeans

coin_discretize

coin_cluster_qr

References#

2.3.5.2.Spectral Clustering Graphs#

Spectral Clustering can also be used to partition graphs via their spectralembeddings. In this case, the affinity matrix is the adjacency matrix of thegraph, and SpectralClustering is initialized withaffinity='precomputed':

>>>fromsklearn.clusterimportSpectralClustering>>>sc=SpectralClustering(3,affinity='precomputed',n_init=100,...assign_labels='discretize')>>>sc.fit_predict(adjacency_matrix)
References#

2.3.6.Hierarchical clustering#

Hierarchical clustering is a general family of clustering algorithms thatbuild nested clusters by merging or splitting them successively. Thishierarchy of clusters is represented as a tree (or dendrogram). The root of thetree is the unique cluster that gathers all the samples, the leaves being theclusters with only one sample. See theWikipedia page for more details.

TheAgglomerativeClustering object performs a hierarchical clusteringusing a bottom up approach: each observation starts in its own cluster, andclusters are successively merged together. The linkage criteria determines themetric used for the merge strategy:

  • Ward minimizes the sum of squared differences within all clusters. It is avariance-minimizing approach and in this sense is similar to the k-meansobjective function but tackled with an agglomerative hierarchicalapproach.

  • Maximum orcomplete linkage minimizes the maximum distance betweenobservations of pairs of clusters.

  • Average linkage minimizes the average of the distances between allobservations of pairs of clusters.

  • Single linkage minimizes the distance between the closestobservations of pairs of clusters.

AgglomerativeClustering can also scale to large number of sampleswhen it is used jointly with a connectivity matrix, but is computationallyexpensive when no connectivity constraints are added between samples: itconsiders at each step all the possible merges.

FeatureAgglomeration

TheFeatureAgglomeration uses agglomerative clustering togroup together features that look very similar, thus decreasing thenumber of features. It is a dimensionality reduction tool, seeUnsupervised dimensionality reduction.

2.3.6.1.Different linkage type: Ward, complete, average, and single linkage#

AgglomerativeClustering supports Ward, single, average, and completelinkage strategies.

../_images/sphx_glr_plot_linkage_comparison_001.png

Agglomerative cluster has a “rich get richer” behavior that leads touneven cluster sizes. In this regard, single linkage is the worststrategy, and Ward gives the most regular sizes. However, the affinity(or distance used in clustering) cannot be varied with Ward, thus for nonEuclidean metrics, average linkage is a good alternative. Single linkage,while not robust to noisy data, can be computed very efficiently and cantherefore be useful to provide hierarchical clustering of larger datasets.Single linkage can also perform well on non-globular data.

Examples

2.3.6.2.Visualization of cluster hierarchy#

It’s possible to visualize the tree representing the hierarchical merging of clustersas a dendrogram. Visual inspection can often be useful for understanding the structureof the data, though more so in the case of small sample sizes.

../_images/sphx_glr_plot_agglomerative_dendrogram_001.png

Examples

2.3.6.3.Adding connectivity constraints#

An interesting aspect ofAgglomerativeClustering is thatconnectivity constraints can be added to this algorithm (only adjacentclusters can be merged together), through a connectivity matrix that definesfor each sample the neighboring samples following a given structure of thedata. For instance, in the swiss-roll example below, the connectivityconstraints forbid the merging of points that are not adjacent on the swissroll, and thus avoid forming clusters that extend across overlapping folds ofthe roll.

unstructuredstructured

These constraint are useful to impose a certain local structure, but theyalso make the algorithm faster, especially when the number of the samplesis high.

The connectivity constraints are imposed via an connectivity matrix: ascipy sparse matrix that has elements only at the intersection of a rowand a column with indices of the dataset that should be connected. Thismatrix can be constructed from a-priori information: for instance, youmay wish to cluster web pages by only merging pages with a link pointingfrom one to another. It can also be learned from the data, for instanceusingsklearn.neighbors.kneighbors_graph to restrictmerging to nearest neighbors as inthis example, orusingsklearn.feature_extraction.image.grid_to_graph toenable only merging of neighboring pixels on an image, as in thecoin example.

Warning

Connectivity constraints with single, average and complete linkage

Connectivity constraints and single, complete or average linkage can enhancethe ‘rich getting richer’ aspect of agglomerative clustering,particularly so if they are built withsklearn.neighbors.kneighbors_graph. In the limit of a smallnumber of clusters, they tend to give a few macroscopically occupiedclusters and almost empty ones. (see the discussion inAgglomerative clustering with and without structure).Single linkage is the most brittle linkage option with regard to this issue.

../_images/sphx_glr_plot_agglomerative_clustering_001.png../_images/sphx_glr_plot_agglomerative_clustering_002.png../_images/sphx_glr_plot_agglomerative_clustering_003.png../_images/sphx_glr_plot_agglomerative_clustering_004.png

Examples

2.3.6.4.Varying the metric#

Single, average and complete linkage can be used with a variety of distances (oraffinities), in particular Euclidean distance (l2), Manhattan distance(or Cityblock, orl1), cosine distance, or any precomputed affinitymatrix.

  • l1 distance is often good for sparse features, or sparse noise: i.e.many of the features are zero, as in text mining using occurrences ofrare words.

  • cosine distance is interesting because it is invariant to globalscalings of the signal.

The guidelines for choosing a metric is to use one that maximizes thedistance between samples in different classes, and minimizes that withineach class.

../_images/sphx_glr_plot_agglomerative_clustering_metrics_005.png../_images/sphx_glr_plot_agglomerative_clustering_metrics_006.png../_images/sphx_glr_plot_agglomerative_clustering_metrics_007.png

Examples

2.3.6.5.Bisecting K-Means#

TheBisectingKMeans is an iterative variant ofKMeans, usingdivisive hierarchical clustering. Instead of creating all centroids at once, centroidsare picked progressively based on a previous clustering: a cluster is split into twonew clusters repeatedly until the target number of clusters is reached.

BisectingKMeans is more efficient thanKMeans when the number ofclusters is large since it only works on a subset of the data at each bisectionwhileKMeans always works on the entire dataset.

AlthoughBisectingKMeans can’t benefit from the advantages of the"k-means++"initialization by design, it will still produce comparable results thanKMeans(init="k-means++") in terms of inertia at cheaper computational costs, and willlikely produce better results thanKMeans with a random initialization.

This variant is more efficient to agglomerative clustering if the number of clusters issmall compared to the number of data points.

This variant also does not produce empty clusters.

There exist two strategies for selecting the cluster to split:

  • bisecting_strategy="largest_cluster" selects the cluster having the most points

  • bisecting_strategy="biggest_inertia" selects the cluster with biggest inertia(cluster with biggest Sum of Squared Errors within)

Picking by largest amount of data points in most cases produces result asaccurate as picking by inertia and is faster (especially for larger amount of datapoints, where calculating error may be costly).

Picking by largest amount of data points will also likely produce clusters of similarsizes whileKMeans is known to produce clusters of different sizes.

Difference between Bisecting K-Means and regular K-Means can be seen on exampleBisecting K-Means and Regular K-Means Performance Comparison.While the regular K-Means algorithm tends to create non-related clusters,clusters from Bisecting K-Means are well ordered and create quite a visible hierarchy.

References#

2.3.7.DBSCAN#

TheDBSCAN algorithm views clusters as areas of high densityseparated by areas of low density. Due to this rather generic view, clustersfound by DBSCAN can be any shape, as opposed to k-means which assumes thatclusters are convex shaped. The central component to the DBSCAN is the conceptofcore samples, which are samples that are in areas of high density. Acluster is therefore a set of core samples, each close to each other(measured by some distance measure)and a set of non-core samples that are close to a core sample (but are notthemselves core samples). There are two parameters to the algorithm,min_samples andeps,which define formally what we mean when we saydense.Highermin_samples or lowerepsindicate higher density necessary to form a cluster.

More formally, we define a core sample as being a sample in the dataset suchthat there existmin_samples other samples within a distance ofeps, which are defined asneighbors of the core sample. This tellsus that the core sample is in a dense area of the vector space. A clusteris a set of core samples that can be built by recursively taking a coresample, finding all of its neighbors that are core samples, finding all oftheir neighbors that are core samples, and so on. A cluster also has aset of non-core samples, which are samples that are neighbors of a core samplein the cluster but are not themselves core samples. Intuitively, these samplesare on the fringes of a cluster.

Any core sample is part of a cluster, by definition. Any sample that is not acore sample, and is at leasteps in distance from any core sample, isconsidered an outlier by the algorithm.

While the parametermin_samples primarily controls how tolerant thealgorithm is towards noise (on noisy and large data sets it may be desirableto increase this parameter), the parametereps iscrucial to chooseappropriately for the data set and distance function and usually cannot beleft at the default value. It controls the local neighborhood of the points.When chosen too small, most data will not be clustered at all (and labeledas-1 for “noise”). When chosen too large, it causes close clusters tobe merged into one cluster, and eventually the entire data set to be returnedas a single cluster. Some heuristics for choosing this parameter have beendiscussed in the literature, for example based on a knee in the nearest neighbordistances plot (as discussed in the references below).

In the figure below, the color indicates cluster membership, with large circlesindicating core samples found by the algorithm. Smaller circles are non-coresamples that are still part of a cluster. Moreover, the outliers are indicatedby black points below.

dbscan_results

Examples

Implementation#

The DBSCAN algorithm is deterministic, always generating the same clusters whengiven the same data in the same order. However, the results can differ whendata is provided in a different order. First, even though the core samples willalways be assigned to the same clusters, the labels of those clusters willdepend on the order in which those samples are encountered in the data. Secondand more importantly, the clusters to which non-core samples are assigned candiffer depending on the data order. This would happen when a non-core samplehas a distance lower thaneps to two core samples in different clusters. Bythe triangular inequality, those two core samples must be more distant thaneps from each other, or they would be in the same cluster. The non-coresample is assigned to whichever cluster is generated first in a pass through thedata, and so the results will depend on the data ordering.

The current implementation uses ball trees and kd-trees to determine theneighborhood of points, which avoids calculating the full distance matrix (aswas done in scikit-learn versions before 0.14). The possibility to use custommetrics is retained; for details, seeNearestNeighbors.

Memory consumption for large sample sizes#

This implementation is by default not memory efficient because it constructs afull pairwise similarity matrix in the case where kd-trees or ball-trees cannotbe used (e.g., with sparse matrices). This matrix will consume\(n^2\)floats. A couple of mechanisms for getting around this are:

  • UseOPTICS clustering in conjunction with theextract_dbscanmethod. OPTICS clustering also calculates the full pairwise matrix, but onlykeeps one row in memory at a time (memory complexity n).

  • A sparse radius neighborhood graph (where missing entries are presumed to beout of eps) can be precomputed in a memory-efficient way and dbscan can be runover this withmetric='precomputed'. Seesklearn.neighbors.NearestNeighbors.radius_neighbors_graph.

  • The dataset can be compressed, either by removing exact duplicates if theseoccur in your data, or by using BIRCH. Then you only have a relatively smallnumber of representatives for a large number of points. You can then provide asample_weight when fitting DBSCAN.

References#

2.3.8.HDBSCAN#

TheHDBSCAN algorithm can be seen as an extension ofDBSCANandOPTICS. Specifically,DBSCAN assumes that the clusteringcriterion (i.e. density requirement) isglobally homogeneous.In other words,DBSCAN may struggle to successfully capture clusterswith different densities.HDBSCAN alleviates this assumption and explores all possible densityscales by building an alternative representation of the clustering problem.

Note

This implementation is adapted from the original implementation of HDBSCAN,scikit-learn-contrib/hdbscan based on[LJ2017].

Examples

2.3.8.1.Mutual Reachability Graph#

HDBSCAN first defines\(d_c(x_p)\), thecore distance of a sample\(x_p\), as thedistance to itsmin_samples th-nearest neighbor, counting itself. For example,ifmin_samples=5 and\(x_*\) is the 5th-nearest neighbor of\(x_p\)then the core distance is:

\[d_c(x_p)=d(x_p, x_*).\]

Next it defines\(d_m(x_p, x_q)\), themutual reachability distance of two points\(x_p, x_q\), as:

\[d_m(x_p, x_q) = \max\{d_c(x_p), d_c(x_q), d(x_p, x_q)\}\]

These two notions allow us to construct themutual reachability graph\(G_{ms}\) defined for a fixed choice ofmin_samples by associating eachsample\(x_p\) with a vertex of the graph, and thus edges between points\(x_p, x_q\) are the mutual reachability distance\(d_m(x_p, x_q)\)between them. We may build subsets of this graph, denoted as\(G_{ms,\varepsilon}\), by removing any edges with value greater than\(\varepsilon\):from the original graph. Any points whose core distance is less than\(\varepsilon\):are at this staged marked as noise. The remaining points are then clustered byfinding the connected components of this trimmed graph.

Note

Taking the connected components of a trimmed graph\(G_{ms,\varepsilon}\) isequivalent to running DBSCAN* withmin_samples and\(\varepsilon\). DBSCAN* is aslightly modified version of DBSCAN mentioned in[CM2013].

2.3.8.2.Hierarchical Clustering#

HDBSCAN can be seen as an algorithm which performs DBSCAN* clustering across allvalues of\(\varepsilon\). As mentioned prior, this is equivalent to finding the connectedcomponents of the mutual reachability graphs for all values of\(\varepsilon\). To do thisefficiently, HDBSCAN first extracts a minimum spanning tree (MST) from the fully-connected mutual reachability graph, then greedily cuts the edges with highestweight. An outline of the HDBSCAN algorithm is as follows:

  1. Extract the MST of\(G_{ms}\).

  2. Extend the MST by adding a “self edge” for each vertex, with weight equalto the core distance of the underlying sample.

  3. Initialize a single cluster and label for the MST.

  4. Remove the edge with the greatest weight from the MST (ties areremoved simultaneously).

  5. Assign cluster labels to the connected components which contain theend points of the now-removed edge. If the component does not have at leastone edge it is instead assigned a “null” label marking it as noise.

  6. Repeat 4-5 until there are no more connected components.

HDBSCAN is therefore able to obtain all possible partitions achievable byDBSCAN* for a fixed choice ofmin_samples in a hierarchical fashion.Indeed, this allows HDBSCAN to perform clustering across multiple densitiesand as such it no longer needs\(\varepsilon\) to be given as a hyperparameter. Insteadit relies solely on the choice ofmin_samples, which tends to be a more robusthyperparameter.

hdbscan_ground_truth

hdbscan_results

HDBSCAN can be smoothed with an additional hyperparametermin_cluster_sizewhich specifies that during the hierarchical clustering, components with fewerthanminimum_cluster_size many samples are considered noise. In practice, onecan setminimum_cluster_size=min_samples to couple the parameters andsimplify the hyperparameter space.

References

[CM2013]

Campello, R.J.G.B., Moulavi, D., Sander, J. (2013). Density-BasedClustering Based on Hierarchical Density Estimates. In: Pei, J., Tseng, V.S.,Cao, L., Motoda, H., Xu, G. (eds) Advances in Knowledge Discovery and DataMining. PAKDD 2013. Lecture Notes in Computer Science(), vol 7819. Springer,Berlin, Heidelberg.Density-Based Clustering Based on HierarchicalDensity Estimates

[LJ2017]

L. McInnes and J. Healy, (2017). Accelerated Hierarchical DensityBased Clustering. In: IEEE International Conference on Data Mining Workshops(ICDMW), 2017, pp. 33-42.Accelerated Hierarchical Density BasedClustering

2.3.9.OPTICS#

TheOPTICS algorithm shares many similarities with theDBSCANalgorithm, and can be considered a generalization of DBSCAN that relaxes theeps requirement from a single value to a value range. The key differencebetween DBSCAN and OPTICS is that the OPTICS algorithm builds areachabilitygraph, which assigns each sample both areachability_ distance, and a spotwithin the clusterordering_ attribute; these two attributes are assignedwhen the model is fitted, and are used to determine cluster membership. IfOPTICS is run with the default value ofinf set formax_eps, then DBSCANstyle cluster extraction can be performed repeatedly in linear time for anygiveneps value using thecluster_optics_dbscan method. Settingmax_eps to a lower value will result in shorter run times, and can bethought of as the maximum neighborhood radius from each point to find otherpotential reachable points.

optics_results

Thereachability distances generated by OPTICS allow for variable densityextraction of clusters within a single data set. As shown in the above plot,combiningreachability distances and data setordering_ produces areachability plot, where point density is represented on the Y-axis, andpoints are ordered such that nearby points are adjacent. ‘Cutting’ thereachability plot at a single value produces DBSCAN like results; all pointsabove the ‘cut’ are classified as noise, and each time that there is a breakwhen reading from left to right signifies a new cluster. The default clusterextraction with OPTICS looks at the steep slopes within the graph to findclusters, and the user can define what counts as a steep slope using theparameterxi. There are also other possibilities for analysis on the graphitself, such as generating hierarchical representations of the data throughreachability-plot dendrograms, and the hierarchy of clusters detected by thealgorithm can be accessed through thecluster_hierarchy_ parameter. Theplot above has been color-coded so that cluster colors in planar space matchthe linear segment clusters of the reachability plot. Note that the blue andred clusters are adjacent in the reachability plot, and can be hierarchicallyrepresented as children of a larger parent cluster.

Examples

Comparison with DBSCAN#

The results from OPTICScluster_optics_dbscan method and DBSCAN are verysimilar, but not always identical; specifically, labeling of periphery and noisepoints. This is in part because the first samples of each dense area processedby OPTICS have a large reachability value while being close to other points intheir area, and will thus sometimes be marked as noise rather than periphery.This affects adjacent points when they are considered as candidates for beingmarked as either periphery or noise.

Note that for any single value ofeps, DBSCAN will tend to have a shorterrun time than OPTICS; however, for repeated runs at varyingeps values, asingle run of OPTICS may require less cumulative runtime than DBSCAN. It is alsoimportant to note that OPTICS’ output is close to DBSCAN’s only ifeps andmax_eps are close.

Computational Complexity#

Spatial indexing trees are used to avoid calculating the full distance matrix,and allow for efficient memory usage on large sets of samples. Differentdistance metrics can be supplied via themetric keyword.

For large datasets, similar (but not identical) results can be obtained viaHDBSCAN. The HDBSCAN implementation is multithreaded, and has betteralgorithmic runtime complexity than OPTICS, at the cost of worse memory scaling.For extremely large datasets that exhaust system memory using HDBSCAN, OPTICSwill maintain\(n\) (as opposed to\(n^2\)) memory scaling; however,tuning of themax_eps parameter will likely need to be used to give asolution in a reasonable amount of wall time.

References#

  • “OPTICS: ordering points to identify the clustering structure.” Ankerst,Mihael, Markus M. Breunig, Hans-Peter Kriegel, and Jörg Sander. In ACM SigmodRecord, vol. 28, no. 2, pp. 49-60. ACM, 1999.

2.3.10.BIRCH#

TheBirch builds a tree called the Clustering Feature Tree (CFT)for the given data. The data is essentially lossy compressed to a set ofClustering Feature nodes (CF Nodes). The CF Nodes have a number ofsubclusters called Clustering Feature subclusters (CF Subclusters)and these CF Subclusters located in the non-terminal CF Nodescan have CF Nodes as children.

The CF Subclusters hold the necessary information for clustering which preventsthe need to hold the entire input data in memory. This information includes:

  • Number of samples in a subcluster.

  • Linear Sum - An n-dimensional vector holding the sum of all samples

  • Squared Sum - Sum of the squared L2 norm of all samples.

  • Centroids - To avoid recalculation linear sum / n_samples.

  • Squared norm of the centroids.

The BIRCH algorithm has two parameters, the threshold and the branching factor.The branching factor limits the number of subclusters in a node and thethreshold limits the distance between the entering sample and the existingsubclusters.

This algorithm can be viewed as an instance or data reduction method,since it reduces the input data to a set of subclusters which are obtained directlyfrom the leaves of the CFT. This reduced data can be further processed by feedingit into a global clusterer. This global clusterer can be set byn_clusters.Ifn_clusters is set to None, the subclusters from the leaves are directlyread off, otherwise a global clustering step labels these subclusters into globalclusters (labels) and the samples are mapped to the global label of the nearest subcluster.

Algorithm description#

  • A new sample is inserted into the root of the CF Tree which is a CF Node. Itis then merged with the subcluster of the root, that has the smallest radiusafter merging, constrained by the threshold and branching factor conditions.If the subcluster has any child node, then this is done repeatedly till itreaches a leaf. After finding the nearest subcluster in the leaf, theproperties of this subcluster and the parent subclusters are recursivelyupdated.

  • If the radius of the subcluster obtained by merging the new sample and thenearest subcluster is greater than the square of the threshold and if thenumber of subclusters is greater than the branching factor, then a space istemporarily allocated to this new sample. The two farthest subclusters aretaken and the subclusters are divided into two groups on the basis of thedistance between these subclusters.

  • If this split node has a parent subcluster and there is room for a newsubcluster, then the parent is split into two. If there is no room, then thisnode is again split into two and the process is continued recursively, till itreaches the root.

BIRCH or MiniBatchKMeans?#

  • BIRCH does not scale very well to high dimensional data. As a rule of thumb ifn_features is greater than twenty, it is generally better to use MiniBatchKMeans.

  • If the number of instances of data needs to be reduced, or if one wants alarge number of subclusters either as a preprocessing step or otherwise,BIRCH is more useful than MiniBatchKMeans.

../_images/sphx_glr_plot_birch_vs_minibatchkmeans_001.png
How to use partial_fit?#

To avoid the computation of global clustering, for every call ofpartial_fitthe user is advised:

  1. To setn_clusters=None initially.

  2. Train all data by multiple calls to partial_fit.

  3. Setn_clusters to a required value usingbrc.set_params(n_clusters=n_clusters).

  4. Callpartial_fit finally with no arguments, i.e.brc.partial_fit()which performs the global clustering.

References#

2.3.11.Clustering performance evaluation#

Evaluating the performance of a clustering algorithm is not as trivial ascounting the number of errors or the precision and recall of a supervisedclassification algorithm. In particular any evaluation metric should nottake the absolute values of the cluster labels into account but ratherif this clustering define separations of the data similar to some groundtruth set of classes or satisfying some assumption such that membersbelong to the same class are more similar than members of differentclasses according to some similarity metric.

2.3.11.1.Rand index#

Given the knowledge of the ground truth class assignmentslabels_true and our clustering algorithm assignments of the samesampleslabels_pred, the(adjusted or unadjusted) Rand indexis a function that measures thesimilarity of the two assignments,ignoring permutations:

>>>fromsklearnimportmetrics>>>labels_true=[0,0,0,1,1,1]>>>labels_pred=[0,0,1,1,2,2]>>>metrics.rand_score(labels_true,labels_pred)0.66

The Rand index does not ensure to obtain a value close to 0.0 for arandom labelling. The adjusted Rand indexcorrects for chance andwill give such a baseline.

>>>metrics.adjusted_rand_score(labels_true,labels_pred)0.24

As with all clustering metrics, one can permute 0 and 1 in the predictedlabels, rename 2 to 3, and get the same score:

>>>labels_pred=[1,1,0,0,3,3]>>>metrics.rand_score(labels_true,labels_pred)0.66>>>metrics.adjusted_rand_score(labels_true,labels_pred)0.24

Furthermore, bothrand_score andadjusted_rand_score aresymmetric: swapping the argument does not change the scores. They canthus be used asconsensus measures:

>>>metrics.rand_score(labels_pred,labels_true)0.66>>>metrics.adjusted_rand_score(labels_pred,labels_true)0.24

Perfect labeling is scored 1.0:

>>>labels_pred=labels_true[:]>>>metrics.rand_score(labels_true,labels_pred)1.0>>>metrics.adjusted_rand_score(labels_true,labels_pred)1.0

Poorly agreeing labels (e.g. independent labelings) have lower scores,and for the adjusted Rand index the score will be negative or close tozero. However, for the unadjusted Rand index the score, while lower,will not necessarily be close to zero:

>>>labels_true=[0,0,0,0,0,0,1,1]>>>labels_pred=[0,1,2,3,4,5,5,6]>>>metrics.rand_score(labels_true,labels_pred)0.39>>>metrics.adjusted_rand_score(labels_true,labels_pred)-0.072

Advantages:

  • Interpretability: The unadjusted Rand index is proportional to thenumber of sample pairs whose labels are the same in bothlabels_pred andlabels_true, or are different in both.

  • Random (uniform) label assignments have an adjusted Rand index score closeto 0.0 for any value ofn_clusters andn_samples (which is not thecase for the unadjusted Rand index or the V-measure for instance).

  • Bounded range: Lower values indicate different labelings, similarclusterings have a high (adjusted or unadjusted) Rand index, 1.0 is theperfect match score. The score range is [0, 1] for the unadjusted Rand indexand [-0.5, 1] for the adjusted Rand index.

  • No assumption is made on the cluster structure: The (adjusted orunadjusted) Rand index can be used to compare all kinds of clusteringalgorithms, and can be used to compare clustering algorithms such as k-meanswhich assumes isotropic blob shapes with results of spectral clusteringalgorithms which can find cluster with “folded” shapes.

Drawbacks:

  • Contrary to inertia, the(adjusted or unadjusted) Rand index requiresknowledge of the ground truth classes which is almost never available inpractice or requires manual assignment by human annotators (as in thesupervised learning setting).

    However (adjusted or unadjusted) Rand index can also be useful in a purelyunsupervised setting as a building block for a Consensus Index that can beused for clustering model selection (TODO).

  • Theunadjusted Rand index is often close to 1.0 even if the clusteringsthemselves differ significantly. This can be understood when interpretingthe Rand index as the accuracy of element pair labeling resulting from theclusterings: In practice there often is a majority of element pairs that areassigned thedifferent pair label under both the predicted and theground truth clustering resulting in a high proportion of pair labels thatagree, which leads subsequently to a high score.

Examples

Mathematical formulation#

If C is a ground truth class assignment and K the clustering, let us define\(a\) and\(b\) as:

  • \(a\), the number of pairs of elements that are in the same set in C andin the same set in K

  • \(b\), the number of pairs of elements that are in different sets in C andin different sets in K

The unadjusted Rand index is then given by:

\[\text{RI} = \frac{a + b}{C_2^{n_{samples}}}\]

where\(C_2^{n_{samples}}\) is the total number of possible pairs in thedataset. It does not matter if the calculation is performed on ordered pairs orunordered pairs as long as the calculation is performed consistently.

However, the Rand index does not guarantee that random label assignments willget a value close to zero (esp. if the number of clusters is in the same orderof magnitude as the number of samples).

To counter this effect we can discount the expected RI\(E[\text{RI}]\) ofrandom labelings by defining the adjusted Rand index as follows:

\[\text{ARI} = \frac{\text{RI} - E[\text{RI}]}{\max(\text{RI}) - E[\text{RI}]}\]
References#

2.3.11.2.Mutual Information based scores#

Given the knowledge of the ground truth class assignmentslabels_true andour clustering algorithm assignments of the same sampleslabels_pred, theMutual Information is a function that measures theagreement of the twoassignments, ignoring permutations. Two different normalized versions of thismeasure are available,Normalized Mutual Information (NMI) andAdjustedMutual Information (AMI). NMI is often used in the literature, while AMI wasproposed more recently and isnormalized against chance:

>>>fromsklearnimportmetrics>>>labels_true=[0,0,0,1,1,1]>>>labels_pred=[0,0,1,1,2,2]>>>metrics.adjusted_mutual_info_score(labels_true,labels_pred)0.22504

One can permute 0 and 1 in the predicted labels, rename 2 to 3 and getthe same score:

>>>labels_pred=[1,1,0,0,3,3]>>>metrics.adjusted_mutual_info_score(labels_true,labels_pred)0.22504

All,mutual_info_score,adjusted_mutual_info_score andnormalized_mutual_info_score are symmetric: swapping the argument doesnot change the score. Thus they can be used as aconsensus measure:

>>>metrics.adjusted_mutual_info_score(labels_pred,labels_true)0.22504

Perfect labeling is scored 1.0:

>>>labels_pred=labels_true[:]>>>metrics.adjusted_mutual_info_score(labels_true,labels_pred)1.0>>>metrics.normalized_mutual_info_score(labels_true,labels_pred)1.0

This is not true formutual_info_score, which is therefore harder to judge:

>>>metrics.mutual_info_score(labels_true,labels_pred)0.69

Bad (e.g. independent labelings) have non-positive scores:

>>>labels_true=[0,1,2,0,3,4,5,1]>>>labels_pred=[1,1,0,0,2,2,2,2]>>>metrics.adjusted_mutual_info_score(labels_true,labels_pred)-0.10526

Advantages:

  • Random (uniform) label assignments have a AMI score close to 0.0 for anyvalue ofn_clusters andn_samples (which is not the case for rawMutual Information or the V-measure for instance).

  • Upper bound of 1: Values close to zero indicate two label assignmentsthat are largely independent, while values close to one indicate significantagreement. Further, an AMI of exactly 1 indicates that the two labelassignments are equal (with or without permutation).

Drawbacks:

  • Contrary to inertia,MI-based measures require the knowledge of the groundtruth classes while almost never available in practice or requires manualassignment by human annotators (as in the supervised learning setting).

    However MI-based measures can also be useful in purely unsupervised settingas a building block for a Consensus Index that can be used for clusteringmodel selection.

  • NMI and MI are not adjusted against chance.

Examples

Mathematical formulation#

Assume two label assignments (of the same N objects),\(U\) and\(V\).Their entropy is the amount of uncertainty for a partition set, defined by:

\[H(U) = - \sum_{i=1}^{|U|}P(i)\log(P(i))\]

where\(P(i) = |U_i| / N\) is the probability that an object picked atrandom from\(U\) falls into class\(U_i\). Likewise for\(V\):

\[H(V) = - \sum_{j=1}^{|V|}P'(j)\log(P'(j))\]

With\(P'(j) = |V_j| / N\). The mutual information (MI) between\(U\)and\(V\) is calculated by:

\[\text{MI}(U, V) = \sum_{i=1}^{|U|}\sum_{j=1}^{|V|}P(i, j)\log\left(\frac{P(i,j)}{P(i)P'(j)}\right)\]

where\(P(i, j) = |U_i \cap V_j| / N\) is the probability that an objectpicked at random falls into both classes\(U_i\) and\(V_j\).

It also can be expressed in set cardinality formulation:

\[\text{MI}(U, V) = \sum_{i=1}^{|U|} \sum_{j=1}^{|V|} \frac{|U_i \cap V_j|}{N}\log\left(\frac{N|U_i \cap V_j|}{|U_i||V_j|}\right)\]

The normalized mutual information is defined as

\[\text{NMI}(U, V) = \frac{\text{MI}(U, V)}{\text{mean}(H(U), H(V))}\]

This value of the mutual information and also the normalized variant is notadjusted for chance and will tend to increase as the number of different labels(clusters) increases, regardless of the actual amount of “mutual information”between the label assignments.

The expected value for the mutual information can be calculated using thefollowing equation[VEB2009]. In this equation,\(a_i = |U_i|\) (the numberof elements in\(U_i\)) and\(b_j = |V_j|\) (the number of elements in\(V_j\)).

\[E[\text{MI}(U,V)]=\sum_{i=1}^{|U|} \sum_{j=1}^{|V|} \sum_{n_{ij}=(a_i+b_j-N)^+}^{\min(a_i, b_j)} \frac{n_{ij}}{N}\log \left( \frac{ N.n_{ij}}{a_i b_j}\right)\frac{a_i!b_j!(N-a_i)!(N-b_j)!}{N!n_{ij}!(a_i-n_{ij})!(b_j-n_{ij})!(N-a_i-b_j+n_{ij})!}\]

Using the expected value, the adjusted mutual information can then be calculatedusing a similar form to that of the adjusted Rand index:

\[\text{AMI} = \frac{\text{MI} - E[\text{MI}]}{\text{mean}(H(U), H(V)) - E[\text{MI}]}\]

For normalized mutual information and adjusted mutual information, thenormalizing value is typically somegeneralized mean of the entropies of eachclustering. Various generalized means exist, and no firm rules exist forpreferring one over the others. The decision is largely a field-by-field basis;for instance, in community detection, the arithmetic mean is most common. Eachnormalizing method provides “qualitatively similar behaviours”[YAT2016]. Inour implementation, this is controlled by theaverage_method parameter.

Vinh et al. (2010) named variants of NMI and AMI by their averaging method[VEB2010]. Their ‘sqrt’ and ‘sum’ averages are the geometric and arithmeticmeans; we use these more broadly common names.

References

[VEB2009]

Vinh, Epps, and Bailey, (2009). “Information theoretic measuresfor clusterings comparison”. Proceedings of the 26th Annual InternationalConference on Machine Learning - ICML ‘09.doi:10.1145/1553374.1553511. ISBN9781605585161.

[VEB2010]

Vinh, Epps, and Bailey, (2010). “Information Theoretic Measuresfor Clusterings Comparison: Variants, Properties, Normalization andCorrection for Chance”. JMLR<https://jmlr.csail.mit.edu/papers/volume11/vinh10a/vinh10a.pdf>

[YAT2016]

Yang, Algesheimer, and Tessone, (2016). “A comparative analysisof community detection algorithms on artificial networks”. ScientificReports 6: 30750.doi:10.1038/srep30750.

2.3.11.3.Homogeneity, completeness and V-measure#

Given the knowledge of the ground truth class assignments of the samples,it is possible to define some intuitive metric using conditional entropyanalysis.

In particular Rosenberg and Hirschberg (2007) define the following twodesirable objectives for any cluster assignment:

  • homogeneity: each cluster contains only members of a single class.

  • completeness: all members of a given class are assigned to the samecluster.

We can turn those concept as scoreshomogeneity_score andcompleteness_score. Both are bounded below by 0.0 and above by1.0 (higher is better):

>>>fromsklearnimportmetrics>>>labels_true=[0,0,0,1,1,1]>>>labels_pred=[0,0,1,1,2,2]>>>metrics.homogeneity_score(labels_true,labels_pred)0.66>>>metrics.completeness_score(labels_true,labels_pred)0.42

Their harmonic mean calledV-measure is computed byv_measure_score:

>>>metrics.v_measure_score(labels_true,labels_pred)0.516

This function’s formula is as follows:

\[v = \frac{(1 + \beta) \times \text{homogeneity} \times \text{completeness}}{(\beta \times \text{homogeneity} + \text{completeness})}\]

beta defaults to a value of 1.0, but for using a value less than 1 for beta:

>>>metrics.v_measure_score(labels_true,labels_pred,beta=0.6)0.547

more weight will be attributed to homogeneity, and using a value greater than 1:

>>>metrics.v_measure_score(labels_true,labels_pred,beta=1.8)0.48

more weight will be attributed to completeness.

The V-measure is actually equivalent to the mutual information (NMI)discussed above, with the aggregation function being the arithmetic mean[B2011].

Homogeneity, completeness and V-measure can be computed at once usinghomogeneity_completeness_v_measure as follows:

>>>metrics.homogeneity_completeness_v_measure(labels_true,labels_pred)(0.67, 0.42, 0.52)

The following clustering assignment is slightly better, since it ishomogeneous but not complete:

>>>labels_pred=[0,0,0,1,2,2]>>>metrics.homogeneity_completeness_v_measure(labels_true,labels_pred)(1.0, 0.68, 0.81)

Note

v_measure_score issymmetric: it can be used to evaluatetheagreement of two independent assignments on the same dataset.

This is not the case forcompleteness_score andhomogeneity_score: both are bound by the relationship:

homogeneity_score(a,b)==completeness_score(b,a)

Advantages:

  • Bounded scores: 0.0 is as bad as it can be, 1.0 is a perfect score.

  • Intuitive interpretation: clustering with bad V-measure can bequalitatively analyzed in terms of homogeneity and completeness tobetter feel what ‘kind’ of mistakes is done by the assignment.

  • No assumption is made on the cluster structure: can be used to compareclustering algorithms such as k-means which assumes isotropic blob shapeswith results of spectral clustering algorithms which can find cluster with“folded” shapes.

Drawbacks:

  • The previously introduced metrics arenot normalized with regards torandom labeling: this means that depending on the number of samples,clusters and ground truth classes, a completely random labeling will notalways yield the same values for homogeneity, completeness and hencev-measure. In particularrandom labeling won’t yield zero scoresespecially when the number of clusters is large.

    This problem can safely be ignored when the number of samples is more than athousand and the number of clusters is less than 10.For smaller samplesizes or larger number of clusters it is safer to use an adjusted index suchas the Adjusted Rand Index (ARI).

../_images/sphx_glr_plot_adjusted_for_chance_measures_001.png

  • These metricsrequire the knowledge of the ground truth classes whilealmost never available in practice or requires manual assignment by humanannotators (as in the supervised learning setting).

Examples

Mathematical formulation#

Homogeneity and completeness scores are formally given by:

\[h = 1 - \frac{H(C|K)}{H(C)}\]
\[c = 1 - \frac{H(K|C)}{H(K)}\]

where\(H(C|K)\) is theconditional entropy of the classes given thecluster assignments and is given by:

\[H(C|K) = - \sum_{c=1}^{|C|} \sum_{k=1}^{|K|} \frac{n_{c,k}}{n}\cdot \log\left(\frac{n_{c,k}}{n_k}\right)\]

and\(H(C)\) is theentropy of the classes and is given by:

\[H(C) = - \sum_{c=1}^{|C|} \frac{n_c}{n} \cdot \log\left(\frac{n_c}{n}\right)\]

with\(n\) the total number of samples,\(n_c\) and\(n_k\) thenumber of samples respectively belonging to class\(c\) and cluster\(k\), and finally\(n_{c,k}\) the number of samples from class\(c\) assigned to cluster\(k\).

Theconditional entropy of clusters given class\(H(K|C)\) and theentropy of clusters\(H(K)\) are defined in a symmetric manner.

Rosenberg and Hirschberg further defineV-measure as theharmonic mean ofhomogeneity and completeness:

\[v = 2 \cdot \frac{h \cdot c}{h + c}\]

References

2.3.11.4.Fowlkes-Mallows scores#

The original Fowlkes-Mallows index (FMI) was intended to measure the similaritybetween two clustering results, which is inherently an unsupervised comparison.The supervised adaptation of the Fowlkes-Mallows index(as implemented insklearn.metrics.fowlkes_mallows_score) can be usedwhen the ground truth class assignments of the samples are known.The FMI is defined as the geometric mean of the pairwise precision and recall:

\[\text{FMI} = \frac{\text{TP}}{\sqrt{(\text{TP} + \text{FP}) (\text{TP} + \text{FN})}}\]

In the above formula:

  • TP (True Positive): The number of pairs of points that are clustered togetherboth in the true labels and in the predicted labels.

  • FP (False Positive): The number of pairs of points that are clustered togetherin the predicted labels but not in the true labels.

  • FN (False Negative): The number of pairs of points that are clustered togetherin the true labels but not in the predicted labels.

The score ranges from 0 to 1. A high value indicates a good similaritybetween two clusters.

>>>fromsklearnimportmetrics>>>labels_true=[0,0,0,1,1,1]>>>labels_pred=[0,0,1,1,2,2]
>>>metrics.fowlkes_mallows_score(labels_true,labels_pred)0.47140

One can permute 0 and 1 in the predicted labels, rename 2 to 3 and getthe same score:

>>>labels_pred=[1,1,0,0,3,3]>>>metrics.fowlkes_mallows_score(labels_true,labels_pred)0.47140

Perfect labeling is scored 1.0:

>>>labels_pred=labels_true[:]>>>metrics.fowlkes_mallows_score(labels_true,labels_pred)1.0

Bad (e.g. independent labelings) have zero scores:

>>>labels_true=[0,1,2,0,3,4,5,1]>>>labels_pred=[1,1,0,0,2,2,2,2]>>>metrics.fowlkes_mallows_score(labels_true,labels_pred)0.0

Advantages:

  • Random (uniform) label assignments have a FMI score close to 0.0 for anyvalue ofn_clusters andn_samples (which is not the case for rawMutual Information or the V-measure for instance).

  • Upper-bounded at 1: Values close to zero indicate two label assignmentsthat are largely independent, while values close to one indicate significantagreement. Further, values of exactly 0 indicatepurely independentlabel assignments and a FMI of exactly 1 indicates that the two labelassignments are equal (with or without permutation).

  • No assumption is made on the cluster structure: can be used to compareclustering algorithms such as k-means which assumes isotropic blob shapeswith results of spectral clustering algorithms which can find cluster with“folded” shapes.

Drawbacks:

  • Contrary to inertia,FMI-based measures require the knowledge of theground truth classes while almost never available in practice or requiresmanual assignment by human annotators (as in the supervised learningsetting).

References#

2.3.11.5.Silhouette Coefficient#

If the ground truth labels are not known, evaluation must be performed usingthe model itself. The Silhouette Coefficient(sklearn.metrics.silhouette_score)is an example of such an evaluation, where ahigher Silhouette Coefficient score relates to a model with better definedclusters. The Silhouette Coefficient is defined for each sample and is composedof two scores:

  • a: The mean distance between a sample and all other points in the sameclass.

  • b: The mean distance between a sample and all other points in thenextnearest cluster.

The Silhouette Coefficients for a single sample is then given as:

\[s = \frac{b - a}{max(a, b)}\]

The Silhouette Coefficient for a set of samples is given as the mean of theSilhouette Coefficient for each sample.

>>>fromsklearnimportmetrics>>>fromsklearn.metricsimportpairwise_distances>>>fromsklearnimportdatasets>>>X,y=datasets.load_iris(return_X_y=True)

In normal usage, the Silhouette Coefficient is applied to the results of acluster analysis.

>>>importnumpyasnp>>>fromsklearn.clusterimportKMeans>>>kmeans_model=KMeans(n_clusters=3,random_state=1).fit(X)>>>labels=kmeans_model.labels_>>>metrics.silhouette_score(X,labels,metric='euclidean')0.55

Advantages:

  • The score is bounded between -1 for incorrect clustering and +1 for highlydense clustering. Scores around zero indicate overlapping clusters.

  • The score is higher when clusters are dense and well separated, whichrelates to a standard concept of a cluster.

Drawbacks:

  • The Silhouette Coefficient is generally higher for convex clusters thanother concepts of clusters, such as density based clusters like thoseobtained through DBSCAN.

Examples

References#

2.3.11.6.Calinski-Harabasz Index#

If the ground truth labels are not known, the Calinski-Harabasz index(sklearn.metrics.calinski_harabasz_score) - also known as the VarianceRatio Criterion - can be used to evaluate the model, where a higherCalinski-Harabasz score relates to a model with better defined clusters.

The index is the ratio of the sum of between-clusters dispersion and ofwithin-cluster dispersion for all clusters (where dispersion is defined as thesum of distances squared):

>>>fromsklearnimportmetrics>>>fromsklearn.metricsimportpairwise_distances>>>fromsklearnimportdatasets>>>X,y=datasets.load_iris(return_X_y=True)

In normal usage, the Calinski-Harabasz index is applied to the results of acluster analysis:

>>>importnumpyasnp>>>fromsklearn.clusterimportKMeans>>>kmeans_model=KMeans(n_clusters=3,random_state=1).fit(X)>>>labels=kmeans_model.labels_>>>metrics.calinski_harabasz_score(X,labels)561.59

Advantages:

  • The score is higher when clusters are dense and well separated, whichrelates to a standard concept of a cluster.

  • The score is fast to compute.

Drawbacks:

  • The Calinski-Harabasz index is generally higher for convex clusters thanother concepts of clusters, such as density based clusters like thoseobtained through DBSCAN.

Mathematical formulation#

For a set of data\(E\) of size\(n_E\) which has been clustered into\(k\) clusters, the Calinski-Harabasz score\(s\) is defined as theratio of the between-clusters dispersion mean and the within-clusterdispersion:

\[s = \frac{\mathrm{tr}(B_k)}{\mathrm{tr}(W_k)} \times \frac{n_E - k}{k - 1}\]

where\(\mathrm{tr}(B_k)\) is trace of the between group dispersion matrixand\(\mathrm{tr}(W_k)\) is the trace of the within-cluster dispersionmatrix defined by:

\[W_k = \sum_{q=1}^k \sum_{x \in C_q} (x - c_q) (x - c_q)^T\]
\[B_k = \sum_{q=1}^k n_q (c_q - c_E) (c_q - c_E)^T\]

with\(C_q\) the set of points in cluster\(q\),\(c_q\) thecenter of cluster\(q\),\(c_E\) the center of\(E\), and\(n_q\) the number of points in cluster\(q\).

References#

2.3.11.7.Davies-Bouldin Index#

If the ground truth labels are not known, the Davies-Bouldin index(sklearn.metrics.davies_bouldin_score) can be used to evaluate themodel, where a lower Davies-Bouldin index relates to a model with betterseparation between the clusters.

This index signifies the average ‘similarity’ between clusters, where thesimilarity is a measure that compares the distance between clusters with thesize of the clusters themselves.

Zero is the lowest possible score. Values closer to zero indicate a betterpartition.

In normal usage, the Davies-Bouldin index is applied to the results of acluster analysis as follows:

>>>fromsklearnimportdatasets>>>iris=datasets.load_iris()>>>X=iris.data>>>fromsklearn.clusterimportKMeans>>>fromsklearn.metricsimportdavies_bouldin_score>>>kmeans=KMeans(n_clusters=3,random_state=1).fit(X)>>>labels=kmeans.labels_>>>davies_bouldin_score(X,labels)0.666

Advantages:

  • The computation of Davies-Bouldin is simpler than that of Silhouette scores.

  • The index is solely based on quantities and features inherent to the datasetas its computation only uses point-wise distances.

Drawbacks:

  • The Davies-Bouldin index is generally higher for convex clusters than otherconcepts of clusters, such as density-based clusters like thoseobtained from DBSCAN.

  • The usage of centroid distance limits the distance metric to Euclideanspace.

Mathematical formulation#

The index is defined as the average similarity between each cluster\(C_i\)for\(i=1, ..., k\) and its most similar one\(C_j\). In the context ofthis index, similarity is defined as a measure\(R_{ij}\) that trades off:

  • \(s_i\), the average distance between each point of cluster\(i\) andthe centroid of that cluster – also known as cluster diameter.

  • \(d_{ij}\), the distance between cluster centroids\(i\) and\(j\).

A simple choice to construct\(R_{ij}\) so that it is nonnegative andsymmetric is:

\[R_{ij} = \frac{s_i + s_j}{d_{ij}}\]

Then the Davies-Bouldin index is defined as:

\[DB = \frac{1}{k} \sum_{i=1}^k \max_{i \neq j} R_{ij}\]
References#

2.3.11.8.Contingency Matrix#

Contingency matrix (sklearn.metrics.cluster.contingency_matrix)reports the intersection cardinality for every true/predicted cluster pair.The contingency matrix provides sufficient statistics for all clusteringmetrics where the samples are independent and identically distributed andone doesn’t need to account for some instances not being clustered.

Here is an example:

>>>fromsklearn.metrics.clusterimportcontingency_matrix>>>x=["a","a","a","b","b","b"]>>>y=[0,0,1,1,2,2]>>>contingency_matrix(x,y)array([[2, 1, 0],       [0, 1, 2]])

The first row of the output array indicates that there are three samples whosetrue cluster is “a”. Of them, two are in predicted cluster 0, one is in 1,and none is in 2. And the second row indicates that there are three sampleswhose true cluster is “b”. Of them, none is in predicted cluster 0, one is in1 and two are in 2.

Aconfusion matrix for classification is a squarecontingency matrix where the order of rows and columns correspond to a listof classes.

Advantages:

  • Allows to examine the spread of each true cluster across predicted clustersand vice versa.

  • The contingency table calculated is typically utilized in the calculation ofa similarity statistic (like the others listed in this document) between thetwo clusterings.

Drawbacks:

  • Contingency matrix is easy to interpret for a small number of clusters, butbecomes very hard to interpret for a large number of clusters.

  • It doesn’t give a single metric to use as an objective for clusteringoptimisation.

References#

2.3.11.9.Pair Confusion Matrix#

The pair confusion matrix(sklearn.metrics.cluster.pair_confusion_matrix) is a 2x2similarity matrix

\[\begin{split}C = \left[\begin{matrix}C_{00} & C_{01} \\C_{10} & C_{11}\end{matrix}\right]\end{split}\]

between two clusterings computed by considering all pairs of samples andcounting pairs that are assigned into the same or into different clustersunder the true and predicted clusterings.

It has the following entries:

\(C_{00}\) : number of pairs with both clusterings having the samplesnot clustered together

\(C_{10}\) : number of pairs with the true label clustering having thesamples clustered together but the other clustering not having the samplesclustered together

\(C_{01}\) : number of pairs with the true label clustering not havingthe samples clustered together but the other clustering having the samplesclustered together

\(C_{11}\) : number of pairs with both clusterings having the samplesclustered together

Considering a pair of samples that is clustered together a positive pair,then as in binary classification the count of true negatives is\(C_{00}\), false negatives is\(C_{10}\), true positives is\(C_{11}\) and false positives is\(C_{01}\).

Perfectly matching labelings have all non-zero entries on thediagonal regardless of actual label values:

>>>fromsklearn.metrics.clusterimportpair_confusion_matrix>>>pair_confusion_matrix([0,0,1,1],[0,0,1,1])array([[8, 0],       [0, 4]])
>>>pair_confusion_matrix([0,0,1,1],[1,1,0,0])array([[8, 0],       [0, 4]])

Labelings that assign all classes members to the same clustersare complete but may not always be pure, hence penalized, andhave some off-diagonal non-zero entries:

>>>pair_confusion_matrix([0,0,1,2],[0,0,1,1])array([[8, 2],       [0, 2]])

The matrix is not symmetric:

>>>pair_confusion_matrix([0,0,1,1],[0,0,1,2])array([[8, 0],       [2, 2]])

If classes members are completely split across different clusters, theassignment is totally incomplete, hence the matrix has all zerodiagonal entries:

>>>pair_confusion_matrix([0,0,0,0],[0,1,2,3])array([[ 0,  0],       [12,  0]])
References#
On this page

This Page