Uniform Concentration Bounds toward a Unified Framework for Robust Clustering
Year: 2021

Abstract:
Recent advances in center-based clustering continue to improve upon the drawbacks of Lloyd’s celebrated k-means algorithm over 60 years after its introduction. Various methods seek to address poor local minima, sensitivity to outliers, and data that are not well-suited to Euclidean measures of fit, but many are supported largely empirically. Moreover, combining such approaches in a piecemeal manner can result in ad hoc methods, and the limited theoretical results supporting each individual contribution may no longer hold. Toward addressing these issues in a principled way, this paper proposes a cohesive robust framework for center-based clustering under a general class of dissimilarity measures. In particular, we present a rigorous theoretical treatment within a Median-of-Means (MoM) estimation framework, showing that it subsumes several popular k-means variants. In addition to unifying existing methods, we derive uniform concentration bounds that complete their analyses, and bridge these results to the MoM framework via Dudley’s chaining arguments. Importantly, we neither require any assumptions on the distribution of the outlying observations nor on the relative number of observations n to features p. We establish strong consistency and an error rate of O(n−1/2) under mild conditions, surpassing the best-known results in the literature. The methods are empirically validated thoroughly on real and synthetic datasets.

heading: 1 Introduction
text: Clustering is a fundamental task in unsupervised learning, which seeks to discover naturally occurring groups within a dataset. Among a plethora of algorithms in the literature, center-based methods remain widely popular, with k-means [35, 33] as the most prominent example. Given n data points X = {Xi : i = 1, . . . , n} ⊂ Rp, k-means seeks to partition the data into k mutually exclusive and exhaustive groups by minimizing the within cluster variance. Representing the cluster centroids Θ = {θ1, . . . ,θk} ⊂ Rp, the k-means problem is formulated as the minimization of the objective function fk-means(Θ) = ∑n i=1 min1≤j≤k d(Xi,θj), where d(·, ·) is a dissimilarity measure. Taking the squared Euclidean distance d(x,y) = ‖x− y‖22 yields the classical k-means formulation. Clustering via k-means is well-documented to be sensitive to initialization [54], prone to getting stuck in poor local optima [57, 56, 17], and fragile against linearly non-separable clusters [39] or ∗Joint first authors contributed equally ‡ Corresponding author
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
in the presence of even a single outlying observation [26]. Researchers continue to tackle these drawbacks of k-means while preserving its simplicity and interpretability [4, 1, 48, 40, 58, 49]. Although often successful in practice, few of these approaches are grounded in rigorous theory or provide finite-sample statistical guarantees. The available consistency results are mostly asymptotic in nature [17, 51, 52], lacking convergence rates or operating under restrictive assumptions on the relation between p and n [42]. For example, the recent large sample analysis [42] of the hard Bregman k-means algorithm [4] obtains an asymptotic error rate of O( √ log n/n) under restrictive assumptions on the relation between n and p. The approach by [48] focuses on a unified framework from an optimization perspective, while [3, 49] focus on different divergence-based methods to better understand the underlying feature-space.
The presence of outliers only further complicates matters. Outlying data are common in real applications, but many of the aforementioned approaches are fragile to deviations from the assumed data generating mechanism. On the other hand, recent work on robust clustering methods [18, 21, 26] do not integrate the practical advances surveyed above, and similarly lack finite-sample guarantees. To bridge this gap, the Median of Means (MoM) literature provides a promising and attractive framework to robustify center-based clustering methods against outliers. MoM estimators are not only insensitive to outliers, but are also equipped with exponential concentration results under the mild condition of finite variance [34, 29, 6, 31, 27]. Recently, near-optimal results for mean estimation [37], classification [30], regression [36, 34], clustering [26, 13], bandits [14] and optimal transport [46] have been established from this perspective.
Under the MoM lens, we propose a unified framework for robust center-based clustering. Our treatment considers a family of Bregman loss functions not restricted to the classical squared Euclidean loss. We explore the exact sample error bounds for a general class of algorithms under mild regularity assumptions that apply to popular existing approaches, which we show to be special cases. The proposed framework allows for the data to be divided into two categories: the set of inliers (I) and the set of outliers (O). The inliers are assumed to be independent and identically distributed (i.i.d.), while absolutely no assumption is required on O, allowing outliers to be unboundedly large, dependent on each other, sampled from a heavy-tailed distribution and so on. Our contributions within the MoM framework make use of Rademacher complexities [5, 7] and symmetrization arguments [53], powerful tools that often find use in the empirical process literature but, in our view, are underexplored in the context of robust clustering.
The paper is organized as follows. In Section 2, we identify a general centroid-based clustering framework which encompasses k-means [35], k-harmonic means [57], and power k-means [56] as special cases, to name a few. We show how this framework is made robust via Median of Means estimation, yielding an array of center-based robust clustering methods. Within this framework, we derive uniform deviation bounds and concentration inequalities under standard regularity conditions through bounding Rademacher complexity by metric entropy via Dudley’s chaining argument in Section 3. The analysis newly reveals the convergence rate for popularly used clustering methods such as k-harmonic means and power k-means, matching the known rate results for k-means, and elegantly carries over to their MoM counterparts. We then implement and empirically assess the resulting algorithms through simulated and real data experiments. In particular, we find that power k-means [56] under the MoM paradigm outperforms the state-of-the-art in the presence of outliers.

