Beyond the Nystrom Approximation: Speeding up Spectral Clustering using Uniform Sampling and Weighted Kernel k-means

1 Introduction

Clustering is one of the fundamental problems in machine learning. Spectral clustering techniques [Von Luxburg, 2007] have been shown to outperform other clustering methods, such as k-means or single-linkage algorithms. However, two significant obstacles to scaling up spectral clustering to large datasets include building an affinity matrix between pairs of data points and computing the eigenvectors of the resulting normalized graph Laplacian, both of which become computationally prohibitive for large data-sets.

Approaches to scaling up spectral clustering typically focus on improving the efficiency of computing the eigenvectors using various matrix approximation schemes. A popular technique for matrix approximation is the Nystr¨om method [Williams and Seeger, 2001] which constructs a low rank approximation by sampling a subset of the columns of the affinity matrix. Various sampling schemes have been proposed in the literature to sample this subset of columns [Kumar et al., 2012; Gittens and Mahoney, 2013], offering a tradeoff between the complexity of the sampling distribution and quality of the approximation. However other than uniform sampling, the practicality of these approaches for massive datasets may be limited [Yan et al., 2009]. The major disadvantage of these methods is that if l columns are sampled, they involve the computation of eigenvectors of a sub-matrix of size l × l, which can itself be computationally expensive when l is large.

Also, the subsequent k-means step used for rounding the cluster assignment has to be performed for all the input points.

The connection between spectral clustering and weighted kernel k-means was first introduced in [Dhillon et al., 2004]. In this paper, we use this connection to propose a simple framework for spectral clustering that samples a subset S of the input points, computes the clusters for S using weighted kernel k-means and uses the resulting centers to compute a clustering for rest of the data points. The proposed framework does not require computation of the entire affinity matrix and does not require the computation of eigenvectors of the selected sub-matrix. Sampling algorithms developed for the Nystr¨om approximation can be leveraged to obtain a tradeoff between accuracy and computational efficiency. In fact, our experiments show that the use of sampling followed by weighted kernel k-means outperforms sampling followed by the Nystr¨om approximation. Since the proposed method uses weighted kernel k-means to compute the clusters, it is easy to parallelize. The contributions of this paper are:

• Prior work on the Nystr¨om approximation such as [Gittens and Mahoney, 2013; Talwalkar and Rostamizadeh, 2010] provide approximation guarantees with respect to various matrix norms, but not with respect to the spectral clustering objective itself. We show that matrix norm errors are not always indicative of the quality of spectral clustering.
• We propose a framework for spectral clustering that applies sampling followed by the weighted kernel k-means algorithm instead of approximating the affinity matrix or the Laplacian.
• Without making any assumptions about the data, we show that when points are sampled uniformly at random without replacement, the following theorem holds.

Theorem 1. Let 0 < δ < 1, α ≥ 1, 0 < β < 1 and ϵ > 0 be approximation parameters. Let A be an α-approximation algorithm for the weighted kernel k-means problem. Given a set of n points, V , suppose we sample a subset S ⊂ V of size s uniformly at random without replacement such that,

where ∆ = maxi,j wiwj‖φ(ai) - φ(aj)‖2 . If we run algorithm A with input S, then for the solution C ∗ obtained, with probability at least 1- δ,

In other words, by choosing sufficient number of samples s, we can obtain centers that provide a good approximation to the optimal cluster centers. As the input size n increases, the term (1/n) vanishes, resulting in a bound for s that is independent of n. To the best of our knowledge, the relation between the number of samples and the spectral clustering objective studied in the proposed work has not been explored in the literature.

2 Preliminaries

2.1 Weighted Kernel k-means

Given a set of n points V = {v1, v2..., vn} with associated weights {w1, w2...wn} and a kernel matrix K, the weighted kernel k-means objective for clusters V1, V2, ...Vk with centers C = {C1, C2, ..Ck} is defined as:

where Ci = ∑vj∈Vi wjφ(vj) ∑vj∈Vi ‖wj‖ . Here φ(v) is the kernel function that maps the point v to a higher dimensional feature space. It was shown in [Dhillon et al., 2004],

Definition 1. Let A be an approximation algorithm for the weighted kernel k-means objective. Let α ≥ 1. A is an αapproximation algorithm if the set of centers, C, returned by A satisfies,

where Copt is the set of optimal centers.

2.2 Spectral Clustering Using Norm-Cuts

