聚類之譜聚類（轉）

從樣本相似性到圖

根據我們一般的理解，聚類是將相似的樣本歸為一類，或者說使得同類樣本相似度盡量高，異類樣本相似性盡量低。無論如何，我們需要一個方式度量樣本間的相似性。常用的方式就是引入各種度量，如歐氏距離、余弦相似度、高斯度量等等。

度量的選擇提現了你對樣本或者業務的理解。比如說如果你要比較兩個用戶對音樂選擇的品味，考慮到有些用戶習慣打高分，有些用戶習慣打低分，那么選擇余弦相似度可能會比歐式距離更合理。

現在我們假設已有的樣本為X={x1,x2,…,xn}, 我們選擇的樣本相似性度量函數為(xi,xj)→s(xi,xj),其中s≥0,s越大表示相似性越高。一般我們選擇的是距離的倒數。據此我們可以構造相似性圖，節點為樣本，節點之間的連邊權重為樣本間的相似性，如圖所示。 

這是一個完全圖，我們的目的是去掉一些邊，使得這個圖變成K個連通子圖。同一個子圖內的節點歸為一類。因此有兩方面考慮：

• 子圖內的連邊權重盡量大，即同類樣本間盡量相似
• 去掉的邊權重盡量小，即異類樣本間盡量不同

一個初步的優化方法是去掉部分權重小的邊，常用的有兩種方式：

• ϵ准則，即去掉權重小於ϵ的所有邊
• k鄰近准則，即每個節點之保留權最大的幾條邊

現在我們得到一個較為稀疏的圖。 

圖與圖的Laplacian矩陣

為了下一步的算法推導，首先介紹圖的Laplacian矩陣，我們記節點i,j連邊的權重為wi,j,如果兩個節點之間沒有連邊，wi,j=0 ，另外wii=0，那么圖對應的Laplacian矩陣為： 

L(G,W)=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜∑j≠1nw1j−w21⋮−wn1−w1,2∑j≠2nw2j⋮−wn2⋯⋯⋱⋯−w1n−w2n⋮∑j≠nnwnj⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜∑j≠1nw1j∑j≠2nw2j⋱∑j≠nnwnj⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟−⎛⎝⎜⎜⎜⎜⎜w11w21⋮wn1w1,2w22⋮wn2⋯⋯⋱⋯w1nw2n⋮wnn⎞⎠⎟⎟⎟⎟⎟=D−W

 

容易看到，矩陣L行和為零,且對於任意的向量f有 

f′Lf=f′Df−f′Wf=∑i=1ndif2i−∑i,j=1nfifjwij=12⎛⎝∑i=1ndif2i−2∑i,j=1nfifjwij+∑j=1ndjf2j⎞⎠=12∑i,j=1nwij(fi−fj)2

 

優化目標

現在我們來推導我們要優化的目標函數。前面說過，我們的目的是去掉一些邊，使得這個圖變成K個連通子圖，我們還希望去掉的邊權重盡量小。為此，假設我們已經把圖分割成立K個連通子圖A1,⋯,AK,那么我們去掉的邊集合為

{ei,j|∃k,st.xi∈Ak and xj∉Ak}

為了方便，引入記號 

W(A,B)=∑i∈A,j∈Bwij

那么 

W(Ak,A¯k)=∑i∈Ak,j∉Akwij

因此去掉的邊的權重和為 

12∑k=1nW(Ak,A¯k)

 

現在的問題就轉換為：找到X的划分A1,⋯,AK,使得上式最小。不幸的是，存在兩個問題：

• 這是個NP難問題，沒有有效算法
• 實際實驗得到的結果常常將單獨的一個樣本分為一類

先來解決第二個問題： 
我們實際希望的是，每個類別中的樣本數要充分大，有兩種調整目標函數的方法：

1. RatioCut，將目標函數改成
   12∑k=1nW(Ak,A¯k)|Ak|
2. Ncut, 將目標函數改成
   12∑k=1nW(Ak,A¯k)vol（Ak）
   
   其中vol(A)=∑i∈Adi

兩種方法都使得某個類樣本量少的時候，對應的目標函數項變大。這里我們以第一種方法為例，第二種是類似的。

現在來解決第二個問題： 
我們碰到NP難問題的時候，通常是考慮近似解，譜聚類也不例外。首先，我們要引入列向量hk=(h1k,⋯,hn,k)′，其中 

hij=⎧⎩⎨1|Aj|√0xi∈Ajxi∉Aj

那么， 

