PaperlayerView exact revision on arXiv

Source-linked AI summary

Fast and Eager k-Medoids Clustering: O(k) Runtime Improvement of the PAM, CLARA, and CLARANS Algorithms

Erich Schubert, Peter J. Rousseeuw

arXiv:2008.05171v2cs.LGcs.AIstat.ML

TL;DR

PAM supports clustering with arbitrary dissimilarities but has high runtime costs. This paper accelerates its SWAP phase through caching and eager swapping, with theoretical and experimental evidence of substantial speedups and comparable quality, while extending the approach to CLARA and CLARANS.

  • Problem

    PAM uses medoids for arbitrary dissimilarities, but its SWAP phase has high runtime costs and k-medoids methods can become trapped in local optima.

  • Method

    The paper modifies PAM by caching partial results to eliminate a nested loop of length k, and studies eager swapping plus integrations with CLARA and CLARANS.

  • Results

    The modified PAM achieves a provable O(k)-fold SWAP speedup, while eager swapping substantially reduces iterations without observed quality loss; related PAM-based methods also benefit.

  • Takeaways & Limitations

    The speedups enable PAM to handle much larger problems, particularly those with large k, and extend to methods such as CLARA and CLARANS.

  • Takeaways & Limitations

    Performing only a few swaps reduces runtime but can produce significantly worse clustering results, reflecting a quality–runtime trade-off in simplified alternatives.

Abstract

from arXiv · show

Clustering non-Euclidean data is difficult, and one of the most used algorithms besides hierarchical clustering is the popular algorithm Partitioning Around Medoids (PAM), also simply referred to as k-medoids clustering. In Euclidean geometry the mean-as used in k-means-is a good estimator for the cluster center, but this does not exist for arbitrary dissimilarities. PAM uses the medoid instead, the object with the smallest dissimilarity to all others in the cluster. This notion of centrality can be used with any (dis-)similarity, and thus is of high relevance to many domains and applications. A key issue with PAM is its high run time cost. We propose modifications to the PAM algorithm that achieve an O(k)-fold speedup in the second ("SWAP") phase of the algorithm, but will still find the same results as the original PAM algorithm. If we relax the choice of swaps performed (while retaining comparable quality), we can further accelerate the algorithm by eagerly performing additional swaps in each iteration. With the substantially faster SWAP, we can now explore faster initialization strategies, because (i) the classic ("BUILD") initialization now becomes the bottleneck, and (ii) our swap is fast enough to compensate for worse starting conditions. We also show how the CLARA and CLARANS algorithms benefit from the proposed modifications. While we do not study the parallelization of our approach in this work, it can easily be combined with earlier approaches to use PAM and CLARA on big data (some of which use PAM as a subroutine, hence can immediately benefit from these improvements), where the performance with high k becomes increasingly important. In experiments on real data with k=100,200, we observed a 458x respectively 1191x speedup compared to the original PAM SWAP algorithm, making PAM applicable to larger data sets, and in particular to higher k.

1. Introduction

The paper motivates k-medoids as a flexible alternative for clustering arbitrary dissimilarities, while introducing PAM's computational cost and this work's speedup-focused extensions.

  • 1. Introduction: K-medoids uses representative data objects rather than means, allowing arbitrary dissimilarities and input domains beyond vector spaces.The medoid is the cluster object minimizing total dissimilarity to other objects.
  • 1. Introduction: PAM is a widely known medoid-based algorithm whose BUILD initialization and SWAP refinement address the NP-hard clustering problem through local search.PAM requires a dissimilarity matrix, which can itself be costly to compute and store.
  • 1. Introduction: The improved algorithm removes SWAP's nested loop of length k, enabling a provable O(k)-fold speedup over original PAM.The paper also studies eager swapping, updated benchmarks, and additional related work.
  • 1. Introduction: The work extends earlier conference research with a provable speedup, eager swapping, simpler initialization, and broader evaluation.It is an extended version of the SISAP’19 conference paper.
  • 1. Introduction: PAM originated for arbitrary dissimilarity matrices, including data that need not satisfy the triangle inequality.Related operations-research formulations include p-medians and facility-location problems.

