Mini-batch kernel $k$-means

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 12:17:04

Abstract: We present the first mini-batch kernel $k$-means algorithm. Our algorithm achieves an order of magnitude improvement in running time compared to the full batch algorithm, with only a minor negative effect on the quality of the solution. Specifically, a single iteration of our algorithm requires only $O(n(k+b))$ time, compared to $O(n^2)$ for the full batch kernel $k$-means, where $n$ is the size of the dataset and $b$ is the batch size. 
  
  We provide a theoretical analysis for our algorithm with an early stopping condition and show that if the batch is of size $\Omega((\gamma / \epsilon)^2\log (n\gamma/\epsilon))$, the algorithm must terminate within $O(\gamma^2/\epsilon)$ iterations with high probability, where $\gamma$ is the bound on the norm of points in the dataset in feature space,
and $\epsilon$ is a threshold parameter for termination. 
Our results hold for any reasonable initialization of centers. When the algorithm is initialized with the $k$-means++ initialization scheme, it
achieves an approximation ratio of $O(\log k)$.

Many popular kernels are normalized (e.g., Gaussian, Laplacian), which implies $\gamma=1$. For these kernels, taking $\epsilon$ to be a constant and $b=\Theta(\log n)$, our algorithm terminates within $O(1)$ iterations where each iteration takes time $O(n(\log n+k))$.

Main paper: 
 introduction: 
 Mini-batch methods are among the most successful tools for handling huge datasets for machine learning. Notable examples include Stochastic Gradient Descent (SGD) and mini-batch k-means [27] . Mini-batch k-means [27] is one of the most popular clustering algorithms used in practice [21] .While k-means is widely used due to it's simplicity and fast running time, it requires the data to be linearly separable to achieve meaningful clustering. Unfortunately, many real-world datasets do not have this property. One way to overcome this problem is to project the data into a high, even infinite, dimensional space (where it is hopefully linearly separable) and run k-means on the projected data using the "kernel-trick".Kernel k-means achieves significantly better clustering compared to k-means in practice. However, its running time is considerably slower. Surprisingly, prior to our work there was no attempt to speed up kernel k-means using a mini-batch approach.Problem statement We are given an input (dataset), X = {x i } n i=1 , of size n and a parameter k representing the number of clusters. A kernel for X is a function K : X × X → R that can be realized by inner products. That is, there exists a Hilbert space H and a map ϕ : X → H such that ∀x, y ∈ X, ⟨ϕ(x), ϕ(y)⟩ = K(x, y). We call H the feature space and ϕ the feature map.In kernel k-means the input is a dataset X and a kernel function K as above. Our goal is to find a set C of k centers (elements in H) such that the following goal function is minimized:1 n x∈X min c∈C ∥c -ϕ(x)∥ 2 .Equivalently we may ask for a partition of X into k parts, keeping C implicit. 1Lloyd's algorithm The most popular algorithm for (non kernel) k-means is Lloyd's algorithm, often referred to as the k-means algorithm [17] . It works by randomly initializing a set of k centers and performing the following two steps: (1) Assign every point in X to the center closest to it. (2) Update every center to be the mean of the points assigned to it. The algorithm terminates when no point is reassigned to a new center. This algorithm is extremely fast in practice but has a worst-case exponential running time [3, 30] .Mini-batch k-means To update the centers, Lloyd's algorithm must go over the entire input at every iteration. This can be computationally expensive when the input data is extremely large. To tackle this, the mini-batch k-means method was introduced by [27] . It is similar to Lloyd's algorithm except that steps (1) and ( 2) are performed on a batch of b elements sampled uniformly at random with repetitions, and in step (2) the centers are updated slightly differently. Specifically, every center is updated to be the weighted average of its current value and the mean of the points (in the batch) assigned to it. The parameter by which we weigh these values is called the learning rate, and its value differs between centers and iterations. The larger the learning rate, the more a center will drift towards the new batch cluster mean.Lloyd's algorithm in feature space Implementing Lloyd's algorithm in feature space is challenging as we cannot explicitly keep the set of centers C. Luckily, we can use the kernel function together with the fact that centers are always set to be the mean of cluster points to compute the distance from any point x ∈ X in feature space to any center c = 1 |A| y∈A ϕ(y) as follows:∥ϕ(x) -c∥ 2 = ⟨ϕ(x) -c, ϕ(x) -c⟩ = ⟨ϕ(x), ϕ(x)⟩ -2⟨ϕ(x), c⟩ + ⟨c, c⟩ = ⟨ϕ(x), ϕ(x)⟩ -2⟨ϕ(x), 1 |A| y∈A ϕ(y)⟩ + ⟨ 1 |A| y∈A ϕ(y), 1 |A| y∈A ϕ(y)⟩,where A can be any subset of the input X. While the above can be computed using only kernel evaluations, it makes the update step significantly more costly than standard k-means. Specifically, the complexity of the above may be quadratic in n [9] .Mini-batch kernel k-means Applying the mini-batch approach for kernel k-means is even more difficult because the assumption that cluster centers are always the mean of some subset of X in feature space no longer holds.In Section 4 we present our algorithm and derive a recursive expression that allows us to compute the distances of all points to current cluster centers (in feature space). Our algorithm implements this by updating a data structure that maintains the inner products between the data and centers in feature space. This means the running time of each iteration of our algorithm is only O(n(b + k)) compared to O(n2 ) for the full-batch algorithm.In Section 5 we go on to provide theoretical guarantees for our algorithm. This is somewhat tricky for mini-batch algorithms due to their stochastic nature, as they may not even converge to a local-minima.To overcome this hurdle, we take the approach of [26] and answer the question: how long does it take mini-batch kernel k-means to terminate with an early stopping condition. Specifically, we terminate the algorithm when the improvement on the batch drops below some user provided parameter, ϵ.Early stopping conditions are very common in practice (e.g., sklearn [21] ).Applying the k-means++ initialization scheme for our initial centers implies we achieve the same approximation ratio, O(log k) in expectation, as the full-batch algorithm. The approximation guarantee of k-means++ is guaranteed already in the initialization phase (Theorem 3.1 in [4] ), and the execution of Lloyd's algorithm following initialization can only improve the solution. We show that w.h.p 2 the global goal function is decreasing throughout our execution which implies that the approximation guarantee remains the same.While our general approach is similar to [26] , we must deal with the fact that H may have an infinite dimension. The guarantees of [26] depend on the dimension of the space in which k-means is executed, which is unacceptable in our case. We overcome this by parameterizing our results by a new parameter γ = max x∈X ∥ϕ(x)∥. We note that for normalized kernels, such as the popular Gaussian and Laplacian kernels, it holds that γ = 1. We show that if the batch size is Ω((γ/ϵ) 2 log(nγ/ϵ)) then w.h.p. our algorithm terminates in O(γ 2 /ϵ) iterations. Our theoretical results are summarised in Theorem 1 (where Algorithm 1 is presented in Section 4). Theorem 1. The following holds for Algorithm 1:1. Each iteration takes O(n(b + k)) time,2. If b = Ω((γ/ϵ) 2 log(nγ/ϵ)) then it terminates in O(γ 2 /ϵ) iterations with high probability, 3. When initialized with k-means++ it achieve a O(log k) approximation ratio in expectation.Our result improves upon [26] significantly when a normalized kernel is used since Theorem 1 doesn't depend on the input dimension. Our algorithm copes better with non linearly separable data and requires a smaller batch size ( O(1/ϵ 2 ) vs O((d/ϵ) 2 )))3 for normalized kernels. This is particularly apparent with high dimensional datasets such as MNIST [16] where the dimension squared is already nearly ten times the number of datapoints.The learning rate we use, suggested in [26] , differs from the standard learning rate of sklearn in that it does not go to 0 over time. Unfortunately, this new learning rate is non-standard and [26] did not present experiments comparing their learning rate to that of sklearn.In Section 6 we evaluate our results experimentally both with the learning rate of [26] and that of sklearn. We also fill the experimental gap left in [26] by evaluating (non-kernel) mini-batch k-means with their new learning rate compared to that of sklearn. To allow a fair empirical comparison, we run each algorithm for a fixed number of iterations without stopping conditions. Our results are as follows:• Mini-batch kernel k-means is significantly faster than full-batch kernel k-means, while achieving solutions of similar quality, which are superior to the non-kernel version.• The learning rate of [26] results in solutions with better quality both for mini-batch kernel k-means and mini-batch k-means. 
