parallelpam

Juan Domingo

Load package

library(parallelpam)

Purpose

The parallelpam package (Domingo (2023b)) is meant as a way to apply the PAM algorithm to quite (or very) big sets of data, such as the results of single-cell sequencing, but can be generally used for any type of data, as long as a distance/dissimilarity matrix can be calculated.

Differently to other packages, its main strength is its ability to perform clustering based on Partitioning Around Medoids (PAM) using a large number of items and doing it in parallel. Memory and speed limitations are reduced by extensive use of C++ programming which allows use of alternative data types (for instance, float vs. double to represent distance/dissimilarity), intermediate disk-storage to avoid memory copy operations whenever possible and use of threads.

Both phases of PAM (initialization with BUILD and optimization) have been parallelized so it you have a multi-core machine, many threads will be launched with a great acceleration of the calculations. This is done automatically, even you are allowed to choose the number of threads yourself for comparison purposes or to allow your machine to execute other things simultaneously.

Also, calculation of the matrix of distances/dissimilarities from the initial data and calculation of silhouette of the resulting clusterization are available, too, and are calculated in parallel in a multi-core machine.

The data are stored in memory reading them from binary files created with the package jmatrix (Domingo (2023a)). To be familiar with them please read the vignette jmatrixpp which is included with this package.