2. Partitioning Around Medoids (PAM) and its Variants

PAM builds an initial medoid set and repeatedly improves it through candidate swaps, but its SWAP phase and some alternative heuristics have important computational or quality trade-offs.

  • 2. Partitioning Around Medoids (PAM) and its Variants: PAM combines BUILD, which selects an initial clustering, with SWAP, which iteratively improves it toward a local optimum.Finding the global optimum of k-medoids is NP-hard.
  • 2. Partitioning Around Medoids (PAM) and its Variants: BUILD greedily selects k medoids by first choosing the smallest distance sum, then repeatedly adding the point that reduces total deviation most.The procedure tracks the best change in total deviation at each selection step.
  • 2. Partitioning Around Medoids (PAM) and its Variants: SWAP considers k · (n − k) candidate replacements and applies the best improving exchange until no further improvement exists.Because the medoid search space is finite, this steepest-descent process converges after finitely many iterations.
  • 2. Partitioning Around Medoids (PAM) and its Variants: Caching nearest and second-nearest medoid distances reduces the main SWAP loop to O(k(n−k)^2) per iteration.Cached values must be updated after each performed swap.
  • 2. Partitioning Around Medoids (PAM) and its Variants: CLARA repeatedly runs PAM on a subsample and assigns remaining objects to their closest medoid, reducing runtime to about O(k^3 + n) when the sample size is O(k).The suggested sample size is n′ = 40 + 2k.
  • 2.2. Alternating k-medoids Algorithm: The alternating k-medoids heuristic can get stuck because new medoids must cover entire current clusters, missing improvements enabled by reassignment during SWAP.Allowing reassignment produces a restricted SWAP variant with little benefit over accelerated SWAP.

3. Finding the Best Swap

The paper accelerates PAM's SWAP phase by exploiting redundant computations, removing the nested medoid loop while preserving exact PAM results in FastPAM1. Eager swapping further accelerates local search, and the resulting faster SWAP enables alternative initialization and FastCLARA integration.

  • Best-swap computation: PAM evaluates every medoid–non-medoid swap, but recomputing total deviation repeatedly causes redundant nearest-medoid calculations.The original computation considers k · (n − k) candidate swaps and repeatedly evaluates assignments.
  • Best-swap computation: The swap-loss calculation uses each object's nearest and second-nearest medoid distances to evaluate candidate swaps in O(1) per object.With nearest(o), d_n(o), and d_s(o) known, ΔTD can be computed in O(n − k) rather than O(nk).
  • FastPAM optimization: FastPAM removes the nested medoid loop by sharing computations across medoids, eliminating the factor k from SWAP runtime.The common cases produce the same result for all non-nearest medoids, allowing a shared accumulator and cached removal losses.
  • FastPAM variants: FastPAM1 implements optimizations A–C and guarantees the same results as original PAM, while FasterPAM additionally performs eager swaps.The combined techniques include nested-loop removal, a shared accumulator, removal-loss precomputation, and eager execution.
  • Eager swapping: Eager variants apply improving swaps immediately and continue until a full data pass finds no improvement, potentially performing many swaps per iteration.Selecting the best of k choices can produce different results, with quality depending on the smoothness of the solution space.
  • Initialization and integration: Fast SWAP changes PAM's bottleneck to BUILD, whose O(kn^2) complexity dominates after SWAP falls from O(k(n − k)^2) to O((n − k)n).For k = 100, SWAP falls from 95% of PAM runtime to 3.7% with FasterPAM, while BUILD rises to 91%.
  • Initialization and integration: The faster SWAP permits worse starting conditions, and the paper recommends uniform random or distance-weighted initialization with FasterPAM.LAB reduced swaps for PAM and FastPAM1, but that benefit was offset by eager swapping.

