A penalized criterion for selecting the number of clusters for K-medians

In what follows, we consider a random variable X that takes values in $\mathbb{R}^{d}$ for some $d\geq 1$. It is well-known that the standard mean of X is not robust to corruptions. Hence, the median is preferred to the mean in robust statistics. The geometric median $m$, also called $L^{1}$-median or spatial median, of a random variable $X\in\mathbb{R}^{d}$ is defined by Haldane, (1948) as follows:

For the 1-dimensional case, the geometric median coincides with the usual median in $\mathbb{R}$. As Euclidean space $\mathbb{R}^{d}$ is strictly convex, the geometric median $m$ exists and is unique if the points are not concentrated around a straight line (Kemperman,, 1987). The geometric median is known to be robust and has a breakdown point of $0.5$.

Let us now consider a sequence of i.i.d copies $X_{1},...,X_{n}$ of $X$. In this paper, we focus on two methods to determine the geometric median. The first one is iterative and consists in considering the fix point estimates (Weiszfeld,, 1937; Vardi and Zhang,, 2000)

with a initial point $\hat{m}_{0}\in\mathbb{R}^{d}$ chosen arbitrarily such that it does not coincide with any of the $X_{i}$ and $\mathcal{X}_{t}=\{i,X_{i}\neq\hat{m}_{t}\}$. This Weiszfeld algorithm can be a flexible technique, but there are many implementation difficulties for massive data in high-dimensional spaces.

An alternative and simple estimation algorithm which can be seen as a stochastic gradient algorithm (Robbins and Monro,, 1951; Ruppert,, 1985; Duflo,, 1997; Cardot et al.,, 2013) and is defined as follows

where $m_{0}$ is an arbitrarily chosen starting point and $\gamma_{j}$ is a step size such that $\forall j\geq 1,\gamma_{j}>0,\sum_{j\geq 1}\gamma_{j}=\infty$ and $\sum_{j\geq 1}\gamma_{j}^{2}<\infty$. Its averaged version (ASG), which is effective for large samples of high dimension data, introduced by Polyak and Juditsky, (1992) and adapted by Cardot et al., (2013), is defined by

One can speak about averaging since $\overline{m}_{j}=\frac{1}{j}\sum_{i=1}^{j}m_{i}$. We note that, under suitable assumptions, both $\hat{m}_{t}$ and $\overline{m}_{t}$ are asymptotically efficient (Vardi and Zhang,, 2000; Cardot et al.,, 2013).

For a positive integer k, a vector quantizer Q of dimension d and codebook size k is a (measurable) mapping of the d-dimensional Euclidean $\mathbb{R}^{d}$ into a finite set of points $\left\{c_{1},...,c_{k}\right\}$ (Linder,, 2000). More precisely, the points $c_{i}\in\mathbb{R}^{d},i=1,...,k$ are called the codepoints and the vector composed of the code points $\left\{c_{1},...,c_{k}\right\}$ is called codebook, denoted by c. Given a d-dimensional random vector X admitting a finite first order moment, the $L^{1}$-distortion of a vector quantizer Q with codebook $c=\left\{c_{1},...,c_{k}\right\}$ is defined by

Let us now consider $X_{1},...,X_{n}$ random vectors $\in\mathbb{R}^{d}$ i.i.d with the same law as X. Then, one can define the empirical $L^{1}$-distortion as :

In this paper, we consider two types of K-medians algorithms : sequential and non sequential algorithm. The non sequential algorithm uses Lloyd-style iteration which alternates between an expectation (E) and maximization (M) step and is precisely described in Algorithm 1:

For $1\leq j\leq k,m_{j}$ is nothing but the geometric median of the points in the cluster $C_{j}$. As $m_{j}$ is not explicit, we will use Weiszfeld (indicated by “Offline”) or ASG (indicated by “Semi-Online”) to estimate it. The Online K-median algorithm proposed by Cardot et al., (2012) based on an averaged Robbins-Monro procedure (Robbins and Monro,, 1951; Polyak and Juditsky,, 1992) is described in Algorithm 2:

The non-sequential algorithms are effective but the computational time is huge compared to the sequential (“Online”) algorithm, which is very fast and only requires $O(knd)$ operations, where n is the sample size, k is the number of clusters and d is dimension. Furthermore, in case of large samples, Online algorithm is expected to estimate the centers of the clusters as well as the non-sequential algorithm Cardot et al., (2012). Then, in case of large sample size, Online algorithm should be preferred and vice versa.

