Clustering Algorithms: A Comparative Approach

In plain words

Clustering methods sort unlabeled data into groups. There are many of them, and they each come with settings that someone has to pick. This paper asks which of seven common clustering methods, all from the R language, gives the best groups. It answers with a large set of made-up datasets whose true groups are known, so every method can be scored against the answer.

The authors test three situations. First, each method runs with its out-of-the-box settings. Second, one setting is changed at a time. Third, all settings are changed together at random. The defaults are not always good: a spectral method does best with defaults, but a few random tweaks help several other methods a lot. The paper is a practical guide for someone who uses clustering but is not a specialist.

TL;DR

• Compares seven R clustering methods on 270 synthetic datasets, scored by four agreement indices (Jaccard, ARI, Fowlkes-Mallows, NMI) §4.2 tab. 1.
• With default settings, the spectral method has the highest average ARI (65.80%) and hierarchical clustering the lowest (24.81%) tab. 2.
• Random search over all parameters helps EM and hierarchical clustering most on the two-class, two-feature set tab. 5. EM and hierarchical also gain strongly on the 10-feature sets tab. 7 tab. 8; most methods gain little on DB10C2F tab. 6.
• It is a benchmark, not a method to implement. Its main use is the data generator and the evaluation protocol, which are underspecified (see Pitfalls).

Problem

Clustering results depend on the method, on its parameters and on the data. Users often run software defaults without tuning, and parameter search is expensive and can overfit §1. The paper asks how the default configurations of common R clustering implementations perform, how sensitive each method is to single parameters, and whether random parameter sampling improves results §1. Earlier comparisons use real data, artificial data or both §2. This paper uses only artificial data with controlled properties, so many scenarios can be run §1 §4.1.

Method

The study has three parts: a dataset generator §4.1, a set of external agreement indices §4.2, and three evaluation protocols §5.

Algorithms. Seven methods, each with one R function tab. 1: k-means (kmeans, stats), clara (clara, cluster), hierarchical (agnes, cluster), EM (mstep/estep, mclust), hcmodel (hc, mclust), spectral (specc, kernlab) and subspace (hddc, HDclassif) tab. 1 §3. The paper describes each method in words §3; no pseudocode is given for any of them, so the R functions are the reference.

Dataset generator §4.1. Based on Hirschberger et al. §4.1:

1. Choose the number of classes C, features F, objects per class Ne and the separation parameter α §4.1.
2. For each class, build a random F×F matrix R = G Gᵀ, which is guaranteed to be positive semi-definite §4.1. The paper says the entries follow a normal distribution after a moment-matching step from Hirschberger et al., but the moments are not given §4.1.
3. Draw Ne points from the class, divide by α, and translate by a shift s_i drawn uniformly from [−1, 1] §4.1.
4. α is tuned per dataset so that no algorithm reaches 0% or 100% accuracy §4.1.

The pseudocode below is our reading of steps 1–3. Parts the paper does not state are marked.

# Artificial dataset generator  [§4.1]
# Input:  C classes, F features, Ne objects per class, separation alpha
# Output: points (C * Ne, F), labels (C * Ne)

def generate_dataset(C, F, Ne, alpha):
    # Build a valid covariance per class: R = G G^T  [§4.1]
    # Entry distribution and moments of G: not stated  [§4.1]
    for each class i in range(C):
        G = random_matrix(F, F)
        covariance_i = G @ G.T                   # (F, F)
        # Sample Ne objects from the class  [§4.1]
        # Mean of the Gaussian: not stated
        points_i = sample_gaussian(covariance_i, Ne)
        # Scale by 1/alpha and shift by s_i ~ U[-1, 1]  [§4.1]
        # Scalar or per-feature shift: not stated
        shift_i = uniform(-1, 1)
        points_i = points_i / alpha + shift_i
    return stack(points), labels

