Year: 2021

Abstract:
We devise coresets for kernel k-MEANS with a general kernel, and use them to obtain new, more efficient, algorithms. Kernel k-MEANS has superior clustering capability compared to classical k-MEANS, particularly when clusters are separable non-linearly, but it also introduces significant computational challenges. We address this computational issue by constructing a coreset, which is a reduced dataset that accurately preserves the clustering costs. Our main result is the first coreset for kernel k-MEANS with a general kernel, that has size independent of n, the number of input points; moreover, our coreset can be constructed in time near-linear in n. This result immediately implies new algorithms for kernel k-MEANS, such as a (1 + )-approximation in time near-linear in n, and a streaming algorithm using space and update time poly(k −1 log n). We validate our coreset on various datasets with different kernels. Our coreset performs consistently well, achieving small errors while using very few points. We show that our coresets can speed up kernel k-MEANS++ (the kernelized version of the widely used k-MEANS++ algorithm), and we further use this faster kernel k-MEANS++ for spectral clustering. In both applications, we achieve up to 1000x speedup while the error is comparable to baselines that do not use coresets.

heading: 1 INTRODUCTION
text: We design the first coresets, and consequently new efficient algorithms, for kernel k-MEANS and related problems, like its generalization kernel (k, z)-CLUSTERING, under general kernels. The kMEANS problem has proved to be fundamental for unsupervised learning in numerous application domains. Vanilla k-MEANS fails to capture sophisticated cluster structures, e.g., when the clusters are separable non-linearly, but this can be tackled by applying kernel methods (Schölkopf et al., 1998; Girolami, 2002). This has led to kernel k-MEANS, where data points are first mapped to a high-dimensional feature space (possibly implicitly via a kernel function), and then clustered in this richer space using a classical k-MEANS.
Formally, a kernel for a dataset X is a function K : X ×X → R+ (intended to measure similarity between elements in X) that can be realized by inner products, i.e., there exist a Hilbert space H and a map ϕ : X → H (called feature space and feature map) such that
∀x, y ∈ X, 〈ϕ(x), ϕ(y)〉 = K(x, y). (1)
In kernel k-MEANS, the input is a dataset X with weight function wX : X → R+ and a kernel function K : X × X → R+ as above, and the goal is to find a k-point center set C ⊆ H that minimizes the objective
costϕ(X,C) = ∑ x∈X wX(x) ·min c∈C ‖ϕ(x)− c‖2. (2)
(An equivalent formulation asks for a k-partitioning of X , keeping C implicit.)
This kernel version has superior clustering capability compared to classical k-MEANS (Zhang & Rudnicky, 2002; Kim et al., 2005), and has proved useful in different application domains, such as pattern recognition (Shawe-Taylor & Cristianini, 2004), natural language processing (Andrews & Fox, 2007), biology (Gönen & Margolin, 2014) and social networks (van Laarhoven & Marchiori, 2016). In fact, kernel k-MEANS is useful also for solving other clustering problems, such as normalized cut and spectral clustering (Dhillon et al., 2004; Ding et al., 2005).
Computational challenges. As observed in previous work (Girolami, 2002; Dhillon et al., 2004), the kernel trick can be applied to rewrite kernel k-MEANS using access only to the kernel K(·, ·) and without computing the very high-dimensional map ϕ explicitly. However, this approach has outstanding computational challenges (compared to classical k-MEANS), essentially because of the kernel trick. Consider the special case where k = 1 and the input is n unweighted points (i.e., 1-MEAN clustering). It is well known that the optimal center c? has a closed form c? := 1n ∑ x∈X ϕ(x). But the kernel trick requires Ω(n2) accesses to K to evaluate costϕ(X, c?),1 while in the classical setting such evaluation needs only O(n) distance computations.
This Ω(n2) barrier can be bypassed by allowing (1 + )-approximation. In particular, let S be a uniform sample of poly( −1) points from X , and let ĉ := 1|S| ∑ x∈S ϕ(x) be its 1-MEAN; then with high probability, costϕ(X, ĉ) ≤ (1 + ) costϕ(X, c?) and evaluating costϕ(X, ĉ) takes only poly( −1)n time. However, this uniform-sampling approach does not generally work for k ≥ 2, because if the optimal clustering is highly imbalanced, a uniform sample is unlikely to include any point from a small cluster. Alternative approaches, such as dimension reduction, were also proposed to obtain efficient algorithms for kernel k-MEANS, but they too do not fully resolve the computational issue. We elaborate on these approaches in Section 1.2.
Our approach. To tackle this computational challenge, we adapt the notion of a coreset (HarPeled & Mazumdar, 2004) to kernel k-MEANS. Informally, a coreset is a tiny reweighted subset of the original dataset on which the clustering cost is preserved within (1 ± ) factor for all candidate centers C ⊆ H. This notion has proved very successful for classical k-MEANS, e.g., to design efficient near-linear algorithms. In our context of kernel k-MEANS, a coreset of size s for an input of size n = |X| has a huge advantage that its k optimal center points can all be represented as linear combinations of only s points in the feature space. Given these k optimal centers (as linear combinations), evaluating the distance between a point ϕ(x) and such a center takes merely O(s2) time, instead ofO(n2), and consequently the objective can be (1+ )-approximated in timeO(s2kn). Moreover, it suffices to use k centers (again as linear combinations) that are (1 + )-approximately optimal for the coreset S.
In addition, coresets are very useful in dealing with massive datasets, since an offline construction of coresets usually generalizes to the streaming setting (Har-Peled & Mazumdar, 2004), distributed computing (Balcan et al., 2013) and dynamic algorithms (Henzinger & Kale, 2020) via the merge-and-reduce method (Har-Peled & Mazumdar, 2004), and existing (offline) algorithms can be efficiently applied to the coreset, instead of to the original dataset, with minor or no modifications.

