這次介紹的是Alex和Alessandro于2014年發表在的Science上的一篇關于聚類的文章[13]，該文章的基本思想很簡單，但是其聚類效果卻兼具了譜聚類(Spectral Clustering)[11,14,15]和K-Means的特點，著實激起了我的極大的興趣，該聚類算法主要是基于兩個基本點：

1. 聚類中心的密度高于其臨近的樣本點的密度
2. 聚類中心與比其密度還高的聚類中心的距離相對較大

基于這個思想，聚類過程中的聚類中心數目可以很直觀的選取，離群點也能被自動檢測出來并排除在聚類分析外。無論每個聚類的形狀是什么樣的，或者樣本點的維度是多少，聚類分析的結果都能令人很滿意。下面我會主要基于這篇文章來詳述該聚類算法的來龍去脈，并簡單回顧下相關的聚類算法。

聚類算法回顧

眾所周知，聚類分析目的在于根據樣本之間的相似性將樣本劃為不同的類簇，但聚類的科學定義貌似在學術界還未達成共識。論文[17]對聚類算法進行了綜述，發現有好多聚類算法我還沒學習。在K-means和K-medoids中，每個類簇都由一組到各自的聚類中心距離最近的數據組成。兩者的目標函數形式為各樣本點到對應的聚類中心的距離之和，經過反復的更新聚類中心和重新為樣本點分配聚類中心的過程直至收斂，如圖1所示。兩者的區別在于，K-means的聚類中心為屬于該類簇的所有樣本點的均值，而K-medoids的聚類中心為該類簇中離所有樣本點的聚類之和最小的樣本點。這兩者聚類算法實現起來都非常簡單，對于緊湊型的和呈超球體狀分布的數據非常適用。但兩者的缺陷也很明顯：

1. 缺乏能確定類簇數目和進行初步劃分的有效機制；
2. 迭代優化的策略無法保證全局最優解；
3. 對離群點和噪聲非常敏感

　　在基于概率密度的聚類算法中，我們會假設各類簇由不同的概率密度函數產生(如圖2)，而每個樣本點則是以不同的權重服從這些概率分布的。很不幸的是，在這類算法中用最大似然估計求解參數往往不可行，只能用迭代求解的方式獲得一個次優解，而期望最大化(Expectation Maximization,EM)是最常用的一個策略。在這類算法中，最典型的莫過于高斯混合模型(Gaussian Mixture Model,GMM)[12]。這類算法的準確度取決于預先定義的概率分布能否很好的擬合訓練數據，但問題在于很多情況下我們無法知曉數據在整體上或者局部上到底近似于什么樣的概率分布。 ?

　　基于局部密度的聚類算法可以很容易地檢測出任意形狀的類簇。在DBSCAN[10]中，需要用戶給定密度閾值和領域半徑作為參數，在領域半徑內的密度小于該閾值的樣本點被視為噪聲點，剩下的密度較高的非連通區域則被分配到不同的類簇中，其偽代碼如下所示。但是選擇合適的密度閾值并不是那么容易的事情，有關的參數估計建議可參見[3]。DBSCAN的優點[3]總結如下：

• 無需預先指定類簇的數目；
• 可以發現任意形狀的類簇，如圖3所示；
• 可以檢測出噪聲點，且對噪聲點魯棒性較強；
• 除了邊界點外，聚類結果(核心點與噪聲點)與樣本點的遍歷順序無關

DBSCAN的缺點[3]總結如下：

• 針對邊界點而言，DBSCAN的聚類結果并非完全確定的。幸運的是這種情況并非頻繁出現，而且對聚類的結果影響很小。如果把邊界點也當成噪聲點處理，那么聚類結果就具有確定性。
• 聚類結果依賴于距離度量規則。最常用的歐式距離在高維空間里由于“維度災難”幾乎無法發揮有效作用，使得設定合適的搜尋半徑更為困難。
• 不適用于密度差異很大的數據集，因為此時各個類簇的搜尋半徑和密度閾值都不相同，使得參數的選取更為困難。

    DBSCAN(D, eps, minPts)