WARNING: you must NOT load jmatrix explicitly. Indeed, you do not need even to install it. All its functions have been included here, too, so calling to library(parallelpam) is enough.`

Workflow

Debug messages

First of all, the package can show quite informative (but sometimes verbose) messages in the console. To turn on/off such messages you can use

# Initially, state of debug is FALSE. Turn it on exclusively for the
# parallelpam part with
ParallelpamSetDebug(TRUE)
#> Debugging for PAM algorithm set to ON.
# There is another parameter, debjmat, to turn on messages about
# binary matrix creation/manipulation. By default is FALSE but turn it on
# if you like with
# ParallelpamSetDebug(TRUE,debjmat=TRUE)

Data load/storage

The first step is to load raw data coming from external sources like the main formats used in single cell experiments which should have been stored as a binary matrix file in jmatrix format. Since this is a separate package, and for purposes of illustration, we will create an artificial matrix for a small problem that fits in R memory with 5000 vectors with 500 dimensions each. Then we will calculate the dissimilarity matrix and finally we will apply the PAM algorithm to it.

# Create the matrix with row names V1 to V5000 and column names d1 to d500
nvec<-5000
ndim<-500
P<-matrix(runif(nvec*ndim),nrow=nvec)
rownames(P)<-paste0("V",1:nvec)
colnames(P)<-paste0("d",1:ndim)
# Write it to disk as a binary file in jmatrix format. Please,
# see vignette jmatrixpp.
JWriteBin(P,"datatest.bin",dtype="float",dmtype="full",
          comment="Synthetic problem data to test PAM")

For your real problem, the input format can be a .csv file. See function CsvToJMat in package scellpam (Domingo (2023c)).

To know details about the generated files do

JMatInfo("datatest.bin")
#> File:               datatest.bin
#> Matrix type:        FullMatrix
#> Number of elements: 2500000
#> Data type:          float
#> Endianness:         little endian (same as this machine)
#> Number of rows:     5000
#> Number of columns:  500
#> Metadata:           Stored names of rows and columns.
#> Metadata comment:  "Synthetic problem data to test PAM"

Calculating the distance/dissimilarity matrix

This is the most computationally intensive part of the process (particularly, for samples with a high number of points and/or high dimensionality) and therefore has been programmed in parallel, taking advantage of the multiple cores of the machine, if available. The funcion is called CalcAndWriteDissimilarityMatrix. Its input and output files (first and second parameters) are of course compulsory. Input file can be a sparse of full binary jmatrix (but obviously, not a symmetric matrix).

WARNING: notice that the vectors to calculate dissimilarities amongst them MUST be stored BY ROWS. This is due to efficiency reasons.

Output dissimilarity matrix will always be a binary symmetric (obviously square) matrix with a number of rows (and columns) equal to the number of rows (in this case, vectors) of the input file. The type of distance/dissimilarity can be L1 (Manhattan distance), L2 (Euclidean distance) or Pearson (Pearson dissimilarity coefficient). The resulting matrix stores only names for the rows, which are the names of the vectors stored as rows in file datatest.bin. If the number of vectors is \(N\), only \(N(N+1)/2\) dissimilarity values are really stored.

A note on the number of threads, valid also for other algorithms that will be explained later:

Possible values for the number of threads are:

Choosing explicitly a number of threads bigger than the number of available cores is allowed, but discouraged and the program emits a warning about it.

With respect to option 0, for small problems (in this case, less than 1000 vectors) the function makes the choice of not using threads, since the overhead of opening and waiting termination is not worth. For bigger problems the number of chosen threads is the number of available cores, or twice this number if the processor is capable of hyperthreading. Nevertheless, this choice may not be the best, depending on your machine, possibly due (I guess) to the memory access conflicts created by the need of keep processor cache coherence. You may have to do some trials with your data in your machine.

Now, let us try it with this small dataset.

CalcAndWriteDissimilarityMatrix("datatest.bin","datatestL2.bin",
                                distype="L2",restype="float",
                                comment="L2 distance for vectors in
 jmatrix file datatest.bin",nthreads=0)
#> Input matrix is a full matrix  with elements of type 'float' and size (5000,500)
#> Read full matrix from file datatest.bin. Its size is [5000 x 500] and it uses 2500000 elements of 4 bytes each with accounts for 9.53674 MBytes.
#> Creating dissimilarity matrix of size (5000x5000)
#> Loading required package: memuse
#> Package memuse is installed. OK.
#>   Memory used by the matrix: 48837 KiB, which is 0.01% of the available memory, which is 259181160 Kib.
#>   That seems OK.
#> Launching up to 128 simultaneous threads.
#> End of dissimilarity matrix calculation (parallel version). Elapsed time: 0.290003 s
#> Output binary file datatestL2.bin written.

The resulting matrix is stored as a binary symmetric matrix of float values, as we can check.

JMatInfo("datatestL2.bin")
#> File:               datatestL2.bin
#> Matrix type:        SymmetricMatrix
#> Number of elements: 25000000 (12502500 really stored)
#> Data type:          float
#> Endianness:         little endian (same as this machine)
#> Number of rows:     5000
#> Number of columns:  5000
#> Metadata:           Stored only names of rows.
#> Metadata comment:  "L2 distance for vectors in
#>  jmatrix file datatest.bin"

Applying PAM

The last step is to take the previously calculated matrix and apply the Partitioning Around Medoids classifier. Function is ApplyPAM. First parameter (name of the file containing the dissimilarity matrix in jmatrix format) and second parameter (k, number of medoids) are compulsory. The names and default values for the rest of parameters are as in this example.

L=ApplyPAM("datatestL2.bin",k=5,init_method="BUILD",max_iter=1000,
           nthreads=0)
#> Reading symmetric distance/dissimilarity matrix datatestL2.bin
#> Package memuse is installed. OK.
#>   Memory used by the matrix: 48837 KiB, which is 0.01% of the available memory, which is 259192064 Kib.
#>   That seems OK.
#>   Matrix is a correct distance/dissimilarity matrix.
#> Starting BUILD initialization method, parallel version with 128 threads.
#> WARNING: all successive messages use R-numbering (from 1) for points and medoids. Substract 1 to get the internal C-numbers.
#> Looking for medoid 1. Medoid 1 found. Point 722. TD=8.823526
#> Looking for medoid 2. Medoid 2 found. Point 4473. 2434 reassigned points. TD=8.718881
#> Looking for medoid 3. Medoid 3 found. Point 2527. 1586 reassigned points. TD=8.663619
#> Looking for medoid 4. Medoid 4 found. Point 4504. 1175 reassigned points. TD=8.627773
#> Looking for medoid 5. Medoid 5 found. Point 3818. 932 reassigned points. TD=8.600748
#> Current TD: 8.600748
#> BUILD initialization method (parallel version) finished. Elapsed time: 0.0311963 s
#> Starting improved FastPAM1 method in parallel implementation with 128 threads.
#> WARNING: all successive messages use R-numbering (from 1) for points and medoids. Substract 1 to get the internal C-numbers.
#> Iteration 0.    Exiting, since DeltaTDst is 0.000000. Final value of TD is 8.600746
#> Optimization method (parallel version) finished. Elapsed time: 0.0061024 s
#> Time summary  (parallel implementation with 128 threads).
#>    Initalization: 0.0311963 s (method BUILD).
#>    Optimization:  0.0061024 s in 0 iterations.
#>    Total time:    0.0372987 s (0 minutes, 0.0372987 seconds).

Parameters init_method (and another optional parameter, initial_med) deserve special comment. The first is the method to initialize the medoids. Its possible values are BUILD, LAB and PREV. The rest of the algorithm make medoid swapping between the points of the initial set made with BUILD or LAB and the rest of points until no swap can reduce the objective function, which is the sum of distances of each point to its closest medoid. But this may fall (and indeed falls) in local minima. If you initialize with BUILD or LAB the optional parameter initial_med cannot be used.

The initialization methods BUILD and LAB are described in the paper from Schubert at al. (Schubert and Rousseeuw (2019)). BUILD is deterministic. LAB uses a sample of the total points to initialize. Obviously, you can run LAB to get different initializations and compare the results.

The returned object is a list with two fields: med and clasif. This will be explained later.

From now on, typical calls to obtain only the initial medoids would be

Lbuild=ApplyPAM("datatestL2.bin",k=5,init_method="BUILD",max_iter=0)
#> Reading symmetric distance/dissimilarity matrix datatestL2.bin
#> Package memuse is installed. OK.
#>   Memory used by the matrix: 48837 KiB, which is 0.01% of the available memory, which is 259254008 Kib.
#>   That seems OK.
#>   Matrix is a correct distance/dissimilarity matrix.
#> Starting BUILD initialization method, parallel version with 128 threads.
#> WARNING: all successive messages use R-numbering (from 1) for points and medoids. Substract 1 to get the internal C-numbers.
#> Looking for medoid 1. Medoid 1 found. Point 722. TD=8.823526
#> Looking for medoid 2. Medoid 2 found. Point 4473. 2434 reassigned points. TD=8.718881
#> Looking for medoid 3. Medoid 3 found. Point 2527. 1586 reassigned points. TD=8.663619
#> Looking for medoid 4. Medoid 4 found. Point 4504. 1175 reassigned points. TD=8.627773
#> Looking for medoid 5. Medoid 5 found. Point 3818. 932 reassigned points. TD=8.600748
#> Current TD: 8.600748
#> BUILD initialization method (parallel version) finished. Elapsed time: 0.0314978 s
#> Time summary  (parallel implementation with 128 threads).
#>    Initalization: 0.0314978 s (method BUILD).
#>    Optimization:  0 s in 0 iterations.
#>    Total time:    0.0314978 s (0 minutes, 0.0314978 seconds).
Llab1=ApplyPAM("datatestL2.bin",k=5,init_method="LAB",max_iter=0)
#> Reading symmetric distance/dissimilarity matrix datatestL2.bin
#> Package memuse is installed. OK.
#>   Memory used by the matrix: 48837 KiB, which is 0.01% of the available memory, which is 259201688 Kib.
#>   That seems OK.
#>   Matrix is a correct distance/dissimilarity matrix.
#> Starting LAB initialization method, serial version.
#> WARNING: all successive messages use R-numbering (from 1) for points and medoids. Substract 1 to get the internal C-numbers.
#> Looking for medoid 1. Medoid 1 found. Point 3388. TD=8.878892
#> Looking for medoid 2. Medoid 2 found. Point 3742. 2386 reassigned points. TD=8.772475
#> Looking for medoid 3. Medoid 3 found. Point 731. 1441 reassigned points. TD=8.723454
#> Looking for medoid 4. Medoid 4 found. Point 165. 1147 reassigned points. TD=8.688979
#> Looking for medoid 5. Medoid 5 found. Point 2823. 1150 reassigned points. TD=8.654373
#> Current TD: 8.654373
#> LAB initialization method (serial version) finished. Elapsed time: 0.00136817 s
#> Time summary  (parallel implementation with 128 threads).
#>    Initalization: 0.00136817 s (method LAB).
#>    Optimization:  0 s in 0 iterations.
#>    Total time:    0.00136817 s (0 minutes, 0.00136817 seconds).
Llab2=ApplyPAM("datatestL2.bin",k=5,init_method="LAB",max_iter=0)
#> Reading symmetric distance/dissimilarity matrix datatestL2.bin
#> Package memuse is installed. OK.
#>   Memory used by the matrix: 48837 KiB, which is 0.01% of the available memory, which is 259203044 Kib.
#>   That seems OK.
#>   Matrix is a correct distance/dissimilarity matrix.
#> Starting LAB initialization method, serial version.
#> WARNING: all successive messages use R-numbering (from 1) for points and medoids. Substract 1 to get the internal C-numbers.
#> Looking for medoid 1. Medoid 1 found. Point 1864. TD=8.886647
#> Looking for medoid 2. Medoid 2 found. Point 3389. 2528 reassigned points. TD=8.765924
#> Looking for medoid 3. Medoid 3 found. Point 4615. 1643 reassigned points. TD=8.707090
#> Looking for medoid 4. Medoid 4 found. Point 4934. 792 reassigned points. TD=8.685280
#> Looking for medoid 5. Medoid 5 found. Point 944. 1011 reassigned points. TD=8.655772
#> Current TD: 8.655772
#> LAB initialization method (serial version) finished. Elapsed time: 0.00131156 s
#> Time summary  (parallel implementation with 128 threads).
#>    Initalization: 0.00131156 s (method LAB).
#>    Optimization:  0 s in 0 iterations.
#>    Total time:    0.00131156 s (0 minutes, 0.00131156 seconds).

As it can be seen, to get and compare different initializations you must set the parameter max_iter to the value 0. In this case no iterations of objective function reduction are performed, and the returned object contains the initial medoids and the classification induced by them. Notice that even looking equal, the results of the latter two calls are different since LAB initializes with a random component (the sample to choose initial medoids is chosen randomly).

You can check that the medoids, stored in Llab1$med and Llab2$med (see more on this below) are in general different.

Now, these results can be used to initialize PAM if you find that any of them contains a specially good set of medoids. This is the role of method PREV that we have mentioned before. A typical call would be

Llab2Final=ApplyPAM("datatestL2.bin",k=5,init_method="PREV",
                    initial_med=Llab2$med)
#> Reading symmetric distance/dissimilarity matrix datatestL2.bin
#> Package memuse is installed. OK.
#>   Memory used by the matrix: 48837 KiB, which is 0.01% of the available memory, which is 259203044 Kib.
#>   That seems OK.
#>   Matrix is a correct distance/dissimilarity matrix.
#> Starting improved FastPAM1 method in parallel implementation with 128 threads.
#> WARNING: all successive messages use R-numbering (from 1) for points and medoids. Substract 1 to get the internal C-numbers.
#> Iteration 0. Medoid at place 4 (point 4934) swapped with point 723; TD-change=-0.024954; TD=8.630825. 1413 reassigned points.
#> Iteration 1. Medoid at place 5 (point 944) swapped with point 4473; TD-change=-0.011845; TD=8.618980. 1356 reassigned points.
#> Iteration 2. Medoid at place 3 (point 4615) swapped with point 2527; TD-change=-0.007798; TD=8.611182. 1271 reassigned points.
#> Iteration 3. Medoid at place 1 (point 1865) swapped with point 4504; TD-change=-0.005358; TD=8.605824. 1302 reassigned points.
#> Iteration 4. Medoid at place 2 (point 3389) swapped with point 3818; TD-change=-0.005081; TD=8.600743. 1200 reassigned points.
#> Iteration 5.    Exiting, since DeltaTDst is 0.000000. Final value of TD is 8.600743
#> Optimization method (parallel version) finished. Elapsed time: 0.0389671 s
#> Time summary  (parallel implementation with 128 threads).
#>    Initalization: 0 s (method PREV).
#>    Optimization:  0.0389671 s in 4 iterations (0.00974178 seconds/iteration).
#>    Total time:    0.0389671 s (0 minutes, 0.0389671 seconds).

The initial set of medoids is taken from the object returned by the former calls.

With respect to that object, as we said it is a list with two vectors. The first one, L$med, has as many components as requested medoids and the second, L$clasif, has as many components as instances.

Medoids are expressed in L$med by its number in the array of vectors (row number in the dissimilarity matrix) starting at 1 (R convention).

L$clasif contains the number of the medoid (i.e.: the cluster) to which each instance has been assigned, according to their order in L$med (also from 1).

This means that if L$clasif[p] is m, the point p belongs to the class grouped around medoid L$med[m]. Let us see it.

# Which are the indexes of the points chosen as medoids?
L$med
#>  V723 V4473 V2527 V4504 V3818 
#>   723  4473  2527  4504  3818
#
# In which class has point 147 been classified?
L$clasif[147]
#> V147 
#>    3
#
# And which is the index (row in the dissimilarity matrix)
# of the medoid closest to point 147?
L$med[L$clasif[147]]
#> V2527 
#>  2527

In this way, values in L$clasif are between 1 and the number of medoids, as we can see:

min(L$clasif)
#> [1] 1
max(L$clasif)
#> [1] 5

They can be used as factors.

Calculating silhouette

It is interesting to filter points based on the degree in which they belong to a cluster. Indeed, cluster refinement can be done getting rid of points far away from any cluster center, or which are at a similar distance of two or more of them.

This is characterized by the silhouette of each point. Three functions deal with this: CalculateSilhouette, FilterBySilhouetteQuantile and FilterBySilhouetteThreshold.

S=CalculateSilhouette(Llab2$clasif,"datatestL2.bin",nthreads=0)
#> Package memuse is installed. OK.
#>   Memory used by the matrix: 48837 KiB, which is 0.01% of the available memory, which is 259279492 Kib.
#>   That seems OK.
#>    Calculating silhouette (parallel version) with 128 threads.
#> 5000 points classified in 5 classes.
#> Finished parallel implementation of silhouette (including dissimilarity matrix load). Elapsed time: 0.0424387 s

The parameters to function CalculateSilhouette are the array of class membership, as returned by ApplyPAM in its clasif field, and the file with the matrix of dissimilarities.

A parallel implementation has been programmed, being nthreads as explained before.

Silhouette is a number in \([-1,1]\); the higher its value, the most centered a point is in its cluster.

The returned object S is a numeric vector with the value of the silhouette for each point, which will be a named vector if the classification vector was named.

This vector can be converted to an object of the class cluster:silhouette with the function NumSilToClusterSil (which needs the vector of classifications, too). This is done so that, if you load the package cluster (Maechler et al. (2022)), plot will generate the kind of silhouette plots included in such package.

If the package cluster is installed you can try to execute this: (Sorry, we can’t try ourselves since we don’t know if cluster is installed in your system and the CRAN check does not allow the use of installed.packages to test it)

Sclus <- NumSilToClusterSil(Llab2$clasif,S)
library(cluster)
plot(Sclus)

Probably the plot does not look very nice with this random data which yields a very bad clustering (since they are not, by its own nature, organized in clusters) but with real data you should see significant things (see package scellpam (Domingo (2023c))).

Once the silhouette is calculated we can filter it by quantile or by threshold. All points under this quantile or threshold will be discarded, except if they are medoids. Parameters are:

As an example,

Lfilt=FilterBySilhouetteQuantile(S,Llab2,"datatest.bin",
                                 "datatestFilt.bin","datatestL2.bin",
                                 "datatestL2Filt.bin",0.2)

If the original matrix contained row and column names, the column names are copied and the row names are transported for those rows that remain. The same happens with respect to rows of the dissimilarity matrix.

Notice that the new dissimilarity matrix could have been calculated from the matrix of filtered counts with CalcAndWriteDissimilarityMatrix but creating it here, simply getting rid of the filtered rows and columns is much faster.

Also, if a medoid is below the silhouette quantile, it will not be filtered out, but a warning message will be shown, since this is a strange situation that may indicate that some of your clusters are not real but artifacts due to a few outliers that are close to each other.

But remember that this was the result of the first step of the PAM algorithm, so probably you will want to make them iterate.

Lfinal=ApplyPAM("datatestL2Filt.bin",k=length(Lfilt$med),
                init_method="PREV",initial_med=Lfilt$med)
#> Reading symmetric distance/dissimilarity matrix datatestL2Filt.bin
#> Package memuse is installed. OK.
#>   Memory used by the matrix: 31257 KiB, which is 0.01% of the available memory, which is 259233972 Kib.
#>   That seems OK.
#>   Matrix is a correct distance/dissimilarity matrix.
#> Starting improved FastPAM1 method in parallel implementation with 128 threads.
#> WARNING: all successive messages use R-numbering (from 1) for points and medoids. Substract 1 to get the internal C-numbers.
#> Iteration 0. Medoid at place 4 (point 3950) swapped with point 3587; TD-change=-0.018652; TD=8.621079. 1035 reassigned points.
#> Iteration 1. Medoid at place 5 (point 752) swapped with point 575; TD-change=-0.010399; TD=8.610681. 985 reassigned points.
#> Iteration 2. Medoid at place 3 (point 3704) swapped with point 2010; TD-change=-0.003868; TD=8.606813. 1017 reassigned points.
#> Iteration 3. Medoid at place 2 (point 2701) swapped with point 3611; TD-change=-0.003396; TD=8.603418. 976 reassigned points.
#> Iteration 4. Medoid at place 1 (point 1478) swapped with point 3054; TD-change=-0.003293; TD=8.600125. 981 reassigned points.
#> Iteration 5.    Exiting, since DeltaTDst is 0.000000. Final value of TD is 8.600125
#> Optimization method (parallel version) finished. Elapsed time: 0.0376265 s
#> Time summary  (parallel implementation with 128 threads).
#>    Initalization: 0 s (method PREV).
#>    Optimization:  0.0376265 s in 4 iterations (0.00940663 seconds/iteration).
#>    Total time:    0.0376265 s (0 minutes, 0.0376265 seconds).

Of course, we might have used simply 5 as number of medoids, k, since this does not change by filtering, but this is to emphasize the fact that ApplyPAM with method PREV requires both parameters to be consistent.

Comparison with other PAM implementations

The user might want to compare this PAM implementation with others provided for instance by packages cluster (Maechler et al. (2022)) or ClusterR (Mouselimis (2023)). In cluster the input is either the data matrix (so the distance matrix is calculated inside the pam function) or directly the distance matrix but as a R vector with the lower-diagonal part of the symmetric matrix ordered by columns (i.e.: column 1 from M(2,1) to M(n,1), followed by column 2 from M(3,2) to M(n,2) and so on, up to M(n,n-1). This is a vector of \(n(n-1)/2\) components. To facilitate such comparison the function GetSubdiag is provided which takes as input the jmatrix file with the distance matrix and returns the vector to be passed to pam in the aforementioned packages.

d = GetSubdiag("datatestL2.bin")

Then, explicit comparison could be done with:

(Sorry, we can’t try ourselves since we don’t know if cluster is installed in your system and the CRAN check does not allow the use of installed.packages to test it)

library(cluster)
clusterpam = pam(d,diss=TRUE,k=5)
print(sort(clusterpam$id.med))
print(sort(L$med))

Similarly, you can check against the ClusterR package. In this package you need the complete dissimilarity matrix to be passed so we have to get it:

# Be patient, this may take some time...
Dm = GetJManyRows("datatestL2.bin",seq(1:nvec))

and then

library(ClusterR)
ClusterRpam = Cluster_Medoids(Dm,clusters=5)
print(sort(ClusterRpam$medoid_indices))
print(sort(L$med))

In all cases we tried with this simple (but random) example results were the same. In other cases with a large number of points some medoids were different in the different implementations but the value of the function to minimize (sum of distances) was always the same, indicating that they were equivalent minima. You can test this with function GetTD as follows:

TDparallelpam = GetTD(L,"datatestL2.bin")

# This is to adapt cluster package output format to ours, since this is what our GetTD function expects...
Lcl = list()
Lcl$med = clusterpam$id.med
Lcl$clasif = clusterpam$clustering
TDcluster = GetTD(Lcl,"datatestL2.bin")

# The same with ClusterR package:
LclR = list()
LclR$med = ClusterRpam$medoid_indices
LclR$clasif = ClusterRpam$clusters
TDClusterR = GetTD(LclR,"datatestL2.bin")

and see that variables TDparallelpam, TDcluster and TDClusterR are equal.

Domingo, Juan. 2023a. Jmatrix: Read from/Write to Disk Matrices with Any Data Type in a Binary Format.
———. 2023b. Parallelpam: Applies the Partitioning-Around-Medoids (PAM) Clustering Algorithm to Big Sets of Data Using Parallel Implementation, If Several Cores Are Available.
———. 2023c. Scellpam: Applying Partitioning Around Medoids to Single Cell with High Number of Cells.
Maechler, Martin, Peter Rousseeuw, Anja Struyf, Mia Hubert, and Kurt Hornik. 2022. Cluster: Cluster Analysis Basics and Extensions. https://CRAN.R-project.org/package=cluster.
Mouselimis, Lampros. 2023. ClusterR: Gaussian Mixture Models, k-Means, Mini-Batch-Kmeans, k-Medoids and Affinity Propagation Clustering. https://CRAN.R-project.org/package=ClusterR.
Schubert, Erich, and Peter J. Rousseeuw. 2019. “Faster k-Medoids Clustering: Improving the PAM, CLARA, and CLARANS Algorithms.” In Similarity Search and Applications, 171–87. Springer International Publishing. https://doi.org/10.1007/978-3-030-32047-8_16.