Spectral Clustering and Cheeger’s Inequality

3 Spectral Clustering and Cheeger’s Inequality 3.1 Clustering Clustering is one of the central tasks in machine learning. Given a set of data points, the purpose of clustering is to partition the data into a set of clusters where data points assigned to the same cluster correspond to similar data points (depending on the context, it could be for example having small distance to each other if the points are in Euclidean space). 3.1.1 k-means Clustering One the most popular methods used for clustering is k-means clustering. Given x1 , . . . , xn ∈ Rp the k-means clustering partitions the data points in clusters S1 ∪ · · · ∪ Sk with centers µ1 , . . . , µk ∈ Rp as the solution to: k X X min kxi − µl k2 . (25) 1 ,...,Sk partition S µ1 ,...,µk l=1 i∈Si Figure 13: The d = 2 example of the use of this method to the example described above, the value of the nodes is given by color coding. For d = 2 it appears to smoothly interpolate between the labeled points. Note that, given the partition, the optimal centers are given by µl = 1 X xi . |Sl | i∈Sl Lloyd’s algorithm [Llo82] (also known as the k-means algorithm), is an iterative algorithm that alternates between • Given centers µ1 , . . . , µk , assign each point xi to the cluster l = argminl=1,...,k kxi − µl k . • Update the centers µl = 1 |Sl | P i∈Sl xi . Unfortunately, Lloyd’s algorithm is not guaranteed to converge to the solution of (25). Indeed, Lloyd’s algorithm oftentimes gets stuck in local optima of (25). A few lectures from now we’ll discuss convex relaxations for clustering, which can be used as an alternative algorithmic approach to Lloyd’s algorithm, but since optimizing (25) is N P -hard there is not polynomial time algorithm that works in the worst-case (assuming the widely believed conjecture P = 6 NP) While popular, k-means clustering has some potential issues: • One needs to set the number of clusters a priori (a typical way to overcome this issue is by trying the algorithm for different number of clusters). • The way (25) is defined it needs the points to be defined in an Euclidean space, oftentimes we are interested in clustering data for which we only have some measure of affinity between different data points, but not necessarily an embedding in Rp (this issue can be overcome by reformulating (25) in terms of distances only). Figure 14: The d = 3 example of the use of this method to the example described above, the value of the nodes is given by color coding. For d = 3 the solution appears to only learn the label −1. Figure 15: The d = 3 example of the use of this method with the extra regularization f T L2 f to the example described above, the value of the nodes is given by color coding. The extra regularization seems to fix the issue of discontinuities. • The formulation is computationally hard, so algorithms may produce suboptimal instances. • The solutions of k-means are always convex clusters. This means that k-means may have difficulty in finding cluster such as in Figure 17. 3.2 Spectral Clustering A natural way to try to overcome the issues of k-means depicted in Figure 17 is by using Diffusion Maps: Given the data
points we construct a weighted graph G = (V, E, W ) using a kernel K , such as 1 2 K (u) = exp 2 u , by associating each point to a vertex and, for which pair of nodes, set the edge weight as wij = K (kxi − xj k) . Figure 16: Examples of points separated in clusters. Recall the construction of a matrix M = D−1 W as the transition matrix of a random walk wij = Mij , Prob {X(t + 1) = j|X(t) = i} = deg(i) where D is the diagonal with Dii = deg(i). The d-dimensional Diffusion Maps is given by   λt2 ϕ2 (i)   (d) .. φt (i) =  , . λtd+1 ϕd+1 (i) where M = ΦΛΨT where Λ is the diagonal matrix with the eigenvalues of M and Φ and Ψ are, respectively, the right and left eigenvectors of M (note that they form a bi-orthogonal system, ΦT Ψ = I). If we want to cluster the vertices of the graph in k clusters, then it is natural to truncate the Diffusion Map to have k − 1 dimensions (since in k − 1 dimensions we can have k linearly separable sets). If indeed the clusters were linearly separable after embedding then one could attempt to use k-means on the embedding to find the clusters, this is precisely the motivation for Spectral Clustering. Algorithm 3.1 (Spectral Clustering) Given a graph G = (V, E, W ) and a number of clusters k (and t), Spectral Clustering consists in taking a (k − 1) dimensional Diffusion Map  t  λ2 ϕ2 (i)   (k−1) .. φt (i) =   . λtk ϕk (i) (k−1) and clustering the points φt tering. (k−1) (1), φt (k−1) (2), . . . , φt (n) ∈ Rk−1 using, for example, k-means clus- Figure 17: Because the solutions of k-means are always convex clusters, it is not able to handle some cluster structures. 3.3 Two clusters We will mostly focus in the case of two cluster (k = 2). For k = 2, Algorithm 3.1 consists in assigning to each vertex i a real number ϕ2 (i) and then clustering the points in the real line. Note in R, clustering reduces to setting a threshold τ and taking S = {i ∈ V : ϕ2 (i) ≤ τ }. Also, it is computationally tractable to try all possible thresholds (there are ≤ n different possibilities). Figure 18: For two clusters, spectral clustering consists in assigning to each vertex i a real number ϕ2 (i), then setting a threshold τ and taking S = {i ∈ V : ϕ2 (i) ≤ τ }. Algorithm 3.2 (Spectral Clustering for two clusters) Given a graph G = (V, E, W ), consider the two-dimensional Diffusion Map i → ϕ2 (i). set a threshold τ (one can try all different possibilities) and set S = {i ∈ V : ϕ2 (i) ≤ τ }. In what follows we’ll give a different motivation for Algorithm 3.2. 3.3.1 Normalized Cut Given a graph G = (V, E, W ), a natural measure to measure a vertex partition (S, S c ) is XX cut(S) = wij . i∈S j∈S c Note however that the minimum cut is achieved for S = ∅ (since cut(∅) = 0) which is a rather meaningless choice of partition. Remark 3.3 One way to circumvent this issue is to ask that |S| = |S c | (let’s say that the number of vertices n = |V | is even), corresponding to a balanced partition. We can then identify a partition with a label vector y ∈ P {±1}n where yi = 1 is i ∈ S, and yi = −1 otherwise. Also, the balanced condition can be written as ni=1 yi = 0. This means that we can write the minimum balanced cut as min cut(S) = S⊂V |S|=|S c | 1X 1 min y T LG y, wij (yi − yj )2 = 4 4 y∈{−1,1}n min y ∈{−1,1}n 1T y=0 i≤j 1T y=0 where LG = D − W is the graph Laplacian.13 . Since asking for the partition to be balanced is too restrictive in many cases, there are several ways to evaluate a partition that are variations of cut(S) that take into account the intuition that one wants both S and S c to not be too small (although not necessarily equal to |V |/2). A prime example is Cheeger’s cut. Definition 3.4 (Cheeger’s cut) Given a graph and a vertex partition (S, S c ), the cheeger cut (also known as conductance, and sometimes expansion) of S is given by h(S) = cut(S) , min{vol(S), vol(S c )} P where vol(S) = i∈S deg(i). Also, the Cheeger’s constant of G is given by hG = min h(S). S⊂V A similar object is the Normalized Cut, Ncut, which is given by Ncut(S) = cut(S) cut(S c ) + . vol(S) vol(S c ) Note that Ncut(S) and h(S) are tightly related, in fact it is easy to see that: h(S) ≤ Ncut(S) ≤ 2h(S). 13 W is the matrix of weights and D the degree matrix, a diagonal matrix with diagonal entries Dii = deg(i). Both h(S) and Ncut(S) favor nearly balanced partitions, Proposition 3.5 below will give an interpretation of Ncut via random walks. Let us recall the construction form previous lectures of a random walk on G = (V, E, W ): wij Prob {X (t + 1) = j|X(t) = i} = = Mij , deg(i) where M = D−1 W . Recall that M = ΦΛΨT where Λ is the diagonal matrix with the eigenvalues λk of M and Φ and Ψ form a biorthogonal system ΦT Ψ = I and correspond to, respectively, the right 1 1 and left eigenvectors of M . Moreover they are given by Φ = D− 2 V and Ψ = D 2 V where V T V = I 1 1 1 1 and D− 2 W D− 2 = V ΛV T is the spectral decomposition of D− 2 W D− 2 . Recall also that M 1 = 1, corresponding to M ϕ1 = ϕ1 , which means that ψ1T M = ψ1T , where 1 ψ1 = D 2 v1 = Dϕ1 = [deg(i)]1≤i≤n . This means that h i deg(i) vol(G) 1≤i≤n is the stationary distribution of this random walk. Indeed it is easy to check that, if X(t) has a certain distribution pt then X(t + 1) has a distribution pt+1 given by pTt+1 = pTt M Proposition 3.5 Given a graph G = (V, E, W ) and a partition (S, S c ) of V , Ncut(S) corresponds to the probability, in the random walk associated with G, that a random walker in the stationary distribution goes to S c conditioned on being in S plus the probability of going to S condition on being in S c , more explicitly: Ncut(S) = Prob {X(t + 1) ∈ S c |X(t) ∈ S} + Prob {X(t + 1) ∈ S|X(t) ∈ S c } , where Prob{X(t) = i} = deg(i) vol(G) . Proof. Without loss of generality we can take t = 0. Also, the second term in the sum corresponds to the first with S replaced by S c and vice-versa, so we’ll focus on the first one. We have: Prob {X(1) ∈ S c |X(0) ∈ S } = = = Prob {X(1) ∈ S c ∩ X(0) ∈ S } Prob {X(0) ∈ S} P P i∈S j∈S c Prob {X(1) ∈ j ∩ X(0) ∈ i} P i∈S Prob {X(0) = i} P P deg(i) wij P P = j∈S c vol(G) deg(i) deg(i) i∈S vol(G) i∈S i∈S P P i∈S = j∈S c wij deg(i) cut(S) . vol(S) Analogously, Prob {X(t + 1) ∈ S|X(t) ∈ S c } = which concludes the proof. cut(S) , vol(S c ) 2 3.3.2 Normalized Cut as a spectral relaxation Below we will show that Ncut can be written in terms of a minimization of a quadratic form involving the graph Laplacian LG , analogously to the balanced partition. Recall that balanced partition can be written as 1 min y T LG y. 4 y∈{−1,1}n 1T y=0 An intuitive way to relax the balanced condition is to allow the labels y to take values in two different real values a and b (say yi = a if i ∈ S and yj = b if i ∈ / S) but not necessarily ±1. We can then use the notion of volume of a set to ensure a less restrictive notion of balanced by asking that a vol (S) + b vol (S c ) = 0, which corresponds to 1T Dy = 0. We also need to fix a scale/normalization for a and b: a2 vol (S) + b2 vol (S c ) = 1, which corresponds to y T Dy = 1. This suggests considering min y∈{a,b}n 1T Dy=0, y T Dy=1 y T LG y. As we will see below, this corresponds precisely to Ncut. Proposition 3.6 For a and b to satisfy a vol (S) + b vol (S c ) = 0 and a2 vol (S) + b2 vol (S c ) = 1 it must be that  1  1 2 2 vol(S c ) vol(S) a= and b = − , vol(S) vol(G) vol(S c ) vol(G) corresponding to yi =      −   vol(S c ) vol(S) vol(G) vol(S) vol(S c ) vol(G) 1 2 1 2 if i ∈ S if i ∈ S c . Proof. The proof involves only doing simple algebraic manipulations together with noticing that vol(S) + vol(S c ) = vol(G). 2 Proposition 3.7 Ncut(S) = y T LG y, where y is given by yi =      −   vol(S c ) vol(S) vol(G) vol(S) vol(S c ) vol(G) 1 2 1 2 if i ∈ S if i ∈ S c . Proof. 1X wij (yi − yj )2 2 i,j XX = wij (yi − yj )2 y T LG y = i∈S j∈S c = = = = = = 1   1 #2 2 2 vol(S c ) vol(S) + wij c vol(S) vol(G) vol(S ) vol(G) i∈S j ∈S c   XX 1 vol(S c ) vol(S) wij + +2 vol(S c ) vol(G) vol(S) c i∈S j ∈S   XX 1 vol(S c ) vol(S) vol(S) vol(S c ) wij + + + vol(G) vol(S) vol(S c ) vol(S) vol(S c ) i∈S j∈S c   XX 1 1 wij + vol(S) vol(S c ) i∈S j∈S c   1 1 + cut(S) vol(S) vol(S c ) Ncut(S). XX " 2 This means that finding the minimum Ncut corresponds to solving min y T LG y s. t. y ∈ {a, b}n for some a and b y T Dy = 1 y T D1 = 0. (26) Since solving (26) is, in general, NP-hard, we consider a similar problem where the constraint that y can only take two values is removed: min y T LG y s. t. y ∈ Rn y T Dy = 1 y T D1 = 0. (27) Given a solution of (27) we can round it to a partition by setting a threshold τ and taking S = {i ∈ V : yi ≤ τ }. We will see below that (27) is an eigenvector problem (for this reason we call (27) a spectral relaxation) and, moreover, that the solution corresponds to y a multiple of ϕ2 meaning that this approach corresponds exactly to Algorithm 3.2. In order to better see that (27) is an eigenvector problem (and thus computationally tractable), 1 1 1 set z = D 2 y and LG = D− 2 LG D− 2 , then (27) is equivalent min z T LG z s. t. z ∈ Rn kzk2 = 1  1 T D 2 1 z = 0. 1 (28) 1 Note that LG = I − D− 2 W D− 2 . We order its eigenvalues in increasing order 0 = λ1 (LG ) ≤ 1 λ2 (LG ) ≤ · · · ≤ λn (LG ). The eigenvector associated to the smallest eigenvector is given by D 2 1 this means that (by the variational interpretation of the eigenvalues) that the minimum of (28) is λ2 (LG ) 1 1 and the minimizer is given by the second smallest eigenvector of LG = I − D− 2 W D− 2 , which is the 1 1 second largest eigenvector of D− 2 W D− 2 which we know is v2 . This means that the optimal y in (27) 1 is given by ϕ2 = D− 2 v2 . This confirms that this approach is equivalent to Algorithm 3.2. Because the relaxation (27) is obtained from (26) by removing a constraint we immediately have that λ2 (LG ) ≤ min Ncut(S). S⊂V This means that 1 λ2 (LG ) ≤ hG . 2 In what follows we will show a guarantee for Algorithm 3.2. Lemma 3.8 There is a threshold τ producing a partition S such that p h(S) ≤ 2λ2 (LG ). This implies in particular that h(S) ≤ p 4hG , meaning that Algorithm 3.2 is suboptimal at most by a square root factor. Note that this also directly implies the famous Cheeger’s Inequality Theorem 3.9 (Cheeger’s Inequality) Recall the definitions above. The following holds: p 1 λ2 (LG ) ≤ hG ≤ 2λ2 (LG ). 2 Cheeger’s inequality was first established for manifolds by Jeff Cheeger in 1970 [Che70], the graph version is due to Noga Alon and Vitaly Milman [Alo86, AM85] in the mid 80s. The upper bound in Cheeger’s inequality (corresponding to Lemma 3.8) is more interesting but more difficult to prove, it is often referred to as the “the difficult part” of Cheeger’s inequality. We will prove this Lemma in what follows. There are several proofs of this inequality (see [Chu10] for four different proofs!). The proof that follows is an adaptation of the proof in this blog post [Tre11] for the case of weighted graphs. Proof. [of Lemma 3.8] We will show that given y ∈ Rn satisfying R(y) := y T LG y ≤ δ, y T Dy and y T D1 = 0. there is a “rounding of it”, meaning a threshold τ and a corresponding choice of partition S = {i ∈ V : yi ≤ τ } such that h(S) ≤ √ 2δ, since y = ϕ2 satisfies the conditions and gives δ = λ2 (LG ) this proves the Lemma. We will pick this threshold at random and use the probabilistic method to show that at least one of the thresholds works. First we can, without loss of generality, assume that y1 ≤ · ≤ yn (we can simply relabel the vertices). Also, note that scaling of y does not change the value of R(y). Also, if y D 1 = 0 adding a multiple of 1 to y can only decrease the value of R(y): the numerator does not change and the denominator (y + c1)T D(y + c1) = y T Dy + c2 1T D1 ≥ y T Dy. This means that we can construct (from y by adding a multiple of 1 and scaling) a vector x such that x1 ≤ ... ≤ xn , xm = 0, and x21 + x2n = 1, and x T LG x ≤ δ, xT Dx where m be the index for which vol({1, . . . , m−1}) ≤ vol({m, . . . , n}) but vol({1, . . . , m}) > vol({m, . . . , n}). We consider a random construction of S with the following distribution. S = {i ∈ V : xi ≤ τ } where τ ∈ [x1 , xn ] is drawn at random with the distribution Z b Prob {τ ∈ [a, b]} = 2|τ |dτ, a where x1 ≤ a ≤ b ≤ xn . It is not difficult to check that  Prob {τ ∈ [a, b]} = b2 − a2 a2 + b2 if a and b have the same sign if a and b have different signs Let us start by estimating E cut(S). 1 XX E cut(S) = E wij 1(S,S c ) cuts the edge (i,j) 2 i∈V j ∈V = 1 XX 2 wij Prob{(S, S c ) cuts the edge (i, j)} i∈V j ∈V Note that Prob{(S, S c ) cuts the edge (i, j)} is x2i − x2j is xi and xj have the same sign and x2i +xj2 otherwise. Both cases can be conveniently upper bounded by |xi − xj | (|xi | + |xj |). This means that 1X E cut(S) ≤ wij |xi − xj | (|xi | + |xj |) 2 i,j s sX 1 X 2 ≤ wij (xi − xj ) wij (|xi | + |xj |)2 , 2 ij ij where the second inequality follows from the Cauchy-Schwarz inequality. From the construction of x we know that X wij (xi − xj )2 = 2xT LG x ≤ 2δxT Dx. ij Also, ! X wij (|xi | + |xj |)2 ≤ ij X wij 2x2i + 2xj2 . = 2 ij X deg(i)x2i   X + 2 deg(j)xj2  = 4xT Dx. i This means that E cut(S) ≤ j √ √ 1√ 2δxT Dx 4xT Dx = 2δ xT Dx. 2 On the other hand, c E min{vol S, vol S } = n X deg(i) Prob{xi is in the smallest set (in terms of volume)}, i=1 c vol(S ) we to break ties, if vol(S) = take the “smallest” set to be the one with the first indices. Note that m is always in the largest set. Any vertex j < m is in the smallest set if xj ≤ τ ≤ xm = 0 and any j > m is in the smallest set if 0 = xm ≤ τ ≤ xj . This means that, Prob{xi is in the smallest set (in terms of volume) = x2j . Which means that c E min{vol S, vol S } = n X deg(i)x2i = xT Dx. i=1 Hence, Note however that because we do not necessarily have √ E cut(S) ≤ 2δ. c E min{vol S, vol S } E cut(S) E min{vol S,vol S c } cut(S) is not necessarily the same as E min{vol S,vol S c } and so, √ cut(S) ≤ 2δ. min{vol S, vol S c } However, since both random variables are positive, √ E cut(S) ≤ E min{vol S, vol S c } 2δ, E or equivalently h √ i E cut(S) − min{vol S, vol S c } 2δ ≤ 0, which guarantees, by the probabilistic method, the existence of S such that √ cut(S) ≤ min{vol S, vol S c } 2δ, which is equivalent to h(S) = which concludes the proof of the Lemma. √ cut(S) ≤ 2δ, c min{vol S, vol S } 2 3.4 Small Clusters and the Small Set Expansion Hypothesis We now restrict to unweighted regular graphs G = (V, E). Cheeger’s inequality allows to efficiently approximate its Cheeger number up to a square root factor. It means in√particular that, given G = (V, E) and φ we can efficiently between the cases where: hG ≤ φ or hG ≥ 2 φ. Can this be improved? Open Problem 3.1 Does there exists a constant c > 0 such that it is N P -hard to, given φ, and G distinguis between the cases 1. hG ≤ φ, and √ 2. hG ≥ c φ? It turns out that this is a consequence [RST12] of an important conjecture in Theoretical Computer Science (see [BS14] for a nice description of it). This conjecture is known [RS10] to imply the UniqueGames Conjecture [Kho10], that we will discuss in future lectures. Conjecture 3.10 (Small-Set Expansion Hypothesis [RS10]) For every  > 0 there exists δ > 0 such that it is N P -hard to distinguish between the cases 1. There exists a subset S ⊂ V with vol(S) = δ vol(V ) such that 2. 3.5 cut(S) vol(S ) cut(S ) vol(S) ≤ , ≥ 1 − , for every S ⊂ V satisfying vol(S) ≤ δ vol(V ). Computing Eigenvectors Spectral clustering requires us to compute the second smallest eigenvalue of LG . One of the most efficient ways of computing eigenvectors is through the power method. For simplicity we’ll consider the case on which we are computing the leading eigenvector of a matrix A ∈ Rn×n with m nonzero entries, for which |λmax (A)| ≥ |λmin (A)| (the idea is easily adaptable). The power method Ay t proceeds by starting with a guess y 0 and taking iterates y t+1 = kAy t k . One can show [KW92] that the
variantes of the power method can find a vector x in randomized time O δ −1 (m + n) log n satisfying xT Ax ≥ λmax (A)(1 − δ)xT x. Meaning that an approximate solution can be found in quasi-linear time.14 One drawback of the power method is that when using it, one cannot be sure, a posteriori, that there is no eigenvalue of A much larger than what we have found, since it could happen that all our guesses were orthogonal to the corresponding eigenvector. It simply guarantees us that if such an eigenvalue existed, it would have been extremely likely that the power method would have found it. This issue is addressed in the open problem below. Open Problem 3.2 Given a symmetric matrix M with small condition number, is there a quasilinear time (on n and the number of non-zero entries of M ) procedure that certifies that M  0. More specifically, the procedure can be randomized in the sense that it may, with some probably not certify that M  0 even if that is the case, what is important is that it never produces erroneous certificates (and that it has a bounded-away-from-zero probably of succeeding, provided that M  0). 14 Note that, in spectral clustering, an error on the calculation of ϕ2 propagates gracefully to the guarantee given by Cheeger’s inequality. The Cholesky decomposition produces such certificates, but we do not know how to compute it in quasi-linear time. Note also that the power method can be used in αI − M to produce certificates that have arbitrarily small probability of being false certificates. Later in these lecture we will discuss the practical relevance of such a method as a tool to quickly certify solution produced by heuristics [Ban15b]. 3.6 Multiple Clusters Given a graph G = (V, E, W ), a natural way of evaluating k-way clusterign is via the k-way expansion constant (see [LGT12]):   cut(S) ρG (k) = min max , S1 ,...,Sk l=1,...,k vol(S) where the maximum is over all choice of k disjoin subsets of V (but not necessarily forming a partition). Another natural definition is ϕG (k) = cut(S) . vol(G) vol(S) min S:vol S≤ k1 It is easy to see that ϕG (k) ≤ ρG (k). The following is known. Theorem 3.11 ([LGT12]) Let G = (V, E, W ) be a graph and k a positive integer
p λk , ρG (k) ≤ O k 2 (29) Also, ρG (k) ≤ O p  λ2k log k . Open Problem 3.3 Let G = (V, E, W ) be a graph and k a positive integer, is the following true? p ρG (k) ≤ polylog(k) λk . (30) We note that (30) is known not to hold if we ask that the subsets form a partition (meaning that every vertex belongs to at least one of the sets) [LRTV12]. Note also that no dependency on k would contradict the Small-Set Expansion Hypothesis above.