In spectral clustering, we are given a graph G = (V,A), which is made up of a set of n vertices V . The affinity matrix A is n × n whose entries represent the similarity between vertices. If V1, V2 are subsets of V , let links(V1, V2) = ∑i∈V1,j∈V2 Aij .

Furthermore, let degree(V1) = links(V1, V ). The graph partitioning problem seeks to partition the graph into k disjoint clusters V1, ..., Vk. A number of different graph partitioning objectives have been proposed and studied. In this paper, we will focus on the normalized cut objective. The goal is to minimize the following objective over all possible clusterings {V1..Vk},

2.3 Relation BetweenWeighted Kernel k-means And Spectral Clustering

To convert a spectral clustering problem to a weighted kernel k-means problem, it was shown in [Dhillon et al., 2004] that we can set W = D and K = σD-1 + D-1AD-1 , where σ = k/(n - k). Here the term σD-1 is added to ensure K is positive definite and does not change the optimal clustering.

We note that a set of centers C for the weighted kernel kmeans problem induces a clustering V1...Vk for the spectral clustering problem. Since both objectives are equivalent, we use W(V,C), NCut(G, {V1..Vk}) and NCut(G,C) interchangeably throughout the remainder of this paper. We also use NCut(G, k) = minC NCut(G, C) to represent the optimal partitioning of G into k clusters.

3 WhyMatrix Norms Are Insufficient For Norm Cuts-based Spectral Clustering

The Nystr¨om approximation has been extensively studied with respect to various matrix error norms, such as the Frobenius norm, trace norm and the spectral norm [Gittens and Mahoney, 2013; Talwalkar and Rostamizadeh, 2010]. In this section, we show that approximation of these norms is neither necessary nor sufficient for producing a good clustering with respect to the norm cuts objective. For simplicity, we focus on the trace norm error. However the examples can be extended to any matrix norm. Also, for this section, we let NCut(A, k) represent the optimal partitioning of the graph with affinity matrix A into k clusters.

Lemma 1. Preserving trace norm is not sufficient for an approximation to be good with respect to the norm cuts objective.

We construct two matrices that are extremely similar in terms of their trace norms, but have a significantly different norm cuts objective. Here, the trace norm of a matrix A is denoted as ‖A‖∗ = ∑ni=1 |λi(A)|, where λi(A) is the ith eigenvalue of A.

In other words, we show that for any value of ϵ > 0, there exists n, k, A and B such that, A and B represent graphs GA, GB with n vertices, and ‖A - B‖∗ ≤ (1 + ϵ)‖A‖∗ and NCut(A, k) = O(k)·NCut(B, k). Note, that the maximum value of NCut(A, k) is k.

Let A and B be the affinity matrices for undirected, unweighted graphs that have n vertices and k and k - 1 equal sized components, respectively. Further, we assume that each component is completely connected with no self loops. Since each component in A is a complete graph on n/k vertices, its eigenvalues are -1 with multiplicity (n/k)-1 and (n/k)-1 with multiplicity 1. Since each component in A is regular, the eigen values of the normalized symmetric Laplacian of A (denoted by L(A)) are given as, λi(L(A)) = 1 - (1/∆)λi(A) where ∆ = (n/k) - 1. In other words, L(A) has eigenvalues n/(n - k) with multiplicity n - k and 0 with multiplicity k. Similarly, L(B) has eigenvalues n/(n - k + 1) with multiplicity n-k+1 and 0 with multiplicity k-1. A more general statement in this context can be found in [Butler, 2007].

Now we compare the trace norm of the difference of the

Algorithm 1 Proposed Framework for Spectral Clustering

Input: Affinity Matrix A, number of clusters k, number of samples s

Output: Matrix ˆ Y ( ˆ Yij ̸= 0 only if point i belongs to cluster j)

Procedure:

sub ← sample(s)

Asub ← A(sub, sub)

Y ←weighted kernel kmeans(Asub, k)

ˆ Y ← diffuse(Y, A, sub)

two Laplacians, L(A) and L(B), as follows,

where, k′ = k - 1. Given a value of ϵ, we can choose appropriate values of n, k such that ‖L(A)‖2∗ - ‖L(B)‖2∗ ≤ (1 + ϵ)‖L(A)‖2∗ . Now we examine the norm cuts ratio when we try to partition these graphs into exactly k clusters. Thus, NCut(A, k) = 0, since the graph has exactly k components. However for B, one of the k - 1 components will have to be split into two equal parts to minimize the norm cuts ratio. This results in a cut that involves n k-1 ( n k-1 +1)/2 edges. Thus the norm cuts ratio for B is given as,