//eps:search radius
//minPts:density threshold
   C = 0
   for each unvisited point P in dataset D
      mark P as visited
      NeighborPts = regionQuery(P, eps)
      if sizeof(NeighborPts) < minPts
         mark P as NOISE
      else
         C = next cluster
         expandCluster(P, NeighborPts, C, eps, MinPts)
          
expandCluster(P, NeighborPts, C, eps, minPts)
   add P to cluster C
   for each point Q in NeighborPts 
      if Q is not visited
         mark Q as visited
         NeighborPts' = regionQuery(Q, eps)
         if sizeof(NeighborPts') >= minPts
            NeighborPts = NeighborPts joined with NeighborPts'
      if Q is not yet member of any cluster
         add Q to cluster C
          
regionQuery(P, eps)
   return all points within P's eps-neighborhood (including P)
  

　　基于均值漂移(Mean-shift)[5,7,9]的聚類算法則無需為搜索半徑和密度閾值的設定而煩惱，不過也面臨bandwidth的選取問題，關于怎么設定bandwidth的研究可參見[8,16]。Mean-sift的基本思路就是從初始點出發，以梯度上升的方式不斷尋找核密度估計函數的局部最大值直至收斂(如圖4(a)所示)，這些駐點代表分布的模式。在基于mean-shift的聚類算法中，依次以每一個樣本點作為mean-shift的起始點，然后將其移至核密度估計函數的某個局部駐點，最后近似收斂到同一個駐點的所有樣本被劃分至同一個類簇，如圖4(b)所示。總體而言，在基于密度的聚類算法中，類簇可被定義為收斂到相同的密度分布函數局部極大值的樣本點的集合。?

基于密度峰值和距離的聚類算法

該聚類算法的假設前提是聚類中心周圍的樣本點的局部密度低于聚類中心的局部密度，并且聚類中心與比其局部密度更高的點之間的距離相對較大。其聚類效果與DBSCAN和mean-shift類似，可以檢測出非球體的類簇。作者號稱可以 自動 找到類簇的數目，雖然文中給了一點相關的尋找聚類數目的思路，但是提供的Matlab代碼中沒有實現該思路，還是需要人工選擇聚類中心，所以在相關評論[2]中“自動”一詞遭到了質疑。與mean-shift類似，聚類中心定義為局部密度最大值點；與mean-shift不同的是，聚類中心是某個特定樣本點，并且無需在核函數定義的空間內針對每個樣本點顯式求解局部密度最大的點。 給定數據集\(\mathcal{S}=\{x_i|x_i\in\mathbb{R}^n,i=1,\cdots,N\}\)，對于每一個樣本點\(x_i\)計算兩個量化值：局部密度值\(\rho_i\)和距離密度更高的樣本點的聚類\(\delta_i\)。\(x_i\)的局部密度\(\rho_i\)定義如為： \begin{equation} \rho_i=\sum_{j=1}^N\chi(d_{ij}-d_c) \end{equation} 其中\(d_c\)為截斷距離(cutoff distance)，其實就是領域的搜索半徑；\(d_{ij}\)為\(x_i\)與\(x_j\)之間的距離；函數\(\chi(x)\)定義為 \begin{equation} \chi(x)=\begin{cases} 1,& \text{if \(x<0\)};\\ 0,& otherwise. \end{cases} \end{equation} 根據這篇文章的評論[2],發現還有兩個密度的度量方法也是很有價值的.第一個是用樣本點與最近的\(M\)個鄰居的距離的均值的負數來描述;另一個就是高斯核函數來度量,會比用截斷距離度量魯棒性更強一些. \begin{equation} \rho(x_i)=-\frac{1}{M}\sum_{j:j\in KNN(x_i)}d_{ij} \label{eq:avg_kernel} \end{equation} \begin{equation} \rho(x_i)=\sum_{j=1}^N\exp(-\frac{d_{ij}^2}{\sigma}) \label{eq:gauss_kernel} \end{equation} 實際上，上述的\(\rho_i\)定義的就是與\(x_i\)之間的距離小于\(d_c\)的樣本點的數目。距離\(\delta_i\)度量\(x_i\)與比其密度高的最近的樣本點之間的距離；如果\(\rho_i\)為最大值，則\(\delta_i\)為與離\(x_i\)最遠樣本之間的距離： \begin{equation} \delta_i=\begin{cases} \underset{j:\rho_j>\rho_i}{\min}(d_{ij}), & \text{if \(\exists j,\rho_j>\rho_i\)};\\ \underset{j}{\max}(d_{ij}), & otherwise. \end{cases} \end{equation} 對于密度值為局部或全局最大的樣本點而言，它們的\(\delta_i\)會比其他樣本點的\(\delta_j\)值要大很多(如圖5所示)，因為前者代表局部密度最大的樣本點之間的距離，而后者代表樣本點與其對應的局部密度最大的樣本點之間的距離。因此，那些\(\delta\)值很大的樣本點也很有可能就是聚類中心。?

　　論文中給出了一個示例，如圖6(a)所示，圖中一共有28個樣本點，樣本點按照密度降序排列。從圖中大致可以觀察到有兩個類簇，剩下的26、27和28號樣本點可被視為離群點。在圖6(b)中，分別以對判定是否為聚類中的最關鍵信息\(\rho\)和\(\delta\)為橫縱坐標繪制決策圖(decision graph)，看到1和10號樣本點位于決策圖的最右上角。9和10號樣本點雖然密度值\(\rho\)非常接近，但是\(\delta\)值卻相差很大；被孤立的26、27和28號樣本點雖然\(\delta\)值較大，但是\(\rho\)值很小。綜上可知，只有\(\rho\)值很高并且\(\delta\)相對較大的樣本點才會是聚類中心。

?

　　在找出聚類中心后，接下來就是將所有剩下的點劃分到比其密度更高且最近的樣本點所屬的類簇中，當然經過這一步之后暫時會為噪聲點也分配到類簇中。在聚類分析中，經常還會進一步分析類簇分配的可靠性。在DBSCAN中，只考慮了密度高于密度閾值的可靠性高一些的樣本點，但是會出現較低密度的類簇被誤認為噪聲的情況。文中取而代之的是為每個類簇引入邊界區域的概念。邊界區域的密度值\(\rho_b\)會根據屬于這個類簇并且與屬于其他類簇的樣本點之間的距離小于\(d_c\)的成員計算出來。對于每個類簇中的所有樣本點，密度值高于\(\rho_b\)的被視為類簇的核心組成部分(cluster core)，剩下的則被視為該類簇的光暈(cluster halo)，類簇光暈中則包含噪聲點。論文中給出了一個聚類的結果，如圖7所示。

?

　　鄰域搜索半徑\(d_c\)到底如何取值呢？\(d_c\)顯然是對聚類結果又影響的，這一點我們僅需要考慮兩個最極端的情形就明白了。如果\(d_c\)太大，那么每個數據點的密度值都近似相等，導致所有數據點被劃分至同一個類簇中；如果\(d_c\)太小，每個類簇包含的樣本點會很少，很有可能出現同一個類簇被分割成好幾部分的情況。另一方面，不同的數據集中數據點之間的密集程度不同，那么想給出一個適合所有數據集的\(d_c\)是不可能的。作者在文中提出，合適的\(d_c\)應該使數據點的平均近鄰數目占整個數據集規模的比例為\(\tau,(\tau=1\%\sim 2\%)\)。如此一來，參數\(\tau\)就獨立于特定數據集了。針對每個數據集，我們都可以尋找一個比較合適的\(d_c\)。結合作者給出的Matlab代碼，分析后可知具體的計算方法如下：取出對稱的距離矩陣的上三角所有的\(M=N(N-1)/2\)個元素，然后對其進行升序排列\(d_1\leq d_2\leq \cdots\leq d_M\)。為了保證平均每個數據點的近鄰點數目所占比例為\(\tau\)，那么只要保證小于\(d_c\)的距離數目所占比例也為\(\tau\)即可，因此取\(d_c=d_{round(\tau M)}\)。 類簇的數目該如何確定呢？作者給Matlab代碼中，聚類中心是需要人工選定的，很多讀者因此質疑文中的"it is able to detect nonspherical clusters and to automatically find the correct number of clusters"，是不是有種被欺騙的感覺。不過作者在文中也給出了一個簡單選擇類簇的數目，雖然我也覺得該方法存在一些問題，但總歸還是給出了解決方案的。由前面解釋的論文的兩個基本立足點可知，聚類中心對應的\(\rho\)和\(\delta\)都是比較大的。作者為每個樣本點\(x_i\)引入\(\gamma_i=\rho_i\delta_i\)，然后將所有的\(\gamma_i\)降序排列后顯示在圖9(a)。 如果分別對\(\rho\)和\(\delta\)先做歸一化處理后會更合理一些，這樣也會使得兩者參與決策的權重相當。 因為如果\(\rho\)和\(\delta\)的不在一個數量級，那么必然數量級小帶來的的影響會很小。 接下來怎么辦呢？作者依然沒有給出具體的解決方案。因為整體而言，\(\gamma\)的值在大多數情況下還是很相近的，差異比較大的就是那幾個聚類中心，我覺得可以從異常檢查(Anomaly Detection)的角度去尋找這個跳躍點。最簡單方法，可以根據相鄰\(\gamma\)的值構建一個高斯分布\(\mathbb{N}(\mu,\sigma^2)\)，根據最大似然參數估計法，該高斯分布的參數只需掃描兩遍\(\gamma\)的值即可，所以模型還是很效率還是很高的。有了這個模型后，我們從后往前掃描\(\gamma\)的值，如果發現某個值的左邊或右邊的累積概率(如圖8的左右兩側藍色區域)小于閾值(比如0.005)時就判定找到了異常的跳躍點，此時就能大致確定類簇的數目了。若想進一步學習如何利用高斯分布進行異常檢測可參見[1]。 我們都知道高斯分布的概率密度函數,可是高斯分布的累積分布函數(Cumulative Distribution Function)不存在初等函數的表達形式,那該如何是好?查找了半天資料,也沒找到如何數值逼近的原理說明,不過搜到了一段用java編寫的基于 Hart Algorithm 近似計算標準正態分布的累積分布函數的代碼[4]。寥寥數行java代碼就搞定了,但是我暫時沒理解為什么這么做是可行的。我將其轉換成了如下的C++代碼,然后將輸出結果和維基百科上的Q函數表[6]中的數據對比分析(注意\(1-Q(x)=\Phi(x)\)),發現結果和預期的一模一樣,簡直把我驚呆了。

?

    double CDFofNormalDistribution(double x)
{
	const double PI=3.1415926;
	double p0=220.2068679123761;
	double p1=221.2135961699311;
	double p2=112.0792914978709;
	double p3=33.91286607838300;
	double p4=6.373962203531650;
	double p5=.7003830644436881;
	double p6=.03326249659989109;

	double q0=440.4137358247552;
	double q1=793.8265125199484;
	double q2=637.3336333788311;
	double q3=296.5642487796737;
	double q4=86.78073220294608;
	double q5=16.06417757920695;
	double q6=1.755667163182642;
	double q7=0.08838834764831844;

	double cutoff=7.071;//10/sqrt(2)
	double root2pi=2.506628274631001;//sqrt(2*PI)

	double xabs=abs(x);

	double res=0;
	if(x>37.0) 
		res=1.0;
	else if(x<-37.0)
		res=0.0;
	else
	{
		double expntl=exp(-.5*xabs*xabs);
		double pdf=expntl/root2pi;
		if(xabs<cutoff)
			res=expntl*((((((p6*xabs + p5)*xabs + p4)*xabs + p3)*xabs+ \
				p2)*xabs + p1)*xabs + p0)/(((((((q7*xabs + q6)*xabs + \
				q5)*xabs + q4)*xabs + q3)*xabs + q2)*xabs + q1)*xabs+q0);
		else
			res=pdf/(xabs+1.0/(xabs+2.0/(xabs+3.0/(xabs+4.0/(xabs+0.65)))));
	}
	if(x>=0.0)
		res=1.0-res;
	return res;
}
  

此外，作者聲稱根據隨機均勻分布生成的數據對應的\(\gamma\)服從冪律分布(Power laws)，但是真正具備聚類中心的數據集是不存在這種情況的。很多現象其實都是近似服從冪律分布的，尤其適用于大多數事件的規模很小但少數事件規模很大的場合，不過作者在此并未給出該定論的出處，所以同樣這一點遭到了很多讀者的質疑。我猜目前只是作者根據一些實驗數歸納出來的，只能說是靠不完全統計得到的經驗，沒有實質性的理論依據。也就是\(\gamma\approx cr^{-k}+\epsilon\)，其中\(r\)為\(\gamma\)的排名序號，那么\(\log\gamma\)和\(\log r\)之間應該近似呈現線性關系，如圖9(b)所示。如果作者的猜測正確的話，我們不妨在聚類前匯出如\(\log\gamma\)和\(\log r\)的關系圖，借此判斷聚類的復雜性，或者說在該數據集上進行聚類的結果可靠性如何。

?

　　 最后,基于這篇文章思想,我最終用C++代碼實現了一個比較完整的聚類算法,并作為我在GitHub上的first repository上傳到了GitHub上面,有需要的請前往 https://github.com/jeromewang-github/cluster-science2014 下載,歡迎大家找出bug和提供修改意見!

References

• [1] Anomaly detection. http://www.holehouse.org/mlclass/15_Anomaly_Detection.html.
• [2] Comments on clustering by fast search and find of density peaks. http://comments.sciencemag.org/content/10.1126/science.1242072.
• [3] Dbscan. http://en.wikipedia.org/wiki/DBSCAN.
• [4] Hart algorithm for normal cdf. http://www.onedigit.org/Home/quantitative-finance/hart-algorithm-for-normal-cdf.
• [5] Mean-shift. http://en.wikipedia.org/wiki/Mean-shift.
• [6] Q-function. http://en.wikipedia.org/wiki/Q-function.
• [7] Dorin Comaniciu and Peter Meer. Mean shift: A robust approach toward feature space analysis. Pattern Analysis and Machine Intelligence, IEEE-Transactions on, 24(5):603–619, 2002.
• [8] Dorin Comaniciu, Visvanathan Ramesh, and Peter Meer. The variable bandwidth mean shift and data-driven scale selection. In Computer Vision, 2001. ICCV 2001. Proceedings. Eighth IEEE International Conference on, volume 1, pages 438–445. IEEE, 2001.
• [9] Konstantinos G. Derpanis. Mean shift clustering. http://www.cse.yorku.ca/~kosta/CompVis_Notes/mean_shift.pdf, 2005.
• [10] Martin Ester, Hans-Peter Kriegel, J ?org Sander, and Xiaowei Xu. A density-based algorithm for discovering clusters in large spatial databases with noise. In Kdd, volume 96, pages 226–231, 1996.
• [11] Andrew Y Ng, Michael I Jordan, Yair Weiss, et al. On spectral clustering: Analysis and an algorithm. Advances in neural information?processing systems, 2:849–856, 2002.
• [12] Douglas Reynolds. Gaussian mixture models. Encyclopedia of Biometrics, pages 659–663, 2009.
• [13] Alex Rodriguez and Alessandro Laio. Clustering by fast search and find of density peaks. Science, 344(6191):1492–1496, 2014.
• [14] Aarti Singh. Spectral clustering. https://www.cs.cmu.edu/~aarti/Class/10701/slides/Lecture21_2.pdf.
• [15] Ulrike Von Luxburg. A tutorial on spectral clustering. Statistics and computing, 17(4):395–416, 2007.
• [16] Jue Wang, Bo Thiesson, Yingqing Xu, and Michael Cohen. Image?and video segmentation by anisotropic kernel mean shift. In Computer?Vision-ECCV 2004, pages 238–249. Springer, 2004.
• [17] Rui Xu, Donald Wunsch, et al. Survey of clustering algorithms. Neural Networks, IEEE Transactions on, 16(3):645–678, 2005.

Clustering by density peaks and distance