heading: 1.1 OUR RESULTS
text: Our main result is the first coreset for kernel k-MEANS with a general kernel, that has size independent of the input size n = |X|; moreover, our coreset can be constructed in near-linear time for small k. (In fact, it generalizes to kernel (k, z)-CLUSTERING, see Section 2 for definitions.) Formally, an -coreset for kernel k-MEANS with respect to weighted dataset X and kernel function K : X × X → R+ is a weighted subset S ⊆ X , such that for every feature space H and feature map ϕ that realize K, as defined in (1),
∀C ⊆ H, |C| = k, costϕ(S,C) ∈ (1± ) · costϕ(X,C). (3)
Previously, only a weak coreset was known for kernel k-MEANS (Feldman et al., 2007), meaning that the objective is preserved only for certain candidate centers (whereas (3) guarantees this for all centers), and that coreset works only for certain kernels (finite-dimensional). While we employ a similar approach, the technical differences make our bottom-line result much stronger.
Throughout, we assume an oracle access to K takes unit time, and therefore our stated running times also bound the number of accesses to K. We denote Õ(f) = O(f · polylog f) to suppress logarithmic factors.
Theorem 1.1 (Informal version of Theorem 3.1). Given n-point weighted dataset X , oracle access to a kernel K : X ×X → R+, integer k ≥ 1 and 0 < < 1, one can construct in time Õ(nk), a
1In fact, evaluating ‖c? − ϕ(u)‖2 for a single point u ∈ X already requires Θ(n2) accesses, since ‖c? − ϕ(u)‖2 = K(u, u)− 2
n ∑ x∈X K(x, u) + 1 n2 ∑ x,y∈X K(x, y).
reweighted subset S ⊆ X of size |S| = poly(k −1), that with high probability is an -coreset for kernel k-MEANS with respect to X and K.
We can employ our coreset to devise a (1 + )-approximation algorithm for kernel k-MEANS, that runs in time that is near-linear in n and parameterized by k. This is stated in Corollary 1.2, whose proof follows by solving k-MEANS on S optimally, using straightforward enumeration over all kpartitions of S. To the best of our knowledge, such a fast (1+ )-approximation for kernel k-MEANS was not known even for k = 2; for example, uniform sampling would fail in cases where the optimal clustering is very imbalanced, as mentioned earlier. Corollary 1.2 (FPT-PTAS). Given n-point weighted dataset X , oracle access to a kernel K : X × X → R+, integer k ≥ 1 and 0 < < 1, one can compute in time O(nk + kpoly(k
−1)), a center set C of k points, each represented as a linear combination of at most poly(k −1) points from ϕ(X), such that with high probability C is a (1 + )-approximation for kernel k-MEANS on X and K. In particular, given such C, one can find for each x ∈ X its closest center in C in time poly(k −1).
In fact, for the purpose of finding near-optimal solutions, it already suffices to preserve the cost for centers coming from span(ϕ(X)) (see Fact 2.1) which is an n-dimensional subspace. However, our definition of coreset in (3) is much stronger, in that the objective is preserved even for centers coming from a possibly infinite-dimensional feature space. This stronger guarantee ensures that the coreset is composable, and thus the standard merge-and-reduce method can be applied. In particular, our coreset implies the first streaming algorithm for kernel k-MEANS. Corollary 1.3 (Streaming kernel k-MEANS). There is a streaming algorithm that given a dataset X presented as a stream of n points, and oracle access to a kernel K : X ×X → R+, constructs a reweighted subset S ⊆ X of poly(k −1) points using poly(k −1 log n) words of space and update time, such that with high probability S is an -coreset for k-MEANS with respect to X and K.
Experiments and other applications. We validate the efficiency and accuracy of our coresets on various data sets with polynomial and Gaussian kernels, which are frequently-used kernels. For every dataset, kernel, and coreset-size that we test, our coreset performs consistently better than uniform sampling which serves as a baseline. In fact, our coreset achieves less than 10% error using only about 1000 points for every dataset.
We also showcase significant speedup to several applications that can be obtained using our coresets. Specifically, we adapt the widely used k-MEANS++ (Arthur & Vassilvitskii, 2007) to the kernel setting, and we compare the running time and accuracy of this kernelized k-MEANS++ with and without coresets. On a dataset of size 105, we observe more than 1000x speedup of k-MEANS++ when using coresets, while achieving a very similar error. Furthermore, this new efficient version of kernelized k-MEANS++ (based on coresets) is applied to solve spectral clustering, using the connection discoverd by Dhillon et al. (2004). Compared to the implementation provided by Scikitlearn (Pedregosa et al., 2011), our algorithm often achieves a better result and uses significantly less time. Hence, our coreset-based approach can potentially become the leading method for solving spectral clustering in practice.