related work: 
 Until recently, mini-batch k-means was only considered with a learning rate going to 0 over time. This was true both in theory [29, 27] and practice [21] . Recently, [26] proposed a new learning which does not go to 0 over time, and showed that if the batch is of size Ω((d/ϵ) 2 ), mini-batch k-means must terminate within O(d/ϵ) iterations with high probability, where d is the dimension of the input, and ϵ is a threshold parameter for termination.A popular approach to deal with the slow running time of kernel k-means is constructing a coreset of the data. A coreset for kernel k-means is a weighted subset of X with the guarantee that the solution quality on the coreset is close to that on the entire dataset up to a (1 + ϵ) multiplicative factor. There has been a long line of work on coresets for k-means an kernel k-means [25, 10, 5] , and the current state-of-the-art for kernel k-means is due to [13] . They present a coreset algorithm with a nearly linear (in n and k) construction time which outputs a coreset of size poly(kϵ -1 ).In [7] the authors only compute the kernel matrix for uniformly sampled set of m points from X. Then they optimize a variant of kernel k-means where the centers are constrained to be linear combinations of the sampled points. The authors do no provide worst case guarantees for the running time or approximation of their algorithm.Another approach to speed up kernel k-means is by computing an approximation for the kernel matrix. This can be done by computing a low dimensional approximation for ϕ (without computing ϕ explicitly) [23, 8, 6] , or by computing a low rank approximation for the kernel matrix [19, 31] .Kernel sparsification techniques construct sparse approximations of the full kernel matrix in subquadratic time. For smooth kernel functions such as the polynomial kernel, [22] presents an algorithm for constructing a (1 + ϵ)-spectral sparsifier for the full kernel matrix with a nearly linear number of non-zero entries in nearly linear time. For the gaussian kernel, [18] show how to construct a weaker, cluster preserving sparsifier using a nearly linear number of kernel density estimation querries.We note that our results are complementary to coresets, dimensionality reduction, and kernel sparsification, in the sense that we can compose our method with these techniques.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: From our experimental results, each iteration of our mini-batch kernel k-means algorithm is still slower than that of mini-batch kmeans. This difference in speed may be alleviated by lazily enumerating the kernel matrix for very large datasets.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