Agreement indices §4.2. Build a and b and c from pair counts between the reference partition U and the found partition V eq. 1–3:

• a: pairs together in both U and V eq. 1.
• b: pairs together in U but split in V eq. 2.
• c: pairs split in U but together in V eq. 3.

Then the four indices:

• Jaccard J = a / (a + b + c) eq. 4.
• Adjusted Rand index (ARI) eq. 5.
• Fowlkes–Mallows FM = a / sqrt((a+b)(a+c)), written in the paper as a·√(a+b)·√(a+c) / ((a+b)(a+c)) eq. 6.
• Normalized mutual information NMI = I(C,T) / sqrt([H(C), H(T)]) eq. 7. See Pitfalls for the notation.

Protocols §5:

• Default parameters: each method runs with its R default settings §5.1. Default K is set to the true number of classes [supp. S1].
• One-dimensional: vary one parameter at a time, others at default §5.2. Summary statistics: mean gain ⟨S⟩ eq. 8, max gain max S eq. 9 and standard deviation ΔS eq. 10.
• Multi-dimensional: for each method draw 500 parameter sets, each parameter uniform in bounds [P_min, P_max] taken from the one-dimensional results §5.3. Summaries are the share of draws that beat the default (p-value), the mean ⟨R⟩, standard deviation ΔR and maximum max R of the gains among improved draws §5.3.

# Random parameter search for one method  [§5.3]
# Input:  bounds (P_min, P_max) per parameter, default score, dataset
# Output: p-value (%), mean/std/max gain over improved draws

def random_search(method, dataset, bounds, default_ari):
    # Draw 500 parameter sets, uniform inside the bounds  [§5.3]
    gains = []
    for draw in range(500):
        params = sample_uniform(bounds)
        # Gain = ARI with the draw minus ARI with defaults  [§5.3]
        ari = score(method(dataset, params), dataset.labels)
        gains.append(ari - default_ari)
    # p-value = share of draws that improve on defaults  [§5.3]
    improved = [g for g in gains if g > 0]
    p_value = 100 * len(improved) / len(gains)
    # <R>, Delta R and max R over the improved draws  [§5.3]
    return p_value, mean(improved), std(improved), max(improved)

Hyperparameters

Dataset settings §4.1 §5.2. Default settings of each algorithm [supp. S1]. The experiments set K to the true number of classes when a method has no default for it [supp. S1].

Results

Author-reported. Setup: 270 synthetic datasets §4.1, seven R implementations tab. 1, four agreement indices §4.2, R default parameters unless stated. The paper does not give the R version.

Average ARI per method, over all datasets, with default parameters tab. 2, ARI block, MAcc column:

Summary of best ARI values, DB2C2F and DB10C2F tab. 9:

Summary of best ARI values, DB2C10F and DB10C10F tab. 10:

Headline claims, as stated by the authors:

• With default parameters, spectral averages best across the four indices; the paper reports a margin "of at least 10%" §5.1. See Pitfalls for the check.
• Kruskal–Wallis: no significant difference between methods for F = 2 (p = 0.07, χ² = 10.26); significant for F = 10 (p = 4.4×10⁻⁶) and F = 50 (p = 1.4×10⁻⁶) §5.1.
• One-dimensional: EM modelName gives ⟨S⟩ = 10.75% on DB2C2F; spectral kpar gives 7.36% §5.2 tab. 3. On DB2C10F, hierarchical method gives 8.76% and EM modelName 18.80% tab. 4.
• Multi-dimensional, DB2C2F: EM p-value 68.1%, mean gain 22.1%; hierarchical mean gain 30.6% tab. 5. DB10C2F: EM p-value 76.5%, k-means 77.7% tab. 6. DB2C10F: EM mean gain 30.1%, hierarchical 25.9% tab. 7. DB10C10F: EM mean gain 17.1%, hierarchical 32.1% tab. 8.
• Changing K: underestimating the number of classes was worse than overestimating it, for all methods §5.2 §6. fig. 5 shows this.