heading: 1.2 COMPARISON TO PREVIOUS APPROACHES
text: The computational issue of kernel k-MEANS is an important research topic and has attracted significant attention. In the following, we compare our result with previous work that is representative of different possible approaches for the problem.
Uniform sampling of data points is a commonly used technique, which fits well in kernel clustering, because samples can be drawn without any access to the kernel. While the coresets that we use rely on sampling, we employ importance sampling, which is non-uniform by definition. Chitta et al. (2011) employs uniform sampling for kernel k-MEANS, but instead of solving kernel k-MEANS on a sample of data points directly (as we do), their method works in iterations, similarly to Lloyd’s algorithm, that find a center set that is a linear combination of the sample. However, it has no worst-case guarantees on the error or on the running time, which could be much larger than Õ(nk). Lastly, this method does not generalize easily to other sublinear settings, such as streaming (as our Corollary 1.3). Ren & Du (2020) analyze uniform sampling for k-MEANS in a Euclidean space,
which could be the kernel’s feature space. In their analysis, the number of samples (and thus running time) crucially depends on the diameter of the dataset and on the optimal objective value, and is thus not bounded in the worse-case. This analysis builds on several earlier papers, for example Czumaj & Sohler (2007) achieve bounds of similar flavor for k-MEANS in general metric spaces.
Another common approach to speed up kernel k-MEANS is to approximate the kernel K using dimension-reduction techniques. In a seminal paper, Rahimi & Recht (2007) proposed a method that efficiently computes a low (namely, O(log n)) dimensional feature map ϕ̃ that approximates ϕ (without computing ϕ explicitly), and this ϕ̃ can be used in downstream applications. Their method is based on designing random Fourier features, and works for a family of kernels that includes the Gaussian one, but not general kernels. This method was subsequently tailored to kernel k-MEANS by Chitta et al. (2012), and another followup work by Chen & Phillips (2017) established worst-case bounds for kernel k-MEANS with Gaussian kernels. Despite these promising advances, we are not aware of any work based on dimension reduction that can handle a general kernel function (as in our approach). In a sense, these dimension-reduction techniques are “orthogonal” to our data-reduction approach, and the two techniques can possibly be combined to yield even better results.
An alternative dimension-reduction approach is low-rank approximation of the kernel matrix K. Recent work by Musco & Musco (2017) and by Wang et al. (2019) presents algorithms based on Nyström approximation to compute a low-rank (namely, O(k/ )) approximation K̃ to the kernel matrix K in time near-linear in n, such that the optimal kernel k-MEANS on K̃ achieves (1 + )- approximation to that on the original kernel K. However, the low-rank K̃ does not immediately imply efficient algorithms for kernel k-MEANS (for instance, Lloyd’s algorithm still requires Θ(n2) queries to K̃ per iteration), while our coreset can be readily combined with off-the-shelf clustering algorithms.