The second term reduces to a constant for a sufficiently large value of n, resulting in a norm cuts ratio of O(k).

Lemma 2. Preserving trace norm is not necessary for an approximation to be good with respect to the norm cuts objective.

Proof. To show that preserving matrix norms is not necessary for an approximation to be good with respect to the norm cuts objective, we consider the case of a block diagonal matrix A with k blocks. Let L(A) be the corresponding normalized Laplacian. The result of norm cuts based spectral clustering will depend on the eigenvectors of L(A) and the relative order of the corresponding eigenvalues. We note that L(A) has the same eigenvectors as L(A)2 , L(A)3 , etc. Thus, L(A), L(A)2 return the same spectral clustering. However ‖L(A)‖∗ - ‖L(A)i ‖∗ can be arbitrarily high.

4 Proposed Framework

The proposed framework for spectral clustering is described in Algorithm 1. The overall procedure consists of

Algorithm 2 Diffuse

Input: Matrix ˆ Y , Affinity Matrix A, set of sampled indices

SOutput: Matrix Y (Yij ̸= 0 only if point i belongs to cluster j)

Procedure:

end if end for

end for

for each i /∈ S do

for c = 1 : k do

three stages. In the first stage, a subset of the input points are sampled according to some distribution. This subset is indicated as sub. In the second stage, the clusters are computed using weighted kernel k-means. However, instead of using the entire affinity matrix, as in [Dhillon et al., 2004], we only use the principal sub-matrix induced by indices in sub. The final step involves using the cluster centers returned by the weighted kernel k-means procedure to assign clusters to the remaining data points. The procedure for this step is described in Algorithm 2.

This is scalable to larger datasets because it does not require explicit computation of the eigenvectors. Also it only requires the computation of a sub-matrix of size n×k, which can result in significant savings in memory. Weighted kernel k-means is extremely easy to parallelize.

5 Analysis

In this section, we present the proof of Theorem 1. We extend the analysis in [Czumaj and Sohler, 2010] to the case of sampling without replacement. Specifically, uniform sampling without replacement allows us to use the sharper bound offered by the Serfling inequality stated below. We refer the reader to [Bardenet et al., 2015] for a more thorough discussion.

Lemma 3. [Serfling, 1974; Bardenet et al., 2015] LetX1..Xs be random variables sampled uniformly without replacement.

Let X = X1 +...+Xn and 0 ≤ Xi ≤ ∆. Then for all γ > 0,

Definition 2. Let β > 0 and α ≥ 1. A set of k centers C is a β-bad, α-approximation if

If not, C is said to be a β-good, α-approximation.

We will need the following lemmas.

Lemma 4. Let S be a subset of V of size s such that,

If an α-approximation algorithm for weighted kernel kmeans A is run on input S, then for the set of centers, CA, obtained, with probability ≥ 1 - δ,

where Copt is the set of optimal centers for V and m = W(V, Copt)/n.

Proof. For each point vi ∈ S, let Xi = wi‖φ(vi) - Ci‖2 where Ci the nearest center to φ(vi) in Copt. Then W(S,Copt) = s∑ i=1 Xi. We also know that E[Xi] = m. Let

Xi be bounded in the interval [0,∆].

The last step is obtained by applying the Serfling inequal-

Let C ′ be the set of k centers in S obtained by replacing each C ∈ Copt by its nearest neighbor in S. From the weakened triangle inequality, it can be shown that W(S,C ′ ) ≤ 4W(S,Copt). This implies with probability at least 1 - δ, S contains k centers with cost at most 4(1 + βα ) sn W(V,Copt). In other words,

If C ′′ is the set of optimum centers for W(S,C), then

Since A is an α-approximation algorithm,

This implies,

The lemma now follows.

Lemma 5. Let S be a subset of V of size s, sampled uniformly at random such that,

Let χ be the set of 12β-bad 4α-approximations of W(V, Copt) then,

Proof. Let s ≥ 4α+5β β k. Let Cb ∈ χ be a 12β-bad 4αapproximation of W(V, Copt). If Cb ⊂ S, let S∗ be the subset of S with Cb removed, such that |S∗ | = s - k. Thus,

The second step holds because the mean remains unchanged at each step when the elements are chosen uniformly at random without replacement and W(S, Cb) = W(S∗ , Cb) when Cb ⊂ S. The third step holds since s ≥ 4α+5β β k.