4. Experiments

The experiments evaluate initialization quality and runtime, algorithm–initialization combinations, practical k-dependent speedup, eager-swapping behavior, approximation, replication, and scalability. The paper expects the O(k) speedup to be directly measurable despite implementation constants.

  • Experimental scope: All evaluated algorithms are local-search methods that may get stuck in local minima; considered changes primarily affect quality, while best-versus-first selection primarily affects performance.This frames the experiment's separation of solution quality from runtime behavior.
  • Research questions: The experiments ask how initialization affects result quality, runtime, and favorable combinations with the clustering algorithms.These are the first three stated research questions.
  • Research questions: The study measures practical speedup over original PAM as a function of k, whose theoretical speedup is O(k).The authors explicitly aim to test whether the theoretical dependence appears in practice.
  • Research questions: The experiments compare eager swapping with optimal swapping by measuring saved iterations and additional swaps until convergence.This is the fifth stated research question.
  • Practical speedup: The proposed O(k)-fold speedup is expected to be immediately measurable because the fixed-array-cache overhead is small.The paper notes that constant factors can still determine the minimum k at which the method becomes beneficial.
  • Experimental scope: The evaluation covers initialization, practical speedup, eager swapping, approximative algorithms, replication on another data set, and scalability with data-set size.Different research questions are assigned to Sections 4.3–4.8.

4.2. Data Sets

The evaluation combines small k-median instances with known optima, larger UCI data sets for scalability, and an ELKI-based Java implementation tested on Xeon servers. The data range includes plant leaves, handwritten digits, and MNIST.

  • Benchmark data: The study uses 40 classic OR-Library k-median problems whose true optima are known, allowing solution gaps to be computed.These instances contain up to 900 cases by the paper's description.
  • Scalability data: Scalability experiments use UCI data sets because the classic optimization problems are too small for that purpose.The selected collections include one-hundred plant species leaves, optical handwritten digits, and MNIST.
  • Scalability data: The plant-leaves data set has 1,600 instances and 100 classes, while the optical-digits data set has 5,620 instances, 64 variables, and 10 natural classes.The paper repeats the leaf experiment on optical digits and MNIST.
  • Scalability data: MNIST has 60,000 instances and 784 pixel variables, a size that regular PAM cannot handle in reasonable time.The 784 variables correspond to pixels in a 28×28 image grid.
  • Implementation: The implementation was developed in Java using the ELKI open-source data-mining toolkit and executed on Intel Xeon E5-2697v2 CPUs at 2.70 GHz.Earlier R-cluster-package replications are omitted as redundant.
  • Implementation: The implementation includes possible additional initialization integration and special-case optimizations for k = 2, but these are not presented as the main evaluation setting.Some closest and second-medoid distances can be reused during initialization; k = 2 permits further optimizations.

4.3. Initialization

Initialization quality and cost vary substantially, but PAM's swapping phase can recover good solutions from poor starts while runtime remains sensitive to initialization. FasterPAM makes cheaper randomized starts more attractive by performing multiple swaps per iteration.

  • Initialization quality: GreedyG and BUILD provide the strongest standalone initialization quality, finding optimum solutions for 6 and 4 problems, respectively.Both are deterministic and costly; GreedyG adds a refinement step beyond direct TD optimization.
  • PAM after initialization: PAM largely removes initialization-quality differences in the final clustering, because swaps can replace poor initial medoids with better alternatives.Initialization primarily affects runtime after PAM optimization.
  • PAM after initialization: GreedyG achieves the fastest PAM runtime despite being the slowest initialization, because it substantially reduces the swaps needed for convergence.Park's initialization has good average quality but the worst runtime, while far-point and overly central heuristics can be worse than random.
  • Random restarts and Alternating: Ten random restarts improve randomized initializations, while the Alternating algorithm remains highly dependent on good starts and can stay close to random quality.The supplied comparison reports Alternating at 69.7% with Park and Jun initialization.
  • FasterPAM initialization: With FasterPAM, LAB's lower swap count is offset by eager swapping, so uniform random or distance-weighted initialization is recommended.Multiple swaps per iteration reduce the impact of poor starting conditions and make random initializations more attractive.