h′kLhk=12∑i,j=1nwij(hkj−hkj)2=12∑xi∈Ak,xj∈Ak¯nwij(1|Ak|−−−√−0)2=12W(Ak,Ak¯)|Ak|

令H=(h1,…,hK),則H′H=I,且 

12∑k=1nW(Ak,A¯k)|Ak|=∑k=1nh′kLhk=tr(H′LH)

現在的問題是找到H使得tr(H′LH)盡量小。這當然還是NP問題。原因在於H的列向量元素只能選擇0或者|Ak|−−−√

 

這里用到的一個trick是放寬H的限制，在更大的變量空間內尋找目標函數的最小值。在這里，我們讓H可以取任意實數，只要滿足H′H=I.那么就是求解優化問題： 

argminH′H=Itr(H′LH)

 

令L=Q′ΛQ,Y=QH 則Y′Y=I, 

tr(H′LH)=tr((QH)′Λ(QH))=tr(Y′ΛY)=tr(YY′Λ)=∑i=1n(YY′)iiλi

 

由於Y′Y=I,有 

0≤(YY′)ii≤1

∑i=1n(YY′)ii=tr(YY′)=tr(Y′Y)=K

可以分析得到tr(H′LH)取最小值時，必有(YY′)ii=0 or 1, 
因此

tr(H′LH)≥∑i=1Kλi

且當Yk=ek時等號成立，對應的hk=(Q′Y)k=Q′ek=Q′k為L的第k小特征值對應的特征向量。

 

最后一步

現在我們得到了放寬限制條件下的優化問題的最優解解h1,⋯hK,如何得到原問題的一個解？

我們知道，如果H滿足原來的限制條件，那么hij≠0表示第i個樣本歸為第j類，但是我們放寬限制條件后得到的H可能一個零都沒有！

譜聚類有意思的地方是選擇了對H的行向量做K-means聚類，我個人認為是基於如下考慮：

1. 對滿足原始限制條件的H，行向量一致等價於類別一致
2. 在原始限制條件下得到的H跟放寬限制條件下得到的H應該比較相近

如此可以推斷在放寬條件下得到的H如果行向量相似，則應該類別相似。因此直接對行向量做k-means聚類即可。

總結

至此，譜聚類的大致步驟就完成了，歸納下主要步驟

1. 計算樣本相似性得到樣本為節點的完全圖
2. 基於ϵ准則或者m鄰近准則將完全圖稀疏化
3. 計算稀疏化后的圖的laplacian矩陣，計算其前K小特征值對應的特征向量作為矩陣H的列
4. 對矩陣H的行向量作k-means聚類
5. 若H的第i行與第j行被聚為一類，則表示xi與xj聚為一類。

代碼例子

左圖是原始數據，右圖是譜聚類結果 

 

import numpy as np
import networkx as nx
import scipy.linalg as llg
from Queue import PriorityQueue
import matplotlib.pylab as plt
import heapq as hp
from sklearn.cluster import k_means

# fake data from multivariate normal distribution
S = np.random.multivariate_normal([1,1], [[0.5,0],[0,0.7]],100)
T = np.random.multivariate_normal([-1,-1], [[0.3,0],[0,0.8]],100)
R = np.random.multivariate_normal([-1,0], [[0.4,0],[0,0.5]],100)
Data = np.vstack([S,T,R])
plt.subplot(1,2,1)
plt.scatter(S.T[0],S.T[1],c='r')
plt.scatter(T.T[0],T.T[1],c='b')
plt.scatter(R.T[0],R.T[1],c='y')

# calc k-nearest neighbors
def min_k(i,k):
    pq = []
    for j in range(len(Data)):
        if i == j:
            continue
        if len(pq) < k:
            hp.heappush( pq,(1/np.linalg.norm(Data[i]-Data[j]), j) )
        else:
            hp.heappushpop( pq, (1/np.linalg.norm(Data[i]-Data[j]), j) )
    return pq

# calc laplacian
L = np.zeros((len(Data),len(Data)))
for i in range(len(Data)):
    for (v,j) in min_k(i, 3):
        L[i,j] = -v
        L[j,i] = -v
L = L + np.diag(-np.sum(L,0)) 

# kmean
(lam, vec) = llg.eigh(L)
A = vec[:,0:3]
kmean = k_means(A,3)

plt.subplot(1,2,2)
plt.scatter(Data.T[0],Data.T[1],c=['r' if v==0 else 'b' if v==1 else 'y' for v in kmean[1]])
plt.show()       

 

 

轉：http://blog.csdn.net/betarun/article/details/51154003