Paper vs code

No official code. The paper gives no repository or project page, and the arXiv metadata we hold has no comments field. Searched: (1) the paper and metadata for links, none found; (2) GitHub repository search for "clustering algorithms comparative R" and for "clustering chcomin" (the corresponding author's account), with no repository by an author; (3) the first 100 repositories of that account, with no comparison repository visible. The account's listing was truncated at the end, so this search is not exhaustive.

Pitfalls

• NMI notation. The paper writes NMI = I(C,T) / sqrt([H(C), H(T)]) eq. 7. The bracket is ambiguous and the paper does not say whether it means the geometric mean of H(C) and H(T). Check the reference the paper cites for NMI §4.2 before implementing eq. 7.
• Gain statistics are mislabelled in the caption. Table captions list ⟨S⟩, max S and ΔS as "average, standard deviation and maximum difference" tab. 3 tab. 4, which reverses the order of max S and ΔS. The definitions in eq. 9 and eq. 10 are authoritative. The text does not say whether the standard deviation is the sample or population form.
• Subspace K is contradictory. The supplement says that where a method has no default K, K is set to the true number of classes [supp. S1]. For hddc, it says K is chosen by maximum BIC in the default interval (1, 10] [supp. S1]. The paper does not say which was used in the experiments.
• Defaults used the true K. Default-parameter results set K to the true number of classes [supp. S1], which an applied user would not know. Results for mis-specified K are reported separately in §5.2 fig. 5.
• Generator details are missing. The moments used to build G, the distribution of the Gaussian mean, the sample size of covariances and the rule for choosing α are all not stated §4.1. The shift s_i is written as one value per class; whether it is a scalar or a vector is not stated §4.1. Reproducing the 270 datasets needs the original code, which is not available.
• Dataset naming. "DB2F" etc. pools all classes and object counts, while a dataset name like DB2C10F fixes C and F and uses Ne = 50 §4.1. Keep the pooled and fixed cases apart when reproducing.
• "At least 10%" claim. The margin over the next method is at least 10 points for NMI (smallest margin 12.49, against k-means) tab. 2. For ARI, Jaccard and FM, the margin over subspace is 7.22, 3.45 and 3.56 points tab. 2. The paper does not state whether margins are in percentage points or relative terms.
• Undefined MAcc column. Table 2 has a last column "MAcc" that the text describes only as the average performance of each algorithm tab. 2 §5.1. Assume it is the mean over the 270 datasets; the text does not confirm this.
• Parameter sweep values and bounds are not listed. The swept values in §5.2 and the bounds [P_min, P_max] used for random sampling in §5.3 are not stated. Random-search results cannot be matched without them.
• Seeds and repetitions. The number of random seeds per parameter set, the R version and the seed handling are not stated. Several methods are random (EM z, k-means nstart, subspace init), so single-run results may vary.
• Text vs table, small mismatch. §5.2 says the subspace losses from min.individuals and models are about −20%. The tables give −26.9 and −22.44 on DB2C10F and DB10C10F tab. 4. Use the tables.
• Hierarchical defaults. The default method is "average" (UPGMA), which produced very unbalanced clusters and poor accuracy on DB10C10F §5.2. Its default ARI on DB2C10F is 0.0 tab. 10.
• Costs the abstract hides. The study ran 500 random parameter sets per method per dataset §5.3, and one-dimensional sweeps over several parameters tab. 3 tab. 4. Reproducing the full study is expensive.

Lineage

• Builds on: Hirschberger et al. dataset generator, with a normal-distribution choice for the covariance entries §4.1. The generator is not a card in this wiki.
• Related approaches named in the paper: a prior work that uses the same generation procedure §4.1; graph-clustering and other clustering comparisons §2. None has a card in this wiki.
• Compared against: the paper compares its seven methods with each other only.
• Superseded by: not stated.