heading: 1.1 Related theoretical analyses of clustering
text: In seminal work, Pollard [43], proved the strong consistency of k-means under a finite second moment assumption, spurring the large sample analysis of unsupervised clustering methods. This result has been extended to separable Hilbert spaces [10, 32] and for related algorithms [15, 51, 52], but these do not provide guarantees on the number of samples required so that the excess risk falls below a given threshold. Towards finding probabilistic error bounds, following research derived uniform concentration results for k-means and its variants [49], sub-Gaussian distortion bounds for the kmedians problem [12], and aO( √ log n/n) bound on k-means with Bregman divergences [41]. More recently, concentration inequalities for k-means under the MoM paradigm have been established [26, 13] under the restriction that sample cluster centroids (Θ̂n in Section 3) are assumed to be bounded. This paper shows how a number of center-based clustering methods can be brought under the same umbrella and can be robustified using a general-purpose scheme. The theoretical analyses of this broad spectrum of methods is conducted via Dudley’s chaining arguments and through the aid
of Rademacher complexity-based uniform concentration bounds. This approach enables us to replace assumptions on the sample cluster centroids by a bounded support assumption of the (inlying) data points, which yields an intuitive way of ensuring the boundedness of the cluster centroids through the obtuse angle property of Bregman divergences (Lemma 3.1 below). In contrast to the prior results, our analyses are not asymptotic in nature, so that the derived bounds hold for all values of the model parameters.