preliminaries: 
 Throughout ℓ i=1 g(Y [i]). We borrow the following notation from [14] and generalize it to Hilbert spaces. For every x, y ∈ H let ∆(x, y) = ∥x -y∥ 2 . We slightly abuse notation and and also write ∆(x, y) = ∥ϕ(x) -ϕ(y)∥ 2 when x, y ∈ X and ∆(x, y) = ∥ϕ(x) -y∥ 2 when x ∈ X, y ∈ H (similarly when x ∈ H, y ∈ X). For every finite tuple S ⊆ X and a vector x ∈ H let ∆(S, x) = y∈S ∆(y, x). Let us denote γ = max x∈X ∥ϕ(x)∥. Let us define for any finite tuple S ⊆ X the center of mass of the tuple ascm(S) = 1 |S| x∈S ϕ(x).We now state the kernel k-means problem using the above notation.Kernel k-means We are given an input X = (x i ) n i=1 and a parameter k. Our goal is to (implicitly) find a tuple C ⊆ H of k centers such that the following goal function is minimized:1 n x∈X min C∈C ∆(x, C).Let us define for every x ∈ X the function f x : H k → R where f x (C) = min C∈C ∆(x, C). We can treat H k as the set of k-tuples of vectors in H. We also define the following function for every tupleA = (a i ) ℓ i=1 ⊆ X: f A (C) = 1 ℓ ℓ i=1 f ai (C).Note that f X is our original goal function. We extensive use of the notion of convex combination: Definition 2. We say that y ∈ H is a convex combination of X if y = x∈X p x ϕ(x), such that ∀x ∈ X, p x ≥ 0 and x∈X p x = 1.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes] Justification: Proofs are complete and correct to the best of our knowledge.Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and cross-referenced.• All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
our algorithm: 
 We present our pseudo-code as Algorithm 1. It requires an initial set of cluster centers such that every center is a convex combination of X. This guarantees that all subsequent centers are also a convex combination of X. Note that if we initialize the centers using the kernel version of k-means++, this is indeed the case.Algorithm 1 proceeds by repeatedly sampling a batch of size b (the batch size is a parameter). For the i-th batch the algorithm (implicitly) updates the centers using the learning rate α i j for center j. Note that the learning rate may take on different values for different centers, and may change between iterations. Finally, the algorithm terminates when the progress on the batch is below ϵ, a user provided parameter. While our termination guarantees (Section 5) require a specific learning rate, it does not affect the running time of a single iteration, and we leave it as a parameter for now.Recursive distance update rule While for (non kernel) k-means the center updates and assignment of points to clusters is straightforward, this is tricky for kernel k-means and even harder for mini-batch kernel k-means. Specifically, how do we overcome the challenge that we do not maintain the centers explicitly?To assign points to centers in the (i + 1)-th iteration, it is sufficient to know ∥ϕ(x) -C j i+1 ∥ 2 for every j. If we can keep track of this quantity through the execution of the algorithm, we are done.Input:• Dataset X = (x i ) n i=1 , batch size b, early stopping parameter ϵ • Initial centers (C j 1 ) k j=1 where C j 1 is a convex combination of X for all j ∈ [k] for i = 1 to ∞ do Sample b elements, B i = (y 1 , . . . , y b), uniformly at random from X (with repetitions)for j = 1 to k do B j i = x ∈ B i | arg min ℓ∈[k] ∆(x, C ℓ i ) = jα j i is the learning rate for the j-th cluster for iteration iC j i+1 = (1 -α j i )C j i + α j i cm(B j i ) if f Bi (C i+1 ) -f Bi (C i ) < ϵ then Return C i+1Let us derive a recursive expression for the distances∥ϕ(x) -C j i+1 ∥ 2 = ⟨ϕ(x), ϕ(x)⟩ -2⟨ϕ(x), C j i+1 ⟩ + ⟨C j i+1 , C j i+1 ⟩. We first expand ⟨ϕ(x), C j i+1 ⟩, ⟨ϕ(x), C j i+1 ⟩ = ⟨ϕ(x), (1 -α j i )C j i + α j i cm(B j i )⟩ = (1 -α j i )⟨ϕ(x), C j i ⟩ + α j i ⟨ϕ(x), cm(B j i )⟩.Then we expand⟨C j i+1 , C j i+1 ⟩, ⟨C j i+1 , C j i+1 ⟩ = ⟨(1 -α j i )C j i + α j i cm(B j i ), (1 -α j i )C j i + α j i cm(B j i )⟩ = (1 -α j i ) 2 ⟨C j i , C j i ⟩ + 2α j i (1 -α j i )⟨C j i , cm(B j i )⟩ + (α j i ) 2 ⟨cm(B j i ), cm(B j i )⟩.Assuming that ⟨C j i , C j i ⟩ and ⟨ϕ(x), C j i ⟩ are known for all j ∈ [k] and for all x ∈ X, we can compute ⟨C j i+1 , C j i+1 ⟩ and ⟨ϕ(x), C j i+1 ⟩ for all j ∈ [k] and x ∈ X, which implies we can compute the distances from any point in the batch to all centers. We now bound the running time of a single iteration of the outer loop in Algorithm 1. Let us denote b j i = B j i and recall that cm(B j i ) = 1 b j i y∈B j i ϕ(y). Therefore, computing ⟨ϕ(x), cm(B j i )⟩ = 1 b j i y∈B j i ⟨ϕ(x), ϕ(y)⟩ requires O(b j i ) time. Similarly, computing ⟨cm(B j i ), cm(B j i )⟩ requires O((b j i ) 2 ) time.Let us now bound the time it requires to compute ⟨ϕ(x), C j i+1 ⟩ and ⟨C j i+1 , C j i+1 ⟩.Assuming we know ⟨ϕ(x), C j i ⟩ and ⟨C j i , C j i ⟩, updating ⟨ϕ(x), C j i+1 ⟩ for all x ∈ X, j ∈ [k] requires O(n(b + k)) time. Specifically, the ⟨ϕ(x), C j i ⟩ term is already known from the previous iteration and we need to computeα j i ⟨ϕ(x), cm(B j i )⟩ for every x ∈ X, j ∈ [k] which requires n j∈[k] b j i = nb time. Finally, updating ⟨ϕ(x), C j i+1 ⟩ for all x ∈ X, j ∈ [k] requires O(nk) time. Updating ⟨C j i+1 , C j i+1 ⟩ requires O(b 2 + kb) time. Specifically, ⟨C j i , C j i ⟩ is known from the previous iteration and computing ⟨cm(B j i ), cm(B j i )⟩ for all j ∈ [k] requires O( j∈[k] (b j i ) 2 ) = O(b 2 ) time. Computing ⟨C j i , cm(B j i )⟩ for all j ∈ [k]requires time O(b) using ⟨ϕ(x), C j i ⟩ from the previous iteration. Therefore, the total running time of the update step (assigning points to new centers) is O(n(b + k)). To perform the update at the (i + 1)-th step we only need ⟨ϕ(x), C j i ⟩, ⟨C j i , C j i ⟩, which results in a space complexity of O(nk). This completes the first claim of Theorem 1.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: Our main claims are theoretical; however, we provide all the information needed to replicate the experimental results.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.• Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closedsource models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