4.4. Run Time Speedup with Increasing k

As k increases, FastPAM and FasterPAM reduce PAM SWAP runtime substantially, with FasterPAM scaling especially well. The improvements also make random initialization competitive when complete runtime includes distance computation and initialization.

  • Scaling with k: FasterPAM's runtime increases only little with k, making it particularly attractive for larger cluster counts.The theoretical speedup of FastPAM and FasterPAM is O(k), and the experiment varies k from 2 to 200.
  • SWAP runtime: For k = 200, FasterPAM's SWAP takes 169 milliseconds, compared with 2.4 seconds for FastPAM1 and 151 seconds for PAM.These measurements use random initialization for the faster methods; PAM uses BUILD, whose time is excluded from the PAM figure.
  • Speedup: FastPAM1 achieves about 0.75·k empirical speedup on the tested data set, while iteration reductions add roughly 2-7× for the eager variants.The most extreme measured SWAP speedup is about 1190× at k = 200 with BUILD initialization.
  • Complete runtime: Complete-runtime comparisons retain the fundamental speedup pattern, although FasterPAM spends substantial time computing the distance matrix.Including distance-matrix computation and initialization makes the overall FasterPAM runtime appear almost constant in the plotted comparison.

4.5. Number of Iterations

PAM generally converges in few iterations, but iteration counts depend on k and initialization. Eager variants reduce iterations by applying multiple swaps per pass, making additional iterations from cheaper starts less costly.

  • PAM iterations: PAM usually needs only a few iterations, although its single-swap-per-iteration design implies iteration growth with k.The number of iterations required by PAM lacks established theoretical results, and worst-case behavior may be superpolynomial.
  • PAM iterations: BUILD initialization often requires fewer than k PAM iterations, leaving many medoids unchanged from their initial values on the smaller ORLib data sets.PAM's iteration count equals its number of swaps plus one final iteration without an improving swap.
  • Eager swapping: EagerPAM and FasterPAM drastically reduce iteration counts by performing multiple swaps per pass over the data set.On the examined data set, EagerPAM's maximum was 13 iterations and its worst average was 9.9.
  • Alternating: The Alternating approach is omitted from the iteration figures because it uses few iterations but produces substantially worse clustering results.The paper compares it separately with subsample-based algorithms such as CLARA.

4.6. Quality

The PAM variants generally find similarly high-quality solutions, while initialization affects individual outcomes and runtime. Eager variants may converge to different local optima, but are expected to have comparable quality.

  • EagerPAM and FasterPAM can converge to different solutions because they perform the first improving swap rather than the optimal swap.The implementations pass the same unit tests as PAM and FastPAM1 is expected to produce the exact same result as PAM.
  • PAM and its variants, including random initialization, find results of similar quality.Deterministic BUILD or GreedyG initialization is usually better than the worst random-initialization solution.
  • Initialization quality varies with k: BUILD works well for some values, while GreedyG or random initialization performs better for others.The experiment suggests considering BUILD and GreedyG for small k < 10, while also trying different starting conditions.

4.7. Optical Digits Dataset

On the Optical Recognition of Handwritten Digits data set, eager methods achieve better quality, while FasterPAM substantially reduces runtime relative to PAM.

  • The data set contains n = 5620 instances, d = 64 variables, and 10 natural classes.
  • FasterPAM is over 10 times faster at k = 10 and over 200 times faster total at k = 100.The latter comparison includes the time needed to compute the distance matrix.
  • Eager methods find better-quality solutions independently of initialization on the Optical Digits data set.The authors relate this to different local optima caused by processing order.
  • There is little benefit from using k-means++ or LAB instead of random initialization in this experiment.