heading: 2 PRELIMINARIES
text: Notation. A weighted set U is a finite set U associated with a weight function wU : U → R+. For such U , let ‖U‖0 be the number of distinct elements in it. For a weight function as above and a subset S ⊆ U , define wU (S) := ∑ u∈S wU (u). For any other map f : U → V (not a weight function), we follow the standard definition f(S) := {f(x) : x ∈ S}. For an integer t ≥ 1, let [t] := {1, . . . , t}. For a real number x and an integer i ≥ 1, let log(i) x be the i-th iterated log of x, i.e., log(1) x = log x and for i ≥ 2 let log(i) x = log(log(i−1) x).
Kernel functions. Let X be a set of n points. A function K : X ×X → R+ is a kernel function if the n × n matrix M such that Mij = K(xi, xj) (where xi, xj ∈ X) is positive semi-definite. Since M is positive semi-definite, there exists a map ϕ from X to some Hilbert space H, such that all x, y ∈ X satisfy K(x, y) = 〈ϕ(x), ϕ(y)〉. This above (existence of a map ϕ of into H) can be extended to infinite X , e.g., X = Rd, by Mercer’s Theorem. The distance between x′, y′ ∈ H is defined as dist(x′, y′) := ‖x′ − y′‖ = √ 〈x′ − y′, x′ − y′〉. Hence, the distance dist(ϕ(x), ϕ(y)) for x, y ∈ X can be represented using K as
dist(ϕ(x), ϕ(y)) = ‖ϕ(x)− ϕ(y)‖ = √ K(x, x) +K(y, y)− 2K(x, y).
We refer to a survey by Ghojogh et al. (2021) for a more comprehensive introduction to kernel functions.
Kernel (k, z)-CLUSTERING. In the the kernel (k, z)-CLUSTERING problem, the input is a weighted data set X of n objects, a kernel function K : X × X → R+, an integer k ≥ 1, and z > 0. The goal is to find a k-point center set C ⊆ H that minimizes the objective
costϕz (X,C) := ∑ x∈X wX(x)(dist(ϕ(x), C)) z, (4)
whereH is an induced Hilbert space of K and ϕ : X → H is its feature map, and dist(ϕ(x), C) := minc∈H dist(ϕ(x), c) = minc∈H ‖ϕ(x)− c‖. The case z = 2 clearly coincides with kernel kMEANS whose objective is (2). The (non-kernel) (k, z)-CLUSTERING problem may be viewed as kernel (k, z)-CLUSTERING with kernel K(x, y) = 〈x, y〉 and identity feature map ϕ(x) = x.
While the feature map ϕ might not be unique, we show below that this kernel (k, z)-CLUSTERING is well defined, in the sense that the optimal value is independent of ϕ. The following two facts are standard and easy to prove. Fact 2.1. For every map ϕ into H, there is an optimal solution C? in which every center point c ∈ C? lies inside span(ϕ(X)), and is thus a linear combination of ϕ(X). Corollary 2.2. The optimal value of (4) can be represented as a function of kernel values K(x, y), and is thus invariant of ϕ.
-Coresets for kernel (k, z)-CLUSTERING. For 0 < < 1, an -coreset for kernel (k, z)CLUSTERING on a weighted dataset X and a kernel function K is a reweighted subset S ⊆ X , such that for every Hilbert spaceH and map ϕ : X → H satisfying (1), we have
∀C ⊆ H, |C| = k, costϕz (S,C) ∈ (1± ) · costϕz (X,C). The case z = 2 clearly coincides with (3).
3 CORESETS FOR KERNEL (k, z)-CLUSTERING
Theorem 3.1. Given n-point weighted dataset X , oracle access to a kernel K : X × X → R+, z ≥ 1, integer k ≥ 1, and 0 < < 1, one can construct in time Õ(nk), a reweighted subset S ⊆ X of size ‖S‖0 = 2O(z) · poly(k −1), that with high constant probability is an -coreset for kernel (k, z)-CLUSTERING with respect to X and K.
At a high level, our prove this theorem by employing recent constructions of coresets for (k, z)CLUSTERING in Euclidean spaces, in which the coreset size is independent of the Euclidean dimension (Sohler & Woodruff, 2018; Feldman et al., 2020; Huang & Vishnoi, 2020; Braverman et al., 2021). However, these coresets are designed for finite-dimensional Euclidean spaces, and are thus not directly applicable to our feature spaceH, which might have infinite dimension. To employ these coreset constructions, we show in Lemma 3.2 that the data points in the feature space H embed into an (n + 1)-dimensional (Euclidean) space, without any distortion to distances between data points and centers. This observation is similar to one previously made by Sohler & Woodruff (2018) for a different purpose. Due to this embedding, it suffices to construct a coreset for the limited setting where centers come only from an (n+ 1)-dimensional space (Corollary 3.3). Lemma 3.2. Let H be a Hilbert space and let X ⊆ H be a subset of n points. Then there exists a map f : H → Rn+1 such that
∀x ∈ X, c ∈ H, ‖x− c‖ = ‖f(x)− f(c)‖.
Proof. Let S = span(X). Then every point c ∈ H can be written (uniquely) as c = c‖+ c⊥, where c‖ ∈ S and c⊥ is orthogonal to S. Thus, ‖c‖2 = ‖c‖‖2 + ‖c⊥‖2. Note that for all x ∈ X , we have x⊥ = 0. Now, for every c ∈ H, let f(c) := (c‖ ; ‖c⊥‖), where we interpret x‖ as an n-dimensional vector. Then for all x ∈ X and c ∈ H, ‖x− c‖2 = ‖x‖ − c‖‖2 + ‖x⊥ − c⊥‖2 = ‖x‖ − c‖‖2 + ‖c⊥‖2 = ‖f(x)− f(c)‖2. The claim follows.
Corollary 3.3. Consider n-point weighted dataset X , kernel function K : X ×X → R+, z ≥ 1, integer k ≥ 1, and 0 < < 1. Suppose that a reweighted subset S ⊆ X satisfies that for every ϕ : X → Rn+1 such that for all x, y ∈ X , 〈ϕ(x), ϕ(y)〉 = K(x, y), the following holds
∀C ⊆ Rn+1, |C| = k, costϕz (S,C) ∈ (1± ) · costϕz (X,C). Then S is an -coreset for kernel (k, z)-CLUSTERING with respect to X and kernel K.
Proof. To verify that S is a coreset with respect to X and K, consider some feature space H′ and feature map ϕ′ be induced by K. Apply Lemma 3.2 to obtain f : H′ → Rn+1, then for all C ⊆ H′, |C| = k, we have ∀Q ⊆ X, costϕ′z (Q,C) = costf◦ϕ ′
z (Q, f(C)), and using the promise about S with ϕ = f ◦ ϕ′,
costϕ ′
z (S,C) = cost ϕ z (S, f(C)) ∈ (1± ) · costϕ(X, f(C)) = (1± ) · costϕ
′ (X,C).
Thus, S is indeed a coreset with respect to X and K.
Algorithm 1 Constructing -coreset for kernel (k, z)-CLUSTERING on data set X with kernel K 1: let X0 ← X , i← 0 2: repeat 3: let i← i+ 1 and i ← /(log(i) ‖X‖0)1/4 4: Xi ← IMPORTANCE-SAMPLING(Xi−1, i) 5: until ‖Xi‖0 does not decrease compared with ‖Xi−1‖0 6: return Xi
Another issue is that some of the existing algorithms, such as Sohler & Woodruff (2018); Feldman et al. (2020), require the input to be an explicit representations of points in ϕ(X), which is very expensive to compute. Fortunately, the importane-sampling-based algorithms of Huang & Vishnoi (2020) and Braverman et al. (2021) are oblivious to the representation of ϕ, and only rely on a distance oracle that evaluates ‖ϕ(x) − ϕ(y)‖ = √ K(x, x) +K(y, y)− 2K(x, y) for data points x, y ∈ X . Now, by Corollary 3.3, executing these algorithms without any modifications (except for plugging in the distance oracle defined by kernel K) yields the desired coreset for kernel (k, z)CLUSTERING. We choose to use the coreset construction of Braverman et al. (2021), which is arguably simpler. We now recall its statement for completeness.
Theorem 3.4 (Braverman et al. (2021)). Given n-point weighted dataset X ⊂ Rm for some integer m, together with z ≥ 1, integer k ≥ 1, and 0 < < 1, one can construct in time Õ(nk), a reweighted subset S ⊆ X of size ‖S‖0 = Õ( −422zk2), that with high constant probability is an -coreset for (k, z)-CLUSTERING with respect to X .
Proof of Theorem 3.1. It follows immediately by combining Corollary 3.3 and Theorem 3.4.