termination guarantee: 
 Section preliminaries We introduce the following definitions and lemmas to aid our proof of the second claim of Theorem 1. Lemma 3. For every y which is a convex combination of X it holds that ∥y∥ ≤ γ.Proof. The proof follows by a simple application of the triangle inequality:∥y∥ = ∥ x∈X p x ϕ(x)∥ ≤ x∈X ∥p x ϕ(x)∥ = x∈X p x ∥ϕ(x)∥ ≤ x∈X p x γ = γ. Lemma 4. For any tuple of k centers C ⊂ H d which are a convex combination of points in X, it holds that ∀A ⊆ X, f A (C) ≤ 4γ 2 .Proof. It is sufficient to upper bound f x . Combining that fact that every C ∈ C is a convex combination of X with the triangle inequality, we have that∀x ∈ X, f x (C) ≤ max C∈C ∆(x, C) = ∆(x, y∈X p y ϕ(y)) = ∥ϕ(x) - y∈X p y ϕ(y)∥ 2 ≤ (∥ϕ(x)∥ + ∥ y∈X p y ϕ(y)∥) 2 ≤ 4γ 2 .We state the following simplified version of an Azuma bound for Hilbert space valued martingales from [20] , followed by a standard Hoeffding bound. Theorem 5 ([20] ). Let H be a Hilbert space and let Y 0 , ..., Y m be a H-valued martingale, such that∀1 ≤ i ≤ m, ∥Y i -Y i-1 ∥ ≤ a i . It holds that P r[∥Y m -Y 0 ∥ ≥ δ] ≤ e Θ δ 2 m i=1 a 2 i . Theorem 6 ([12]). Let Y 1 , ..., Y m be independent random variables such that ∀1 ≤ i ≤ m, E[Y i ] = µ and Y i ∈ [a min , a max ]. Then P r 1 m m i=1 Y k -µ ≥ δ ≤ 2e -2mδ 2 /(amax-amin) 2 .The following lemma provide concentration guarantees when sampling a batch. Lemma 7. Let B be a tuple of b elements chosen uniformly at random from X with repetitions. For any fixed tuple of k centers, C ⊆ H which are a convex combination of X, it holds that:P r[|f B (C) -f X (C)| ≥ δ] ≤ 2e -bδ 2 /2γ 2 .Proof. Let us write B = (y 1 , . . . , y b ), where y i is a random element selected uniformly at random from X with repetitions. For every such y i define the random variable Z i = f yi (C). These new random variables are IID for any fixed C. It also holds that ∀i ∈[b], E[Z i ] = 1 n x∈X f x (C) = f X (C) and that f B (C) = 1 b x∈B f x (C) = 1 b b i=1 Z i .Applying the Hoeffding bound (Theorem 6) with parameters m = b, µ = f X (C), a max -a min ≤ 4γ 2 (due to Lemma 4) we get that:P r[|f B (C) -f X (C)| ≥ δ] ≤ 2e -bδ 2 /2γ 2 .For any tuple S ⊆ X and some tuple of cluster centers C = (C ℓ ) ℓ∈[k] ⊂ H, C implies a partition (S ℓ ) ℓ∈[k] ofk j=1 ∆(S j , C j ) ≤ k j=1 ∆(S j , C j ).Recall that C j i is the j-th center in the beginning of the i-th iteration of Algorithm 1 and(B ℓ i ) ℓ∈[k] is the partition of B i induced by C i . Let (X ℓ i ) ℓ∈[k]be the partition of X induced by C i . We now have the tools to analyze Algorithm 1 with the learning rate of [26] . Specifically, we assume that the algorithm executes for at least t iterations, the learning rate is α j i = b j i /b, where b j i = B j i , and the batch size is b = Ω((γ/ϵ) 2 log(nt)). We show that the algorithm must terminate within t = O(γ/ϵ) steps w.h.p. Plugging t back into b, we get that a batch size of b = Ω((γ/ϵ) 2 log(nγ/ϵ)) is sufficient.Proof outline We note that when sampling a batch it holds w.h.p that f Bi (C i ) is close to f Xi (C i ) (Lemma 7). This is due to the fact that B i is sampled after C i is fixed. If we could show that f Bi (C i+1 ) is close f Xi (C i+1 ) then combined with the fact that we make progress of at least ϵ on the batch we can conclude that we make progress of at least some constant fraction of ϵ on the entire dataset.Unfortunately, as C i+1 depends on B i , getting the above guarantee is tricky. To overcome this issue we define the auxiliary valueC j i+1 = (1 -α j i )C j i + α j i cm(X j i ). This is the j-th center at step i + 1 if we were to use the entire dataset for the update, rather than just a batch. Note that this is only used in the analysis and not in the algorithm. Note that C i+1 only depends on C i and X and is independent of B i (i.e., we can fix its value before sampling B i ). As C i+1 does not depend on B i we use C i+1 instead of C i+1 in the above analysis outline. We show that for our choice of learning rate it holds that C i+1 , C i+1 are sufficiently close, which implies that f X (C i+1 ), f X (C i+1 ) and f Bi (C i+1 ), f Bi (C i+1 ) are also sufficiently close. That is, C i+1 acts as a proxy for C i+1 . Combining everything together we get our desired result.We start with the following useful observation, which will allow us to use Lemma 3 to bound the norm of the centers by γ throughout the execution of the algorithm.Observation 9. If ∀j ∈ [k], C j 1 is a convex combination of X then ∀i > 1, j ∈ [k], C j i , C j i are also a convex combinations of X.Let us state the following useful lemma from [14] . Although their proof is for Euclidean spaces, it goes through for Hilbert spaces. We provide the proof in the appendix for completeness. Lemma 10 ([14]). For any set S ⊆ X and any C ∈ H it holds that ∆(S, C) = ∆(S, cm(S)) + |S| ∆(C, cm(S)).We use the above to prove the following useful lemma. Lemma 11. For any S ⊆ X and C, C ′ ∈ H which are convex combinations of X, it holds that: |∆(S, C ′ ) -∆(S, C)| ≤ 2γ |S| ∥C -C ′ ∥.= |-2⟨cm(S), C ′ ⟩ + ⟨C ′ , C ′ ⟩ + 2⟨cm(S), C⟩ -⟨C, C⟩| = |2⟨cm(S), C -C ′ ⟩ + ⟨C ′ -C, C ′ + C⟩| = |⟨C -C ′ , 2cm(S) -(C ′ + C)⟩| ≤ ∥C -C ′ ∥∥2cm(S) -(C ′ + C)∥ ≤ 4γ∥C -C ′ ∥.Where in the last transition we used the Cauchy-Schwartz inequality, the triangle inequality, and the fact that C, C ′ , cm(S) are convex combinations of X and therefore their norm is bounded by γ. Now we show that due to our choice of learning rate, C j i+1 and C j i+1 are sufficiently close. Lemma 12. It holds w.h.p that ∀i ∈[t], j ∈ [k], ∥C j i+1 -C j i+1 ∥ ≤ ϵ 20γ . Proof. Note that C j i+1 -C j i+1 = α j i (cm(B j i ) -cm(X j i )). Let us fix some iteration i and center j. To simplify notation, let us denote:X ′ = X j i , B ′ = B j i , b ′ = b j i , α ′ = α j i . Although b ′ is a random variable,As what follows holds for all values of b ′ it also holds without conditioning due to the law of total probabilities.For the rest of the proof, we assume b ′ > 0 (if b ′ = 0 the claim holds trivially). Let us denote by {Y ℓ } b ′ ℓ=1 the sampled points in B ′ . Note that a randomly sampled element from X is in B ′ if and only if it is in X ′ . As batch elements are sampled uniformly at random with repetitions from X, conditioning on the fact that an element is in B ′ means that it is distributed uniformly overX ′ . Note that ∀ℓ, E[ϕ(Y ℓ )] = 1 |X ′ | x∈X ′ ϕ(x) = cm(X ′ ) and E[cm(B ′ )] = 1 b ′ b ′ ℓ=1 E[ϕ(Y ℓ )] = cm(X ′ ).Z r = r ℓ=1 (ϕ(Y ℓ ) -E[ϕ(Y ℓ )]). Note that Z 0 = 0, and when r > 0, Z r = r ℓ=1 ϕ(Y ℓ ) -r • cm(X ′ ). It is easy to see that this is a martingale:E[Z r | Z r-1 ] = E[ r ℓ=1 ϕ(Y ℓ ) -r • cm(X ′ ) | Z r-1 ] = Z r-1 + E[ϕ(Y r ) -cm(X ′ ) | Z r-1 ] = Z r-1 .Let us now bound the differences∥Z r -Z r-1 ∥ = ∥ϕ(Y r ) -cm(X ′ )∥ ≤ ∥ϕ(Y r )∥ + ∥cm(X ′ )∥ ≤ 2γ.: P r[∥Z b ′ -Z 0 ∥ ≥ δ] ≤ e Θ( δ 2 γ 2 b ′ ). Let us now divide both sides of the inequality by b ′ and set δ = b ′ ϵ 20γα ′ . We getP r[∥cm(B ′ ) -cm(X ′ )∥ ≥ ϵ 20γα ′ ] = P r[∥ 1 b ′ b ′ ℓ=1 ϕ(Y ℓ ) -cm(X ′ )∥ ≥ ϵ 20γα ′ ] ≤ e Θ( b ′ ϵ 2 (γα ′ ) 2 ) .Using the fact that α ′ = b ′ /b together with the fact that b = Ω((γ/ϵ) 2 log(nt)) (for an appropriate constant) we get that the above is O(1/ntk). Finally, taking a union bound over all t iterations and all k centers per iteration completes the proof.We can now bound the goal function when cluster centers are close using Lemma 12. Lemma 13. Fix some A ⊆ X. It holds w.h.p that ∀i ∈[t], f A (C i+1 ) -f A (C i+1 ) ≤ ϵ/5. Proof. Let S = (S ℓ ) ℓ∈[k] , S = (S ℓ ) ℓ∈[k] be the partitions induced by C i+1 , C i+1 on A. Let us expand the expression f A (C i+1 ) -f A (C i+1 ) = 1 |A| k j=1 ∆(S j , C j i+1 ) -∆(S j , C j i+1 ) ≤ 1 |A| k j=1 ∆(S j , C j i+1 ) -∆(S j , C j i+1 ) ≤ 1 |A| k j=1 4γ S j ∥C j i+1 -C j i+1 ∥ ≤ 1 |A| k j=1 S j ϵ/5 = ϵ/5.Where the first inequality is due to Observation 8, the second is due Lemma 11 and finally we use Lemma 12 together with the fact that k j=1 S j = |A|. Using the same argument we also get that f A (C i+1 ) -f A (C i+1 ) ≤ ϵ/5, which completes the proof.Let us state the following useful lemma. Lemma 14. It holds w.h.p that for every i ∈EQUATIONEQUATIONEQUATION)Proof. The first two inequalities follow from Lemma 13. The last two are due to Lemma 7 by setting δ = ϵ/5, B = B i :P r[|f Bi (C) -f X (C)| ≥ δ] ≤ 2e -bδ 2 /2γ 2 = e -Θ(bϵ 2 /γ 2 ) = e -Ω(log(nt)) = O(1/nt).Where the last inequality is due to the fact that b = Ω((γ/ϵ) 2 log(nt)) (for an appropriate constant). The above holds for either C = C i or C = C i+1 . Taking a union bound over all t iterations we get the desired result.Putting everything together We wish to lower bound f X (C i ) -f X (C i+1 ). We write the following, where the ± notation means we add and subtract a term:f X (C i ) -f X (C i+1 ) = f X (C i ) ± f Bi (C i ) -f X (C i+1 ) ≥ f Bi (C i ) -f X (C i+1 ) -ϵ/5 = f Bi (C i ) ± f Bi (C i+1 ) -f X (C i+1 ) -ϵ/5 ≥ f Bi (C i+1 ) -f X (C i+1 ) + 4ϵ/5 = f Bi (C i+1 ) ± f Bi (C i+1 ) ± f X (C i+1 ) -f X (C i+1 ) + 4ϵ/5 ≥ ϵ/5.Where the first inequality is due to inequality (3) in Lemma 14 (f X (C i ) -f Bi (C i ) ≥ -ϵ/5), the second is due to the stopping condition of the algorithm (f Bi (C i ) -f Bi (C i+1 ) > ϵ), and the last is due to the remaining inequalities in Lemma 14. The above holds w.h.p over all of the iterations of the algorithms. We conclude that when b = Ω((γ/ϵ) 2 log(γn/ϵ)), w.h.p. the algorithm terminates within t = O(γ 2 /ϵ) iterations. This complete the proof of the second claim of Theorem 1.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes]Justification: We provide code, data, and instructions on how to reproduce our experimental results in the supplemental material. Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experiments: 
 We evaluate our mini-batch algorithm on the following datasets:• MNIST: The MNIST dataset [16] has 70,000 grayscale images of handwritten digits (0 to 9), each image being 28x28 pixels. When flattened, this gives 784 features. • PenDigits: The PenDigits dataset [1] has 10992 instances, each represented by an 16-dimensional vector derived from 2D pen movements. The dataset has 10 labelled clusters, one for each digit. • Letters: The Letters dataset [28] has 20,000 instances of letters from 'A' to 'Z', each represented by 16 features. The dataset has 26 labelled clusters, one for each letter. • HAR: The HAR dataset [2] has 10,299 instances collected from smartphone sensors, capturing human activities like walking, sitting, and standing. Each instance is described by 561 features. The dataset has 6 labelled clusters, corresponding to different types of physical activities.We compare the following algorithms: (full batch) kernel k-means, mini-batch kernel k-means with the learning rate of [26] , mini-batch kernel k-means with the learning rate of sklearn, minibatch (non-kernel) k-means with the learning rate of [26] , mini-batch (non-kernel) k-means with the learning rate of sklearn. We evaluate our results with batch sizes: 1024, 256, 64, and 16. We execute every algorithm for 200 iterations. For the kernel variants we apply the Gaussian kernel: K(x, y) = e -∥x-y∥ 2 /κ , where the κ parameter is set using the heuristic of [31] followed by some manual tuning (exact values appear in the supplementary materials). We repeat every experiment 10 times and present the average Adjusted Rand Index (ARI) [11, 24] and Normalized Mutual Information (NMI) [15] scores for every dataset. All experiments were conducted on a MacBook Pro equipped with an M2 Max chip and 96 GB of RAM. We present partial results in Figure 1 and the full results in the appendix. Error bars in the plot measure the standard deviation.Figure 1 : Our results for a batch size of size 1024. We use the β prefix to denote the algorithm uses the learning rate of [26] .Proof of Lemma 10Proof. Justification: We prove our theoretical claims are true and show that the performance of our algorithm is competitive using experiments.∆(S, C) = x∈S ∆(x, C) = x∈S ⟨x -C, x -C⟩ = x∈S ⟨(x -cm(S)) + (cm(S) -C), (x -cm(S)) + (cm(S) -C)⟩ = x∈S ∆(x,Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper.• The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: Full experimental details can be found within the code.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: We explain how we report error bars using sample standard deviation.• The answer NA means that the paper does not include experiments. • The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean.• It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We give the runtime of each experiment and the type of machine they were performed on.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Justification: We have read the code of ethics and believe that our research conforms to all the requirements.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] A common variant of the above is when every x ∈ X is assigned a weight wx ∈ R + and we aim to minimize x∈X wx • minc∈C ∥c -ϕ(x)∥. Everything that follows, including our results, can be easily generalized to the weighted case. We present the unweighted case to improve readability.2 This is usually taken to be 1 -1/n p for some constant p ≥ 1. For our case, it holds that p = 1, however, this can be amplified arbitrarily by increasing the batch size by a multiplicative constant factor.Where O hides factors that are logarithmic in d, n, 1/ϵ. 
 
 Appendix:  