4.8. Scalability Experiments

The scalability experiments show that FasterPAM remains quadratic in n because of the distance matrix, but greatly reduces clustering time and enables larger k when that matrix is feasible.

  • The main scalability bottleneck is the distance matrix’s O(n2) computation and memory, rather than clustering itself.The method’s quadratic runtime and memory requirements make it unsuitable for big data as-is.
  • FasterPAM is only slightly more expensive than computing the distance matrix after runtime normalization by n2.
  • For k = 100 and n = 35000, FasterPAM averages 743 seconds versus 21 hours for PAM, with an 879× speedup excluding matrix computation.About 655 seconds, or 88%, of FasterPAM’s runtime is spent computing the distance matrix.
  • FasterCLARA can scale linearly in n during final assignment when computing the distance matrix is prohibitive, but usually produces worse results.With k = 10, CLARA and Alternating are about 30% worse than the benchmark’s best results; with k = 100, CLARA is barely better than random medoids.
  • FastCLARANS achieves better quality than CLARANS at similar runtime by considering a k times larger search space.The authors recommend FastCLARANS mainly for inexpensive distance functions such as low-dimensional Euclidean distance.

5. Outlook

The paper connects k-medoids with facility location and identifies opportunities to transfer methods between the areas. Future work targets sparse instances with restricted consumer–supplier connections.

  • The authors identify a close relationship between classic k-medoids clustering and facility location problems.They note that the implementation is already prepared for the bichromatic case with separate consumers and possible locations.
  • Facility-location research may help cluster data without requiring the parameter k beforehand.Capacity constraints and facility opening costs are described as less obvious to integrate into k-medoids.
  • Future work will optimize sparse instances where not every consumer can be serviced from every possible supplier location.The paper gives power-network planning as an example in which connections are restricted by the road network and physical limitations.

6. Conclusions

The proposed FasterPAM modifications provide provable O(k)-fold acceleration through caching and simpler eager swapping, with evidence of reduced iterations without quality loss. These improvements also benefit related PAM-based methods and support larger data sets with higher k.

  • A provable O(k)-fold PAM speedup eliminates recomputation through cached partial results.The modification avoids recomputing partial results during swapping.
  • Eagerly executing the first improvement found substantially reduces the number of iterations without observed quality loss.The paper combines theoretical arguments with experimental evidence for this result.
  • Compared with earlier work, the speedup is now provable, while eager execution is simpler and more effective than FastPAM2.The authors also report that initialization is simpler in the improved version.
  • Uniform random initialization is often fastest with FasterPAM and attains high quality, whereas LAB and distance-weighted initialization show no systematic improvement.Distance-weighted initialization nevertheless yields similar performance and remains useful.
  • CLARA, CLARANS, and parallel or distributed variants benefit because they use PAM or employ similar swapping methods.The corresponding modifications can be applied to these methods.
  • Implementations in ELKI, Rust, Python, and R are intended to make the improvements usable on larger data sets with higher k.The availability of the implementations is presented as supporting broader use of the methods.

Appendix A. Proof of Restructured Equation

The appendix proves that the restructured swap-loss equation is equivalent to the original formulation by decomposing both into per-object contributions. A four-case analysis shows that each contribution matches, so summing over all objects establishes equality.

  • The proof decomposes both equations into individual loss contributions ∆(x_o, m_i, x_c) for each object x_o.This decomposition is used to prove equivalence for all x_o.
  • Four cases based on distance comparisons and nearest-medoid identity reduce the restructured expression to the corresponding original contribution.The cases cover whether the candidate is closer than the nearest or second-nearest medoid and whether it is the current nearest medoid.
  • Because every per-object contribution matches, summing over all objects proves equality between the two equations.The appendix concludes the argument with the resulting equality.
Loading 2008.05171v2…