heading: 3.1 DESCRIPTION OF CORESET CONSTRUCTION ALGORITHMS
text: Next, we present our full algorithm in Algorithm 1 (which depends on the subroutines defined in Algorithm 2 and 3). While it essentially tailors previous work to our kernel clustering setting, we provide full details for completeness. The following notation is needed. For a subset C ⊆ H and data point x ∈ X , define NNC(x) := arg min{dist(ϕ(x), y) : y ∈ C} as the nearest neighbor of x in C with respect to the distances in the feature space (breaking ties arbitrarily). Thus NNC(·) defines a |C|-partition ofX , and the cluster that x belongs to (with respect to C) is denoted C(x) := {x′ ∈ X : NNC(x′) = NNC(x)}. Algorithm 1 is the main procedure for constructing the coreset, and its loop iteratively executes another importance-sampling-based coreset construction (Algorithm 2). Informally, each invocation of IMPORTANCE-SAMPLING constructs a coreset Xi from the current coreset Xi−1, to reduce the number of distinct elements in Xi−1 to roughly log ‖Xi−1‖0. The procedure ends when such size reduction cannot be done any more, at which point the size of the coreset reaches the bound in Theorem 3.4, which is independent of n.
In fact, subroutine IMPORTANCE-SAMPLING already constructs a coreset, albeit it is of a worse size that depends on log ‖X‖0. This subroutine is based on the well-known importance sampling approach that was proposed and improved in a series of works (cf. Langberg & Schulman (2010); Feldman & Langberg (2011); Feldman et al. (2020)). Its first step is to compute an importance score σx for every data point x ∈ X (lines 1–2), and then draw independent samples from X with probability proportional to σx (lines 3–4). The final coreset is formed by reweighting the sampled points (lines 5–6). Roughly speaking, the importance score σx measures the relative contribution of x to the objective function in the worst-case, which here means the maximum over all choices of the center set. It can be computed from an O(log k)-approximate solution for kernel (k, z)CLUSTERING on X , which is generated by the Dz-sampling subroutine (Algorithm 3), a natural generalization of the D2-sampling introduced by Arthur & Vassilvitskii (2007).
We stress that our algorithm description uses the feature vectors ϕ(x) for clarity of exposition. These vectors do not have to be provided explicitly, because only distances between them are required, and thus all steps can be easily implemented using the kernel trick, and the total time (and number of accesses to the kernel function K) is only Õ(nk).
Algorithm 2 IMPORTANCE-SAMPLING(X , ) 1: let C? ← Dz -SAMPLING(X) 2: for each x ∈ X , let σx ← wX(x) · ( (dist(x,C?))z
costϕz (X,C?) + 1wX(C?(x)) ) 3: for each x ∈ X , let px ← σx∑
y∈X σy
4: draw N = O( −422zzk2 log2 k log ‖X‖0) i.i.d. samples from X , using probabilities (px)x∈X 5: let D be the sampled set, and for each x ∈ D let wD(x)← wX(x)pxN 6: return weighted set D
Algorithm 3 Dz -SAMPLING(X) . the feature vectors ϕ(·) are mentioned for clarity and not needed for implementation
1: let x be a uniform random point from X , and initialize C ← {ϕ(x)} 2: for i = 1, . . . , k − 1 do 3: draw one sample x ∈ X , using probabilities wX(x) · (dist(ϕ(x),C)) z
costϕz (X,C)
4: let C ← C ∪ {ϕ(x)} 5: end for 6: return C