Decision: While the reviewers generally recommended for acceptance, none were especially enthusiastic, and–as the authors acknowledged in their response–the reviewers had consistent concerns about theoretical novelty w.r.t. Schwartzman. In light of this concern, it might be a good idea to focus on polishing the experiments to further highlight the advantages of the approach.

Reviews: review: paper_summary: paper_summary: The article presents the first mini-batch kernel k-means algorithm, which significantly improves running time compared to the full batch kernel $k$-means with only a minor negative effect on solution quality. The proposed algorithm runs in $O(n(k+b))$ time per iteration, as opposed to $O(n^2)$ for the full-batch version. The authors provide theoretical guarantees for the algorithm's performance, demonstrating that it terminates within $O(\gamma^2/\epsilon)$ iterations with high probability when the batch size is $\Omega((\gamma/\epsilon)^2 \log(n\gamma/\epsilon))$. Experimental results confirm the efficiency and effectiveness of the algorithm.
strength_weakness: strength_weakness: Approximation Quality: While the solution quality is comparable to the full-batch version, the approximation ratio depends on the batch size and initialization, which may not always guarantee optimal clustering.

Parameter Sensitivity: The performance heavily relies on parameters such as batch size and learning rate, which need careful tuning.

Complexity in Implementation: Implementing the recursive distance update and maintaining inner products can be intricate, potentially increasing the implementation complexity.