heading: 2 Problem Setting and Proposed Method
text: We consider the problem of partitioning a set of n data points X = {Xi : i = 1, . . . , n} ⊂ Rp into k mutually exclusive clusters. In a center-based clustering framework, we represent the jth cluster by its centroid θj ∈ Rp for each j ∈ {1, . . . , k}. To quantify the notion of “closeness", we allow the dissimilarity measure to be any Bregman divergence. Recall any differentiable, convex function φ : Rp → R generates the Bregman divergence dφ : Rp×Rp → R≥0 (R≥0 denoting the set of non-negative reals) defined as
dφ(x,y) = φ(x)− φ(y)− 〈∇φ(y),x− y〉.
For instance, φ(u) = ‖u‖22 generates the Euclidean distance. Without loss of generality, one may assume φ(0) = ∇φ(0) = 0. In this paradigm, clustering is achieved by minimizing the objective
1
n n∑ i=1 Ψα (dφ(Xi,θ1), . . . , dφ(Xi,θk)) := fΘ(X). (1)
Here Ψα : Rk≥0 → R≥0 is a component-wise non-decreasing function (such as a generalized mean) of the dissimilarities {dφ(X,Θj)}kj=1 which satisfies Ψ(0) = 0. The hyperparameter α ∈ A ⊆ Rq is specified by the user, and we will additionally assume that Ψ is Lipschitz. For intuition, we begin by showing how this setup includes several popular clustering methods. Examples: Suppose φ(u) = ‖u‖22 and Ψ(x) = ( ∑k j=1 x −1 j ) −1. Then the objective (1) becomes 1 n ∑n i=1( ∑k j=1 ‖Xi − θj‖ −2 2 ) −1, which is the objective function of k-harmonic means clustering [57]. Now consider other generalized means: take Ψs(x) = Ms(x) where we denote the power mean Ms(x) = (k −1∑k i=1 x s i )
1/s for s ∈ (−∞,−1]. Then objective (1) coincides with the recent power k-means method [56], 1n ∑n i=1Ms ( ‖Xi − θ1‖22, . . . , ‖Xi − θk‖22 ) . When Ψ(x) = min1≤j≤k xj ,
(1) recovers Bregman hard clustering 1n ∑n i=1 min1≤j≤k dφ(Xi,θj) proposed in [4] for any valid φ, while the special case of φ(u) = ‖u‖22 yields the familiar Euclidean k-means problem [35]. In what follows, we derive concentration bounds that establish new theoretical guarantees such as consistency and convergence rates for clustering algorithms in this framework. These analyses lead us to a unified, robust framework by embedding this class of methods within the Median of Means (MoM) estimation paradigm. Via elegant connections between the properties of MoM estimators and Vapnik-Chervonenkis (VC) theory, our MoM estimators too will inherit uniform concentration inequalities from the preceding analysis, extending convergence guarantees to the robust setting.
Median of Means Instead of directly minimizing the empirical risk (1), MoM begins by partitioning the data into L sets B1, . . . , BL which each contain exactly b many elements (discarding a few observations when n is not divisible by L). The partitions can be assigned uniformly at random, or can be shuffled throughout the algorithm [30]. MoM then entails a robust version of the estimator defined under (1) by instead minimizing the following objective with respect to Θ:
MoMnL(Θ) = Median
( 1
b ∑ i∈B1 fΘ(Xi), . . . , 1 b ∑ i∈BL fΘ(Xi)
) . (2)
Intuitively, MoM estimators are more robust than their ERM counterparts since under mild conditions, only a subset of the partitions is contaminated by outliers while others are outlier-free. Taking the median over partitions negates the influence of these spurious partitions, thus reducing the effect of outliers. Formal analysis of robustness via breakdown points is also available for MoM estimators [28, 44]; we will make use of the nice concentration properties of MoM estimators in Section 3.4.
Algorithm 1 MoM Clustering via Adagrad
Input: X ∈ Rn×p, k, L, fΘ(·), η, . Output: The cluster centroids Θ. Initialization: Randomly partition {1, . . . , n} into L many partitions of equal length. Randomly choose k points without replacement from X to initialize Θ0. repeat
Step 1: Find `t ∈ {1, . . . , L}, such that MoMnL(Θt) = 1b ∑ i∈B`t fΘt(Xi).
Step 2: g(t)j ← 1b ∑ i∈B`t
∇θjfΘ(Xi) Step 3: Update Θ by θ(t+1)j ← θ (t) j −
η√ + ∑(t) t′=1 ‖g (t′) j ‖22 g (t) j .
until objective (2) converges
Furthermore, optimizing (2) is made tractable via gradient-based methods. We advocate the Adagrad algorithm [19, 22], whose updates are given by
θ (t+1) j ← θ (t) j − η√ + ∑(t) t′=1 ‖g (t′) j ‖22 g (t) j
for hyperparameter > 0, learning rate η > 0, and g(t)j denoting a subgradient of MoM n L(Θ) at Θt. That is, if `t denotes the median partition at step t, then
∇Θj MoMnL(Θt) = 1
b ∑ i∈B`t ∇θjfΘ(Xi)|Θ=Θt .
For any Ψα differentiable, ∇θjfΘ(x) = ∂jΨα ( dφ(x,θ (t) 1 ), . . . , dφ(x,θ (t) k ) ) ∇θdφ(x,θ) by the chain rule. As a concrete illustration of the general formula, consider the power k-means objective fΘ(x) = Ms(‖x− θ1‖22, . . . , ‖x− θk‖22): upon partial differentiation, we obtain
∇θjfΘ(x) = 2
k 1 k k∑ j′=1 ‖x− θj′‖2s2 1/s−1 ‖x− θj‖2(s−1)2 (θj − x). As another example, the classic k-means objective fΘ(x) = minθ∈Θ ‖x − θ‖22 requires the subgradient∇θjfΘ(x) = 2(θj−x)1{j ∈ J (Θ,x)}, where J (Θ,x) = argmin1≤j≤k ‖x−θj‖22 and 1{·} denotes the indicator function. A pseudocode summary appears in Algorithm 1. This method tends to find clusters efficiently: for instance, we incur O(npk) per-iteration complexity applied to both the k-means and power k-means instances, which matches that of their original algorithms without robustifying via MoM [35, 33, 56]. Of course, here we update θj by an adaptive gradient step rather than by a closed form expression.
As emphasized by an anonymous reviewer, one should note that the algorithms under the proposed framework may converge to sub-optimal local solutions under the proposed optimization scheme. As is the case for k-means and its variants, this possibility arises due to the non-convexity of the objective functions (1) and (2). A complete theoretical understanding from an optimization perspective for such methods is not yet fully developed and global results are in general notoriously difficult to obtain. Having said that, it is worth noting that recent empirical analysis show that techniques such as annealing [56, 17] may be effective to circumvent this difficulty. Indeed, our experimental analysis (see Section 4) suggests that a robust version of power k-means using this annealing technique best overcomes this difficulty even in the presence of outliers.

heading: 3 Theoretical Analysis
text: Here we analyze properties of the proposed framework (1), with complete details of all proofs in the Supplement. Denoting M the set of all probability measures P with support on [−M,M ]p, i.e. P ([−M,M ]p) = 1, we first make standard assumptions that the data are i.i.d. with bounded components [9, 17, 42].
A1. X1 . . . ,Xn i.i.d.∼ P such that P ∈M.
Let Pn be the empirical distribution based of the data X1 . . . ,Xn. That is, for any Borel set A, Pn(A) = 1 n ∑n i=1 1{Xi ∈ A}. For notational simplicity, we write µf := ∫ fdµ. Appealing to the Strong Law of Large Numbers (SLLN) [2], PnfΘ a.s.−−→ PfΘ. Suppose Θ̂n be a (global) minimizer of (1) and Θ∗ be the global minimizer of PfΘ. Since the functions PnfΘ and PfΘ are close to each other as n become large, we can expect that their respective minimizers, Θ̂n and Θ∗ are also close to one another. To show that Θ̂n converges to Θ∗ as n → ∞, we consider bounding the uniform deviation supΘ |PnfΘ − PfΘ|. Towards establishing such bounds, we will posit two regularity assumptions on Ψα(·) and φ(·), beginning with a τα,k-Lipschitz condition on Ψα. A2. For all α ∈ A and any x,y ∈ Rk≥0, we have |Ψα(x)−Ψα(y)| ≤ τα,k‖x− y‖1.
We also assume a weaker form of a standard condition that the gradient∇φ(·) is Hp-Lipschitz [50]; unlike their work, note we do not additionally require strong convexity of φ. A3. There exists Hp ≥ 0 such that ‖∇φ(x)−∇φ(y)‖2 ≤ Hp‖x− y‖2 for all x,y ∈ [−M,M ]p.
These conditions are mild, and can be seen to hold for all of the aforementioned popular special cases. For instance, taking φ(u) = ‖u‖2 yields the squared Euclidean distances, and we see that ∇φ(u) = 2u so that A3 is satisfied with constant Hp = 2. Now, let Ψ(x) = min1≤j≤k xj as in classical k-means, and denote j∗ ∈ argmin1≤j≤k yj . For any vectors x,y ∈ R p ≥0,
Ψ(x)−Ψ(y) = min 1≤j≤k xj − min 1≤j≤k yj = min 1≤j≤k xj − yj∗ ≤ xj∗ − yj∗ ≤ ‖x− y‖1.
Thus, Ψ is clearly non-negative and componentwise non-decreasing, and satisfies A2. Similarly, if we take Ψs(x) = Ms(x), the conditions are met for power k-means: it is again non-negative and non-decreasing in its components, and satisfies A2 with τα,k = k−1/s due to [8].

heading: 3.1 Bounds on Θ̂n and Θ∗
text: Toward proving that Θ̂n converges to Θ∗, we will need to show that both Θ̂n and Θ∗ lie in [−M,M ]k×p. To this end, we first establish the obtuse angle property for Bregman divergences. Lemma 3.1. Let C be a convex set and suppose PC(θ) be the projection of θ onto C, with respect to the Bregman divergence dφ(·, ·), i.e. PC(θ) = arg minx∈C dφ(x,θ) (assuming it exists). Then,
dφ(x,θ) ≥ dφ(x, PC(θ)) + dφ(PC(θ),θ), for all x ∈ C.
We next show that to minimize PnfΘ or PfΘ, it is enough to restrict the search to [−M,M ]k×p. Lemma 3.2. Let A2 hold, and Q ∈M. Let dφ : Rp×Rp → R≥0 be a Bregman divergence. Then for any Θ ∈ Rk×p, there exists Θ′ ∈ [−M,M ]k×p, such that QfΘ′ ≤ QfΘ.
Since we can restrict our attention to [−M,M ]k×p to minimize QfΘ, we have the following: Corollary 3.1. Let Q ∈ M and dφ : Rp×Rp → R≥0 be a Bregman divergence. If Θ0 = arg minΘ∈Rk×p ∫ fΘdQ, then Θ0 ∈ [−M,M ]k×p.
Now note under A1 both P and Pn have support contained in [−M,M ]k×p. The following corollary is thus implied by replacing Q by P and Pn. Corollary 3.2. Under A1—3, both Θ̂n,Θ∗ ∈ [−M,M ]k×p.
Now that we have bounded Θ̂n and Θ∗ in a compact set, the following section supplies probabilistic bounds on the uniform deviation 2 supΘ∈[−M,M ]k×p |PnfΘ − PfΘ| via metric entropy arguments.

heading: 3.2 Concentration Inequality and Metric Entropy Bounds via Rademacher Complexity
text: We have proven that Θ̂n,Θ∗ ∈ [−M,M ]k×p. To bound the difference |PfΘ̂n −PfΘ∗ |, we observe |PfΘ̂n − PfΘ∗ | = PfΘ̂n − PfΘ∗ = PfΘ̂n − PnfΘ̂n + PnfΘ̂n − PnfΘ∗ + PnfΘ∗ − PfΘ∗
≤ PfΘ̂n − PnfΘ̂n + PnfΘ∗ − PfΘ∗ ≤ 2 sup Θ∈[−M,M ]k×p |PnfΘ − PfΘ|. (3)
Thus, to bound |PfΘ̂n − PfΘ∗ |, it is enough to prove bounds on supΘ∈[−M,M ]k×p |PnfΘ − PfΘ|. The main idea here is to bound this uniform deviation via Rademacher complexity, and in turn bound the Rademacher complexity itself [20, 38]. Let F = {fΘ : Θ ∈ [−M,M ]k×p}, and denote the F-norm [2] between two probability measures µ and ν as ‖µ− ν‖F = supf∈F | ∫ fdµ− ∫ fdν|. We recall the definition of Rademacher complexity and covering number as follows. Definition 1. Let i’s be i.i.d. Rademacher random variables independent of X = {X1, . . . ,Xn}, i.e. P( i = 1) = P( i = −1) = 0.5, The population Rademacher complexity of F is defined as Rn(F) = E supf∈F 1n ∑n i=1 if(Xi) , where the expectation is over both and X . Definition 2. (δ-cover and Covering Number) Let (X, d) be a metric space. The set Xδ ⊆ X is said to be a δ-cover of X if for all x ∈ X , ∃ x′ ∈ Xδ, such that d(x, x′) ≤ δ. The δ-covering number of X w.r.t. d, denoted by N(δ;X, d), is the size of the smallest δ-cover of X with respect to d.
The following Lemma gives a bound for the δ-covering number of F under the supremum norm. The main idea here is to use the Lipschitz property of fΘ and then to find a cover of the search space for Θ, i.e. [−M,M ]k×p. This then automatically transcends to a cover of F under the sup-norm. Lemma 3.3. Let N(δ;F , ‖ · ‖∞) be the δ-covering number of F under ‖ · ‖∞. Then, under A1—3,
N(δ;F , ‖ · ‖∞) ≤ ( max {⌊ 8M2τα,kHpkp
δ
⌋ , 1 })kp .
To bound the Rademacher complexity, we will also need to bound the diameter of F under ‖ · ‖∞: Lemma 3.4. Let diam(F) = supf,g∈F ‖f − g‖∞. Then, under A1—3, diam(F) ≤ 8τα,kHpM2kp.
We are now ready to state the bound on the Rademacher complexity Rn(F) in Theorem 3.1. We provide a sketch of the argument below, with the complete proof details available in the Supplement. Theorem 3.1. Under A1—3,Rn(F) ≤ 48 √ πτα,kHpM 2(kp)3/2n−1/2.
Proof Sketch The main idea of the proof is to use Dudley’s chaining argument [20, 55] to bound the Rademacher complexity in terms of an integral involving the metric entropy logN(δ;F , ‖ · ‖∞). Let ∆ = 8HpM2k1−1/sp. Using the chaining approach, one can show that
Rn(F) ≤ 12√ n ∫ ∆ 0 √ logN( ;F , ‖ · ‖∞)d ≤ 12√ n ∫ ∆ 0 √ kp log ( max { ∆ , 1 }) d (4)
= 12√ n ∫ ∆ 0 √ kp log ( ∆ ) d = 48 √ πτα,kHpM 2(kp)3/2n−1/2.
Here the second inequality (4) follows from Lemma 3.3. Before we proceed, we state the following Lemma, which puts an uniform bound on ‖f‖∞, f ∈ F . Lemma 3.5. For all x ∈ [−M,M ]p and Θ ∈ [−M,M ]k×p,
0 ≤ Ψα(dφ(x,θ1), . . . , dφ(x,θk)) ≤ 4τα,kHpM2pk.
Now that we have established a bound on the Rademacher complexityRn(F), we are now ready to establish a uniform concentration inequality on ‖Pn − P‖F , supf∈F |Pnf − Pf | in the following Theorem, proven in the Supplement. Theorem 3.2. Under A1—3, with probability at least 1− δ, the following holds.
‖Pn − P‖F ≤ 96 √ πτα,kHpM 2(kp)3/2n−1/2 + 4τα,kHpM 2pk
√ log(2/δ)
2n .
The result in Theorem 3.2, reveals a non-asymptotic bound on |PfΘ̂n − PfΘ∗ |: Corollary 3.3. Under A1—3, with probability at least 1− δ, the following holds.
|PfΘ̂n − PfΘ∗ | ≤ 192 √ πτα,kHpM 2(kp)3/2n−1/2 + 8τα,kHpM 2pk
√ log(2/δ)
2n .
Proof. From equation (3), we know that |PfΘ̂n − PfΘ∗ | ≤ 2‖Pn − P‖F . The corollary follows immediately by application of Theorem 3.2.
3.3 Inference for Fixed p: Strong and √ n Consistency
We now discuss the classical domain where p is kept fixed [43, 17, 51] and show that the results above imply strong and √ n-consistency, mirroring some known results for existing methods such as k-means. We first solidify the notion of convergence of the centroids Θ̂n to Θ∗, following Pollard [43]. Since centroids are unique only up to label permutations, our notion of dissimilarity
diss(Θ1,Θ2) = min M∈Pk ‖Θ1 −MΘ2‖F
is considered over Pk the set of all k×k real permutation matrices, where ‖·‖F denotes the Frobenius norm. Now, we say that the sequence Θn converges to Θ if limn→∞ diss(Θn,Θ) = 0. We begin by imposing the following standard identifiablity condition [43, 51, 17] on P for our analysis. A4. For all η > 0, there exists > 0, such that PfΘ > PfΘ∗ + whenever diss(Θ,Θ∗) > η.
We now investigate the strong consistency properties of Θ̂n, and also investigate the rate at which |PfΘ̂n − PfΘ∗ | converges to 0. Theorem 3.3 states that indeed strong consistency holds, with convergence rate O(n−1/2). Note that this rate is faster than that found previously in [42]. Before we proceed, recall we say that Xn = OP (an) if the sequence Xn/an is tight [2]. Theorem 3.3. (Strong consistency and √ n-consistency) If p is kept fixed then under A1—4, Θ̂n
a.s.−−→ Θ∗ under P . Moreover, |PfΘ̂n − PfΘ∗ | = OP (n −1/2).
Note: Here An is OP (an) means that {An/an}n∈N is stochastically bounded or tight.

heading: 3.4 Inference Under the MoM Framework
text: Theorem 3.3 and the bounds we present above already establish novel statistical results that pertain to methods under (1) such as power k-means in the “uncontaminated" setting. In this section, we now extend these findings to the MoM setup under (2) in order to understand their behavior in the presence of outliers. Recall that in this setting, the data are partitioned into L equally sized blocks; without loss of generality, we take n = L · b. We denote the set of all inliers by {Xi}i∈I and outliers by {Xi}i∈O. Now let Θ̂ (MoM) n denote the minimizer of (2): towards establishing the error rate at which |Pf Θ̂
(MoM) n
− PfΘ∗ | goes to 0, we assume the following:

heading: A5. {Xi}i∈I
text: i.i.d.∼ P with P ∈M. A6. There exists η > 0 such that L > (2 + η)|O|.
We remark A5 is identical to A1, but imposed only on the inlying observations. A6 ensures at least half of the L partitions are free of outliers; note this is much weaker than requiring L > 4|O| as is done in recent work [30]. Importantly, we emphasize that no distributional assumptions regarding the outlying observations are made, allowing them to be unbounded, generated from heavy-tailed distributions, or dependent among each other. Proofs of the following results appear in the Supplement.
As a point of departure, we first establish that Θ̂ (MoM) n ∈ [−M,M ]k×p. Lemma 3.6. Let A2—3 and A5—6 hold. Then for any Θ ∈ Rk×p, there exists Θ′ ∈ [−M,M ]k×p, such that MoMnL(Θ ′) ≤ MoMnL(Θ).
Again, we may restrict the search space for finding Θ̂ (MoM) n in [−M,M ]k×p due to Lemma 3.6. Corollary 3.4. Under A2—3 and A5—6, Θ̂ (MoM) n ∈ [−M,M ]k×p.
Similarly to Section 3, we derive a uniform bound on supΘ∈[−M,M ]k×p |MoMnL(fΘ)− PfΘ| and then bound |Pf
Θ̂ (MoM) n − PfΘ∗ | in turn. For brevity we define δ := 2/(4 + η)− |O|/L, and use the notation “.” to suppress the absolute constants. The uniform deviation bound is as follows, with complete proof appearing in the Supplement.
Theorem 3.4. Under A2—3 and A5—6, with probability at least 1− 2e−2Lδ2 , the following holds.
sup Θ∈[−M,M ]k×p
|MoMnL(fΘ)− PfΘ| . τα,kHp max
{ kp √ L
n , (kp)3/2 √ |I| n } .
We give a brief outline of our proof for Theorem 3.4 below, with details appearing in the Supplement. Proof sketch We begin by noting that if supΘ∈[−M,M ]k×p ∑L `=1 1 {(P − PB`)fΘ > } > L 2 , then supΘ∈[−M,M ]k×p(PfΘ −MoMnL(fΘ)) > , where PB` denotes the empirical distribution of {Xi}i∈B` . This implies it is suffices to bound the quantity
P ( sup
Θ∈[−M,M ]k×p L∑ `=1 1 {(P − PB`)fΘ > } > L 2
) .
Introducing the function ϕ(t) = (t − 1)1{1 ≤ t ≤ 2} + 1{t > 2}, we begin by bounding outlier-free partitions, making use of the inequalities 1{t ≥ 2} ≤ ϕ(t) ≤ 1{t ≥ 1}. We then proceed by bounding supΘ∈[−M,M ]k×p ∑L `=1 1 {(P − PB`)fΘ > } by
the sum ξ1 + ξ2 + |O|, where ξ1 = supΘ∈[−M,M ]k×p ∑ `∈L Eϕ ( 2(P−PB` )fΘ ) and ξ2 =
supΘ∈[−M,M ]k×p ∑ `∈L [ ϕ ( 2(P−PB` )fΘ ) − Eϕ ( 2(P−PB` )fΘ )] . We bound ξ1 by appealing to
Hoeffding’s inequality, while we show that ξ2 can be bounded by applying the bounded difference inequality together with the result of Theorem 3.1 to bound the resulting Rademacher complexity.
The corollary below follows from this uniform bound, giving a non-asymptotic control over the difference |Pf
Θ̂ (MoM) n
− PfΘ∗ | in terms of the model parameters.
Corollary 3.5. Under A2—3 and A5—6, with probability at least 1− 2e−2Lδ2 , the following holds.
|Pf Θ̂
(MoM) n
− PfΘ∗ | . τα,kHp max { kp √ L
n , (kp)3/2 √ |I| n } .

heading: 3.5 Inference for Fixed k and p Under the MoM Framework
text: We now focus our attention back to the classical setting where the numbers of clusters and features remain fixed. To show Θ̂ (MoM) n is consistent for Θ
∗, we need to impose conditions such that the RHS of the bounds presented in Corollary 3.5 decrease to 0 as n→∞. We state the required conditions as follows.
A7. The number of partitions L = o(n), and L→∞ as n→∞.
These conditions are natural: as n grows, so too must L in order to maintain a proportion of outlierfree partitions. On the other hand, L must grow slowly relative to n to ensure each partition can be assigned sufficient numbers of datapoints. We note that A7 implies |O| = o(n), an intuitive and standard condition [30, 46, 41] as outliers should be few by definition.
In the following corollary, we focus on the squared Euclidean distance for center-based clustering under the MoM framework. We show that the (global) estimates obtained from MoM k-means, MoM power k-means etc. are consistent. We stress that the obtained convergence rate, as a function of n, L and |I|, does not depend on the choice of Ψα(·), as long as it satisfies A2, i.e. is Lipschitz. In particular, replacing Ψα(·) with min1≤j≤k xj , Ms(x) and ( ∑k j=1 x −1 j ) −1, the rate in Corollary 3.6 applies to robustified MoM versions of k-means, power k-means and k-harmonic means alike.
Corollary 3.6. Suppose φ(u) = ‖u‖22 and Ψα(·) satisfy A2. Then under A2—3 and A5—7,
|Pf Θ̂
(MoM) n
− PfΘ∗ | = OP ( max { L1/2n−1/2, n−1 √ |I| })
and Pf Θ̂
(MoM) n
P−→ PfΘ∗ . Moreover, whenever A4 additionally holds, we have Θ̂ (MoM) n P−→ Θ∗.
Remark These results imply that any MoM center-based algorithm under our paradigm admits a convergence rate of O ( max { L1/2n−1/2, n−1 √ |I| }) , when equipped with squared Euclidean
distance. Note that as L ≥ 1, max { L1/2n−1/2, n−1 √ |I| } = Ω ( n−1/2 ) . Thus, the convergence
rates for MoM variants in our framework are generally slower than their ERM counterparts, for which the rate is O(n−1/2). This is unsurprising as MoM operates on outlier-contaminated data; there is “no free lunch" in trading off robustness for rate of convergence. However, if the number of partitions L grows slowly relative to n (say, L = O(log n) so that |O| = O(log n)), then the convergence rates for MoM estimation become comparable to the ERM counterparts at Õ(n−1/2).

heading: 4 Empirical Studies and Performance
text: To validate and assess our proposed framework, we now turn to an empirical comparison of the proposed and peer clustering methods. We evaluate clustering quality under the Adjusted Rand Index (ARI) [23], with values ranging between 0 and 1 and 1 denoting a perfect match with the ground truth. Though it is not feasible to exhaustively survey center-based clustering methods, we consider a broad range of competitors, comparing to k-means [35], Partition Around Medoids (PAM) [25], k-medians [24], Robust kmeans++ (RKMpp) [18], Robust Continuous Clustering (RCC) [45], Bregman Trimmed k-means (BTKM) [21], MoM k-means (MOMKM) [26] and a novel MoM variant of power k-means (MOMPKM) implied under the proposed framework.
We consider two thorough simulated experiments below, with additional simulations and largescale real data results in the Supplement. While the extended comparisons are omitted for space considerations, they convey the same trends as the studies below. In all settings, we generate data in p = 5 dimensions, varying the number of clusters and the outlier percentages. True centers are spaced uniformly on a grid with θk = k−110 and observations are drawn from Gaussians around their ground truth centers with variance 0.1. Because we generate Gaussian data, we focus here on the Euclidean case, and do not consider other Bregman divergences in the present empirical study. In all experiments, we take the number of partitions L to be roughly double the number of outliers, and set the hyperparameters η and α to be 1.02 and 1 respectively by default. Results are averaged over 20 random restarts, and all competitors are initialized and tuned according to the standard implementations described in their original papers.
Experiment 1: increasing the number of clusters The first experiment assesses performance as the number of true clusters grows, while keeping the proportion of outliers fixed at 25%. Datasets are generated with the number of clusters k varying between 3 and 100. For each setting, we create balanced clusters drawing 30 points from each true center. The 25% outliers are generated from a uniform distribution with support on the range of the inlying observations. We repeat this data generation process 30 times under each parameter setting.
The average ARI values at convergence, plotted against the number of clusters along with error bars (± standard deviations), are shown in the left panel of Figure 1. We see that the robustified version of power k-means implied by our framework (labeled MOMPKM) achieves the best performance here. This may be unsurprising as the recent power k-means method was shown to significantly reduce sensitivity to local minima (which tend to increase with k), while MoM further protects the algorithm from outlier influence.
Experiment 2: increasing outlier percentage Following the same data generation process in Experiment 1, we now fix k = 20 while varying the outlier percentage from 0% to 50%. For each parameter setting, we again replicate the experiment 30 times. Average ARI values comparing the inlying observations to their ground truths are plotted in the right panel of Figure 1. Not only does this study convey similar trends as Experiment 1, but we see that competing methods continue to deteriorate with increasing outliers as one might expect, while MOMPKM remarkably remains relatively stable.
We see that the ERM-based methods such as BTKM, MOMKM, and PAM struggle when there is large number of clusters. Similarly, methods such as k-medians and RKMpp often stop short at poor local optima, which is quickly exacerbated by outliers despite initialization via clever seeding techniques. These phenomena are consistent with what has been reported in the literature for their non-robust counterparts [56, 17]. Overall, the empirical study suggests that a robust version of power k-means clustering under the MoM framework shows promise to handle several data challenges at once, in line with our theoretical analysis.

heading: 5 Discussion
text: In this paper, we proposed a paradigm for center-based clustering that unifies a suite of center-based clustering algorithms. Under this view, developed a simple yet efficient framework for robust versions of such algorithms by appealing to the Median of Means (MoM) philosophy. Using gradient-based methods, the MoM objectives can be solved with the same per-iteration complexity of Lloyd’s k-means algorithm, largely retaining its simplicity. Importantly, we derive a thorough analysis of the statistical properties and convergence rates by establishing uniform concentration bounds under i.i.d. sampling of the data. These novel theoretical contributions demonstrate how arguments utilizing Rademacher complexities and Dudley’s chaining arguments can be leveraged in the robust clustering context. As a result, we are able to obtain error rates that do not require asymptotic assumptions, nor restrictions on the relation between n and p. These findings recover asymptotic results such as strong consistency and √ n-consistency under classical assumptions.
As shown in the paper, the robustness of MoM estimators comes at the cost of slower convergence rates compared to their ERM counterparts. We emphasize that there is no “median-of-means magic", and that the efficacy of MoM depends on the interplay between the partitions and the outliers. If the number of partitions circumvents the impact of the outliers, the performance of MoM clustering estimates under our framework scales with the block size b as 1/ √ b = √ L/n. Since L can be
chosen to be approximately 2|O|, the obtained error rate is roughly O( √ |O|/n). If |O| scales
proportionately with n, however, the error bound of O( √ |O|/n) becomes meaningless. For our consistency results to hold, it is crucial that |O| = o(n), which in turn allows us to choose L that satisfies condition A7. If |O| = O(nβ), for some 0 < β < 1, the error rate is O(n(β−1)/2). This suggests possible future research directions in improving these ERM rates via finding so called “fast rates" under additional assumptions [11, 55]. Moreover, it will be fruitful to extend the results for noise distributions that satisfy only moment conditions. Recent works in convex clustering [47, 16] have considered sub-Gaussian models to obtain error rates. The work by [10] and recent work by [26] obtain similar error rates under finite second-moment conditions using assumptions that the cluster centroids Θ̂n are bounded, and it may be possible to extend our approach using local Rademacher complexities [7] to relax the bounded support assumption. One can also seek lower bounds on the approximation error or can explore high-dimensional robust center-based clustering under the proposed paradigm. Finally, we have not explored the implementation and empirical performance of Bregman versions of our MoM estimator, for instance with application to data arising from mixtures of exponential families other than the Gaussian case. These interesting directions remain open for future work.

Review summaries:
The authors consider the minimization of a Median-of-Means (MoM) objective, which is suitable to yield minimizers of a centroid-based clustering which is more robust to outliers. The type of the centroid-based clustering is defined over the employed measurement of distances, which can be any  Bregman divergence. The considered MoM objective is hence suitable to encompass various $k$-means formalizations.  The authors formulate an Adagrad method to minimize the MoM objective, and analyze properties of the global minimizer of the MoM objective with respect to robustness to outliers. They apply concentration inequalities to compare the MoM global minimizer to the global minimizer of the expected value of the generalized distance of points to their centroid. Since outliers do not influence the expected value, a small distance of the MoM minimizer to the one of the expected value indicates robustness to outliers.   The authors provide an experimental analysis of the power $k$-means objective, optimized by the proposed optimization scheme on synthetic data and one real-world dataset.

The paper provides a novel method for center-based clustering of i.i.d. points in $\mathbb{R}^d$ under a Bregman divergence. The key idea of the framework is to employ a median of means estimator defined over the objective function, i.e., the data set is partitioned into sets of equal size and the median of the objective values of the elements of the partition is minimized to obtain the clustering. The aim of the paper is to develop a clustering scheme that is more robust against outliers than common $k$-means and its variants/extensions and admits theoretical guarantees. To compute the clustering the well-known Adagrad gradient descent algorithm is used.  The main contribution of the paper is a rigorous theoretical analysis, proving consistency and convergence, i.e., the minimizer of the sample approaches the minimizer of the population as the sample size tends to infinity, at a bounded rate of order $(\sqrt{n})^{-1}$. The efficacy of the novel clustering method is experimentally evaluated and compared to other state of the art clustering algorithms. 

The paper looks at the Median-of-Means estimation framework for clustering and shows that it covers a range of k-means variants. They attempt to unify different frameworks of center-based clustering and introduce an algorithm that uses an adaptive gradient-based MoM algorithm. Under i.i.d. sampling of the data, they provide the convergence rates and the error rates obtained do not require any restrictions on the relation between the dimensions of the input points and the number of input points.