heading: 4 EXPERIMENTS
text: We validate the empirical performance of our coreset for kernel k-MEANS against various datasets, and show that our coresets can significantly speedup a kernelized version of the widely used k-MEANS++ (Arthur & Vassilvitskii, 2007). In addition, we apply this new coreset-based kernelized k-MEANS++ to spectral clustering (via a reduction devised by Dhillon et al. (2004)), showing that it outperforms the well-known Scikit-learn solver in both running time and objective value.
Experimental setup. Our experiments are conducted on standard clustering datasets that consist of vectors in Rd, and we use the RBF kernel (radial basis function kernel, also known as Gaussian kernel) and polynomial kernels as kernel functions. An RBF kernel KG is of the form KG(x, y) := exp ( −‖x−y‖2σ2 ) , where σ > 0 is a parameter, and a polynomial kernel KP is of the formKP (x, y) := (〈x, y〉+c)d where c and d are parameters. Table 1 summarizes the specifications of datasets and our choice of the parameters for the kernel function. We note that the parameters are dataset-dependent, and that for Twitter and Census1990 dataset we subsample to 105 points since otherwise it takes too long to run for some of our inefficient baselines. Unless otherwise specified, we use a typical value k = 5 for the number of clusters. All experiments are conducted on a PC with Intel Core i7 CPU and 16 GB memory, and algorithms are implemented using C++.