Potential Issues

Stochastic Nature: The inherent stochasticity of mini-batch algorithms can lead to variations in performance, and convergence to local minima is not guaranteed.

Parameter Initialization: Poor initialization of cluster centers can significantly affect the algorithm's performance and convergence speed.

Data Dependence: The effectiveness of the algorithm may vary depending on the dataset characteristics, such as the distribution and dimensionality of the data points.
questions: questions: See Weaknesses above
limitations: limitations: Yes
paper_summary: paper_summary: In this paper, the authors propose the first mini-batch kernel k-means clustering algorithm. It is a variant of Lloyd's algorithm that was introduced by Sculley that takes a batch of random b points instead of the full set of points and a weighted avaerage with the current centers while updating the centers. This paper attempts to translate this idea in the *kernel* k-means setting. The resulting algorithm has the same approximation guarantee as the original k-means but it terminates faster and consumes less time per iteration.

Their analysis follows the recipe of Scwartzman who used an early stopping condition when the improvement on the batch drops below some user-provided parameter. The main challenge in the kernel setting is that the underlying Hilbert space could be large or even infinite-dimensional. This is prohibitive as Scwartzman's bound on the number of iterations depends on the dimension. The authors bypass this by instead giving a bound on the Hilbert norm of the points which can be bounded in practice for example using normalized kernels.
strength_weakness: strength_weakness: New ideas are limited. Mostly an adaptation of Scwartzman's work in the kernel setting
questions: questions: None.
limitations: limitations: None.
paper_summary: paper_summary: This paper proposes the first mini-batch kernel $k$-means algorithm, which significantly reduces running time compared to the previous kernel $k$-means methods relying on the full datasets. With the proposed mini-batch kernel $k$-means algorithm, each iteration can be executed in time $O(n(k+b))$, improving the complexity of $O(n^2)$ for fully-batch methods. The authors also provide theoretical guarantees, ensuring that the algorithm can terminate (reach a convergence) within $O(\gamma^2/\epsilon)$ iterations with high probability, where $\gamma$ is the bound on the norm of points in the feature space. When initialized with the $k$-means++ seeding method, the algorithm achieves an $O(logk)$-approximation. Experimental  evaluations confirm that the mini-batch kernel k-means algorithm performs significantly faster than its full-batch counterpart while maintaining solution quality.
strength_weakness: strength_weakness: The techniques used in this paper are largely based on the work of [1]. Mini-batch $k$-means method is not new for clustering problem. The main contribution of this paper is to combine the idea of mini-batch $k$-means with the kernel $k$-means versions. It should be noted that the theoretical bounds given in this paper are not entirely novel.