In Section 3, we proved that the penalty shape is $a\sqrt{\frac{k}{n}}$ where $a$ is a constant to calibrate. To find the constant $a$, we will use the data-based calibration algorithm for penalization procedures that is explained at the end of section 3. This data-driven slope estimation method is implemented in CAPUSHE (CAlibrating Penalty Using Slope HEuristics) (Brault et al.,, 2011) which is available in the R package capushe https://cran.r-project.org/package=capushe. This proposed slope estimation method is made to be robust in order to preserve the eventual undesirable variations of criteria. More precisely, for a certain number of clusters $k$, the algorithm may be trapped by a local minima, which could create a “bad point” for the slope heuristic. The slope heuristic has therefore been designed to be robust to the presence of such points.

In what follows, we consider a random variable $X$ following a Gaussian Mixture Model with $k=6$ classes where the mixture density function is defined as

with, $\pi_{j}=\frac{1}{k}\hskip 8.53581pt\forall 1\leq j\leq k,\hskip 8.53581pt\mu_{j}\sim\mathcal{U}_{10}$ where $\mathcal{U}_{10}$ is the uniform law on the sphere of radius 10 and,

In what follows, we consider $n=3000$ i.i.d realizations of $X$ and $d=5$. We first focus on some visualization of our slope method.

To estimate $a\approx 2\hat{S}$ in the penalty function, it is sufficient to estimate $\hat{S}$, which is the slope of the red curve in Figure 1. As shown in Figure 1, the regression slope is estimated using the last 21 points, as it behaves like an affine function when $k$ is large. In Figure 2 (left), two possible elbows are observed in the curve. Consequently, the elbow method suggests considering either 5 or 6 as the number of clusters. We would prefer to choose 5 since the elbow point at 5 is more pronounced compared to the one at 6. Therefore, this method is not ideal in this case.

Figures 3 to 5 represent the data as curves, which we call “profiles” (the x-label corresponds to the coordinates, and the y-label to the values of the coordinates), gathered by cluster and with the centers of the groups represented in red. We also show the first two principal components of the data using robust principal component analysis components (RPCA) (Cardot and Godichon-Baggioni,, 2017). In Figure 3, we focus on the clustering obtained with the K-medians algorithm (“Offline” version) for non contaminated data. In each cluster, the curves are close to each other and also close to the median, and the profiles differ from one cluster to another, meaning that our method separated well the 6 groups. In order to visualize the robustness of the proposed method, we considered contaminated data with the law $Z=(Z_{1},...,Z_{5})$ where $Z_{i}$ are i.i.d, with $Z_{i}\sim\mathcal{T}_{1}$ where $\mathcal{T}_{1}$ is a Student law with one degree of freedom. Applying our method for selecting the number of clusters for K-medians algorithms, we selected the corrected number of clusters. Furthermore, the obtained groups, despite the presence of some outliers in each cluster, are coherent. Nevertheless, in the case of K-means clustering, the method found non homogeneous clusters, i.e. the method assimilates some far outliers as single clusters (see Figure 5). It’s important to note that, in the case of contaminated data (Figures 4 and 5), we only represented $95\%$ of the data to better visualize them. Then, in Figure, 5, Clusters 5, 7, 8, 11 and 12 are not visible since they are “far” outliers.

In what follows, we focus on the choice of the number of clusters and compare our results with different methods. For this, we generated some basic data sets in three different scenarios (see Fischer, (2011)) :
(S1) $4$ clusters in dimension $3$ : The data are generated by Gaussian mixture centered at $(0,0,0)$, $(0,2,3)$, $(3,0,-1)$, and $(-3,-1,0)$ with variance equal to the identity matrix. Each cluster contains 500 data points.
(S2) $5$ clusters in dimension $4$ : The data are generated by Gaussian mixture centered at $(0,0,0,0)$, $(3,5,-1,0)$, $(-5,0,0,0)$, $(1,1,6,-2)$ and $(1,-3,-2,5)$ with variance equal to the identity matrix. Each cluster contains 500 data points.
(S3) $3$ clusters in dimension $2$ : The data are generated by a Student mixture centered at $(0,0)$, $(0,6)$ and $(5,3)$ with 2 degree of freedom. Each cluster contains 500 data points.

We applied three different methods for determining the number of clusters : the proposed slope method, Gap Statistic and Silhouette method. For each method, we use four clustering algorithms : K-medians (“Online”, “Semi-Online”, “Offline”) and K-means. For each scenario, we contaminated our data with the law $Z=(Z_{1},...,Z_{d})$ where $Z_{i}$ are i.i.d, with $Z_{i}\sim\mathcal{T}_{1}$ where $\mathcal{T}_{1}$ is a Student law with 1 degree of freedom. Then, we evaluate our method for the different methods and scenarios by considering:

• •

  N : The number of times we get the correct value of cluster in 50 repeated trials without contaminated data.

• •

  $\bar{k}$ : The average of number of clusters obtained over 50 trials without contaminated data.

• •

  $N_{0.1}$ : The number of times we get the correct value of cluster in 50 repeated trials with $10\%$ of contaminated data.

• •

  $\bar{k}_{0.1}$ : The average of number of clusters obtained over 50 trials with $10\%$ of contaminated data.

In case of well separated clusters as in the scenario (S2), the gap statistics method and silhouette method give competitive results. Nevertheless, for closer clusters, the slope method works much better than gap statistics and silhouette method as in the scenario (S1). The gap statistics method only works in scenario 2 and is ineffective in the presence of contamination. In closer cluster scenarios, it often predicts $1$ as the number of clusters. The silhouette method performs moderately well in scenario 2 and very well in scenario 3, but it is globally not as competitive as the slope method, especially in cases of contaminated data. In scenarios 1 and 2 with slope method, Offline, Semi-Online, Online and K-means give better results but in cases of contamination, K-means crashes completely while the other three methods seem to be not too much sensitive. Furthermore, on non-Gaussian data (scenario 3), the K-means method does not work at all. In such cases, K-median clustering is often preferred over K-means clustering.

Overall, in every scenario, Offline, Semi-Online, Online K-medians with the slope method give very competitive results and in the case where the data are contaminated, they clearly outperform other methods, especially the Offline method.

We now focus on the impact of contaminated data on the selection of the number of clusters in K-medians clustering, particularly in higher dimensions. We compare our method with Gap Statistic and SigClust (Liu et al.,, 2008; Huang et al.,, 2015; Bello et al.,, 2023) in the Offline setting, as it yields competitive results, as noted in the previous section. Concerning SigClust, it is a method which enables to test whether a sample comes from a single Gaussian or several in high dimension. Then, starting from $k=k_{0}$, we test for all possible pairs of clusters whether the fusion of the two clusters comes from a single Gaussian or not. If the test rejected the hypothesis that the combined cluster is a single Gaussian for all fittings, the same procedure is repeated for $k+1$. If there is a fitting for which the test is not rejected, it is considered that these two clusters should be merged, and the procedure is stopped. The optimal number of clusters is then determined as $k_{\text{opt}}=k-1$. It is important to note that we did not compare with Gap Statistics, as it is computationally expensive, especially in high dimensions.

In this aim, we generate data using a Gaussian mixture model with 10 classes in 100 and 200 dimensions, where the centers of the classes are randomly generated on a sphere with radius $10$, and each class contains 100 data points. The data is contaminated with the law $Z=(Z_{1},...,Z_{d})$, where $d$ is the dimension (100 or 200 for each scenario), $Z_{i}$ are i.i.d, with two possible scenarios:

1. 1.

   $Z_{i}\sim\mathcal{T}_{1}$,

2. 2.

   $Z_{i}\sim\mathcal{T}_{2}$.

Here, $\mathcal{T}_{m}$ is the Student law with m degrees of freedom. In what follows, let us denote by $\rho$ the proportion of contaminated data. In order to compare the different clustering results, we focus on the Adjusted Rand Index (ARI) (Rand,, 1971; Hubert and Arabie,, 1985) which is a measure of similarity between two clusterings and which relies on taking into account the right number of correctly classified pairs. We evaluate, for each scenario, the average number of clusters obtained over 50 trials and the average ARI evaluated only on uncontaminated data.

We observe that in the case of non-contamination, we obtain similar results across all methods. However, in the presence of contamination, our method consistently performs well, while others struggle to identify an appropriate number of clusters. With a Student distribution contamination of 1 degree of freedom, our method excels in terms of both the number of clusters and the ARI. The results with a Student distribution contamination of 2 degrees of freedom are comparable to those obtained using the Silhouette method.

In summary, our method demonstrates remarkable robustness in the face of contaminated data, making it a strong choice for clustering in higher dimensions. The comparison with the Silhouette, Gap Statistic, and SigClust in the offline setting reaffirms the effectiveness of our approach, especially when computational efficiency is a critical factor in high-dimensional data.

We will first briefly discuss the data we used for clustering, which was provided by Califrais, a company specializing in developing environmentally responsible technology to optimize logistics flows on a large scale. Our goal is to build a Recommender System that is designed to suggest items individually for each user based on their historical data or preferences. In this scenario, the clustering algorithms can be employed to identify groups of similar customers where each group consists of customers share similar features or properties. It is crucial to perform robust clustering in order to develop an effective Recommender System.