Let Yi denote the distance of point vi ∈ S from its closest center in Cb. Since Cb ∈ χ, E[Yi] ≥ (4α + 12β)m. Thus,

Similar to the previous lemma, we have used the Serfling bound to obtain (2). Thus it follows that,

We note, Pr[Cb ⊂ S] ≤ (s/n)k and |χ| ≤ nk . Using this and (3) in (1), we have,

which is upper bounded by δ when,

Next, we use these lemmas to prove Theorem 1.

Proof. Let β∗ be some constant. Let S ⊆ V be a subset of s points chosen uniformly at random such that,

Then from Lemma 3, we have, if C ⊆ S and C ∈ χ, with probability at least 1- δ,

On the other hand, if we run algorithmA for set S, then by Lemma 2, the resulting set of centers satisfies, with probability at least 1- δ,

This implies with probability at least 1 - 2δ the set C is a (12β)-good (4α)-approximation. In order to remove the dependence on the number of samples on m,

• If m ≤ ϵ, we simply set β = (1/√2)ϵ/m. In this case we will obtain the following bound on s.

• If m > ϵ, we set β = 3β∗ to obtain,

The theorem now follows.

6 Experiments

We now present experimental results to show that the use of sampling followed by weighted kernel k-means outperforms sampling followed by the Nystr¨om approximation. We evaluated the performance of two variants of the proposed framework. The first variant uses uniform sampling without replacement to sample a subset of the input points (denoted as wkk-unif). The second variant samples points based on the norm of the corresponding column in the affinity matrix (denoted as wkk-col). These results were compared against approaches described in [Dhillon et al., 2004] (denoted as kulis) and [Yu and Shi, 2003] (denoted as shi). Both of these methods use the entire affinity matrix and do not perform any sampling. We also compare against the Nystr¨om approximation when the columns are sampled uniformly at random (denoted as nys-unif) and sampled based on their column norm (denoted as nys-col).

6.1 Synthetic Data

We evaluated the norm cuts ratio and the computation time on six commonly used synthetic datasets [Bouneffouf and Birol, 2015], described in Table 1, and repeated our evaluations 10 times. We measured the clustering quality of each algorithm using the average accuracy across different datasets.

The results are shown in Figure 1. The first row shows the norm cuts ratio obtained by various approaches. It can be seen that norm cuts ratio for wkk-unif is comparable to kulis. This supports the result presented in Theorem 1. Furthermore, the second row in Figure 1 shows the time taken by each approach. This highlights the speedup obtained by using wkkunif and wkk-col over shi and kulis. For all the datasets, the approaches based on combining sampling and the weighted kernel k-means algorithm outperform approaches based on the Nystr¨om approximation both in terms of accuracy and computation time.

6.2 Image Segmentation

One of the most popular applications of spectral clustering is image segmentation. In this section, we describe results obtained on an image segmentation benchmark [Yu and Shi, 2003]. The affinity matrix and the final discretization were computed using the approach of [Yu and Shi, 2003].

Experimental results in Figure 2 show that wkk-unif and wkk-col obtain the same accuracy as the Nystr¨om based approaches, while taking significantly less time. In cases where only 10% of the columns were sampled, wkk-unif and wkkcol achieved a lower norm cuts ratio compared to nys-unif and nys-col. The baseline methods shi and kulis use the entire affinity matrix. Thus they can achieve lower norm cuts ratio, but have a significantly higher computation time.

In order to study how the proposed approach scales to larger graphs, we vary the size of the input images. The results are shown in Figure 3. For shi, the time taken scales cubically with the size of the graph. Hence it is unsuitable for large graphs. In contrast, kulis, wkk-unif and wkk-col scale almost linearly with the size of the graph. Both wkk-unif and wkk-col outperform kulis in terms of execution time. They also outperform the Nystr¨om-based approaches both in terms of accuracy and computation time. These results show that uniform sampling followed by the weighted kernel k-means algorithm scales well in practice.

7 Conclusion

This paper presented a framework for spectral clustering based on combining sampling and the weighted kernel kmeans algorithm. If the points are sampled uniformly at randomwithout replacement, we show that for large datasets, the number of samples required, is independent of the input size and depends only on the number of clusters and the diameter of the set of points in the kernel space. Experiments show that approaches based on the proposed framework outperform approaches based on the Nystr¨om approximation.