There are some technical issues in the proofs (details see questions),  potentially undermining the theoretical guarantees.


[1] Gregory Schwartzman. Mini-batch $k$-means terminates within $O(d/\epsilon)$ iterations. ICLR 2024.
questions: questions: **Question 1:** For Mini-batch $k$-means methods, the centers are only updated based on a small sample of the whole datasets without calculating the exact clustering cost of the whole dataset. Thus, why the authors claim that "the execution of Lloyd's algorithm following initialization can only improve the solution and the why the approximation guarantee can remain the same for $k$-means++ method?

**Question 2:** In line 8 of Algorithm 1, I think the condition of returning the center set should be "if $f_{B_i}(C_{i}) - f_{B_i}(C_{i+1})) < \epsilon$, then return $C_{i+1}$", since $f_{B_i}(C_{i})$ is intuitivelly larger than than $f_{B_i}(C_{i+1})) $.

**Question 3:** In the proof of Lemma 7, since $a_{max} - a_{min} \le 4\gamma^2$, the probability bound by using Hoeffding Inequality should be $P_r[|f_B(C) - f_X(C)| \ge \delta] \le 2e^{-b\delta^2 / 8\gamma^4}$ instead of $2e^{-b\delta^2/2\gamma^2}$. If I am wrong, please point out.