The dataset includes information on 508 customers, including nine features that represent the total number of products purchased in each of the following categories: Fruits, Vegetables, Dairy products, Seafood, Butcher, Deli, Catering, Grocery, and Accessories and equipment. Therefore, we have a sample size of $n=508$ and a dimensionality of $d=9$. To apply clustering, we will determine the appropriate number of clusters using the proposed method. Before applying our method, we normalize our data using RobustScaler. This removes the median and scales the data according to the Interquartile Range, which is the range between the 1st quartile and the 3rd quartile.

We plotted the profiles of the clusters obtained using our Slope method, Silhouette and Gap Statistic in Figures 6, 7, and 8. We observe that our method indicates 5 clusters, while the Gap Statistic suggests 3 clusters, and Silhouette suggests 2 clusters. Regarding the Silhouette method, the second cluster obtained is not homogeneous, as seen in Figure 7. We obtain 3 clusters with the Gap Statistic method. The important thing to note is that the Gap Statistic method separates the second cluster obtained by Silhouette into two clusters (Cluster 2 and Cluster 3). However, in the third cluster of Gap Statistic (Figure 8), homogeneity is still not achieved. In Figure 6, it can be seen that the clusters generated by our slope method are more or less homogeneous. To establish a connection with the simulations conducted in Section 4.2, for example, in scenario (S1), we observed that Silhouette and Gap Statistics failed to find the correct number of clusters when the clusters are closer. This is reflected here, as the behavior of clients does not change significantly, resulting in close clusters. To provide an overview of our clusters, the first cluster represents customers who regularly consume products from all categories. The third cluster consists of customers who frequently engage with catering products. Clusters 2, 4, and 5 correspond to customers who consume significant amounts of Butcher, Deli, and Catering products at different levels, as depicted in the figure 6.

The proposed penalized criterion, calibrated with the help of the slope heuristic method, consistently gives competitive results for selecting the number of clusters in K-medians, even in the presence of outliers, outperforming other methods such as Gap Statistics, Silhouette, and SigClust. Notably, our method demonstrates excellent performance even in high dimensions. Among the three K-medians algorithms, Offline, Semi-Online, and Online, their performances are generally analogous, with Offline being slightly better. However, for large sample sizes, one may prefer the Online K-medians algorithm in terms of computation time. As discussed in Section 2, it is recommended to use the Offline algorithm for moderate sample sizes, the Semi-Online algorithm for medium sample sizes, and the Online algorithm for large sample sizes. In our real-life data illustration, our proposed method consistently produces more robust clusters and a more suitable number of clusters compared to other methods.

In conclusion, our paper presents a robust and efficient approach for selecting clusters in K-medians, demonstrating superior performance even in challenging scenarios. The findings provide practical recommendations for algorithm selection based on sample size, reinforcing the applicability of our proposed method in real-world clustering scenarios.

The authors wish to thank Califrais for providing the real-life data and Raphaël Carpintero Perez for the data preprocessing work.

First, we provide some definitions and lemmas that are useful to prove Theorems 3.1 and 3.2.

Definitions :

• •

  Let $(S,p)$ be a totally bounded metric space. For any $F\subset S$ and $\epsilon>0$ the $\epsilon$-covering number $N_{p}(F,\epsilon)$ of F is defined as the minimum number of closed balls with radius $\epsilon$ whose union covers $F$.

• •

  Let $(S,p)$ be a totally bounded metric space. For any $F\subset S$,
  $\text{diam}(F)=\sup\left\{p(x,y):x,y\in F\right\}$.

• •

  A Family $\left\{T_{s}:s\in S\right\}$ of zero-mean random variables indexed by the metric space (S, p) is called subgaussian in the metric p if for any $\lambda>0$ and $s,s^{\prime}\in S$ we have

  $$\mathbb{E}\left[e^{\lambda(T_{s}-T_{s^{\prime}})}\right]\leq e^{\frac{\lambda^{2}p(s,s^{\prime})^{2}}{2}}.$$

• •

  The Family $\left\{T_{s}:s\in S\right\}$ is called sample continuous if for any sequence $s_{1},s_{2}..,\in S$ such that $s_{j}\rightarrow s\in S$ we have $T_{s_{j}}\rightarrow T_{s}$ with probability one.

The proof of the Theorem 3.1 is inspired by the proof of Theorem 3 in Linder, (2000).

Theorem 3.2 is an adaptation of Theorem 8.1 in Massart, (2007) and Theorem 2.1 in Fischer, (2011).