heading: 4.1 SIZE AND EMPIRICAL ERROR TRADEOFF
text: Our first experiment evaluates the empirical error versus coreset size. In our coreset implementation, we simplify the construction in Algorithm 1 by running the importance sampling step only once instead of running it iteratively, and it turns out this simplification still achieves excellent performance. As in many previous implementations, instead of setting and solving for the number of samples N in the IMPORTANCE-SAMPLING procedure (Algorithm 2), we simply set N as a parameter to directly control the coreset size. We construct the coreset with thisN and evaluate its error by drawing 500 random center sets C (each consisting of k points) from the data set, and evaluate the maximum empirical error, defined as
̂ := max C∈C | costϕ(S,C)− costϕ(X,C)| costϕ(X,C) . (5)
This empirical error is measured similarly to the definition of coreset, except that it is performed on a sample of center sets. To make the measurement stable, the empirical error is evaluated independently 100 times and the average is reported.
The tradeoff of the coreset size versus empirical error is shown in Figure 1, where we also compare with a baseline that constructs the coreset using uniform sampling. These experiments show that our coreset performs consistently well on all the datasets and kernels. Furthermore, our coreset admits a similar error curve regardless of dataset and kernel function – for example, one can get within 10% error using a coreset of only 1000 points, – which is perfectly justified by our theory that the size of -coreset only depends on and k. Comparing with the uniform-sampling baseline, our coreset generally has superior performance, especially when the coreset size is small. We also observe that the uniform sampling suffers a larger variance compared with our coreset.
4.2 SPEEDING UP KERNELIZED k-MEANS++
k-MEANS++ (Arthur & Vassilvitskii, 2007) is a widely-used algorithm for k-MEANS, and it could easily be adopted to solve kernel k-MEANS by using the kernel trick, however as mentioned earlier this would take Ω(n2) time. We use our coreset to speedup this kernelized k-MEANS++ algorithm, by first computing the coreset and then running kernelized k-MEANS++ on the coreest; this yields an implementation of kernelized k-MEANS++ whose running time is near-linear (in n).
In Figure 2, we demonstrate the running time and the error achieved by kernelized k-MEANS++ with and without coresets, experimented with varying coreset sizes. We measure the relative error of our coreset-based kernelized k-MEANS++ by comparing the objective value it achieves with that of vanilla (i.e., without coreset) kernelized k-MEANS++. These experiments show that the error decreases significantly as the coreset size increases, and it stabilizes around size N = 100, achieving merely < 5% error. Naturally, the running time of our coreset-based approach increases with the coreset size, but even the slowest one is still several orders of magnitude faster than vanilla kernelized k-MEANS++.