**Question 4:** In Theorem 5, I think there is an negative sign in the probability bound, i.e., $e^{\Theta(\delta^2 / \sum_{i=1}^{m}a_i^2)}$ should be $e^{\Theta(-\delta^2 / \sum_{i=1}^{m}a_i^2)}$.

**Question 5:** In the proof of Lemma 12, I think the sample size is not enough for obtaining a success probability of around $\Omega(1 - 1/ntk)$.  $b$ should be at least $b = \Omega((\gamma^2/\epsilon)^2log(nt))$.

**Question 6:** In experiments, how to choose the kernel functions to achieve the desired performance? 

**Question 7:** The technical proofs of termination bounds given in this paper are quite similar to that for Mini-batch $k$-means proposed in ICLR 2024 (see [1] above). It seems that the only differences between these two work is the initial bound for the clustering cost difference, one is $d$ for ICLR 2024 and one is the maximum norm distance for this work. Can the authors give detailed discussions on the differences between these two work and the main technical contribution of the proposed kernel $k$-means method? 


Minor Comments:

In page 6, "The following lemma provide concentration guarantees when sampling a batch" should be "The following lemma provides concentration guarantees when sampling a batch".
limitations: limitations: Although this paper mainly gives theoretical results for clustering problems, it lacks discussions on broader impact as required by the NeurIPS guidelines.
paper_summary: paper_summary: The authors present the first mini-batch algorithm for kernel k-means. The algorithm itself is simple and works the way one would expect mini-batch kernel k-means to work. The authors improve the running time of an iteration of kernel k-means from $O(n^2)$ to $O(n(k+b))$ for the mini-batch version of the algorithm. Additionally, they show that using a specific learning rate function, there is an upper bound on the number of iterations of the algorithm.

The main challenge in the design of this algorithm is to keep track of the intermediate centers as storing the points in feature space is infeasible, as they are updated iteratively as in Lloyd's algorithm. For this, the authors design a recursive update rule to keep track of the quantity $\| \phi(x) - C_i^j \|^2$ for each iteration $i$ and each center $j$. They show that in a new iteration this quantity can be updated by considering the distance of each point in the dataset to the centers of mass of the clusters in the mini-batch and the previous centers.

Finally, the authors provide an experimental study of the mini-batch algorithm on four datasets and compare it to a non-kernel mini-batch algorithm and the full kernel k-means algorithm.
strength_weakness: strength_weakness: - The techniques, while elegant are not particularly novel in terms of theory. The proofs mostly follow from analyzing the inner product terms in the k-means formulation. 
- A number of the proofs in the main body of the paper could have been moved to the appendix, as they do not give the reader more of an understand of the big picture and are very detail specific.
- There is no discussion of the experimental results.
- While the authors state the paper is mostly theoretical, I believe this algorithm is mostly interesting to practicioners and therefore a more thorough focus on the experimental evaluation with more parameters, additional datasets and thorough discussion would strengthen the paper in my eyes.
questions: questions: Throughout the paper you assume that the computation of the kernel function $\langle \phi(x), \phi(y)\rangle$ can be done in constant time. Is there a particular reason for this? It appears to me as though many kernel functions would take time $\Theta(d)$ to compute.

In the plot for the experimental results, it is hard for me to understand the influence of the batch size on the running time, can you elaborate why we don't see bigger changes in many of the plots, and why there is such a large jump in the har dataset going from 256 to 1024?
limitations: limitations: The authors included a checklist in the appendix of the paper, but have not discussed practical limitations in the main body of the paper.

Comments: title: Read the rebuttal
comment: I have read the rebuttal from the authors. I don't have any suggestions for practical improvements. I'll maintain my score.
Thank you for your clarification and running the additional experiments.