heading: 4.3 SPEEDING UP SPECTRAL CLUSTERING
text: In the spectral clustering problem, the input is a set of n objects X and an n × n similarity matrix A that measures the similarity between every pair of elements in X , and the goal is to find a k-
partition of X such that a certain objective function with respect to A is minimized. Dhillon et al. (2004) shows a way to write spectral clustering as a (weighted) kernel k-MEANS problem, which is eventually used to produce a spectral clustering. Specifically, let D be an n × n diagonal matrix such that Dii = ∑ j∈[n]Aij . Then according to Dhillon et al. (2004), spectral clustering can be written as a weighted kernel k-MEANS problem with weights wi := Dii and kernel function K := D−1AD−1, provided that A is positive semidefinite (which could be viewed as a kernel). We use this reduction, and plug in the abovementioned coreset-based kernelized k-MEANS++ as the solver for kernel k-MEANS. We experiment on the subsampled Twitter dataset with varying number of points, and we use the polynomial and RBF kernels as the similarity matrix A.
However, we would need Θ(n2) time if we evaluate Dii naively. To resolve this issue, we draw a uniform sample S from [n], and use D̂ii := n|S| ∑ j∈S Aij as an estimate for Dii. The accuracy of
D̂ is justified by previous work on kernel density estimation (Joshi et al., 2011, Theorem 5.2), and for our application we simply set |S| = 1000, which achieves good accuracy. We compare our implementation with the spectral clustering solver from the well-known Scikitlearn library (Pedregosa et al., 2011) as a baseline. The experimental results, reported in Figure 3, show that our approach has more than 1000x of speedup already for moderately large datasets (n = 40000). This difference might be partially caused by efficiency issues of the Python language used for the Scikit-learn implementation (recall that our implementation is in C++), but we actually see that the asymptotic growth of our algorithm’s running time is also much better than that of the Scikitlearn baseline. This suggests that our improvement in running time is fundamental, and not only due to the programming language. We also observe that our approach yields better objective values than Scikit-learn. One possible reason is that Scikit-learn might be conservative in using more iterations to gain better accuracy, because of the expensive computational cost that we do not suffer.

Review summaries:
The paper claims to present the first coreset for kernel k-means and other (k,z) kernel clusterings given oracle access to the kernel. The size of coreset is independent of the number of points and construction time is near linear in k. The authors use the coreset to get a (1+\epsilon) approximation for the kernel k-means. The theory is supported with experimental results.