1. 程式人生 > >K均值聚類的理解和實現

K均值聚類的理解和實現

目錄

1. 距離的測度

1.1 歐式距離

1.2 馬氏距離

1.2.1 利用馬氏距離對資料進行歸一化

1.2.2 利用馬氏距離進行分類

2. K均值的基本理論

2.1 K均值的原理和實現

2.2 K均值的缺點

2.3 K均值改進

3. 演算法實現

3.1 獲取樣本

3.2 協方差逆陣方根的計算方法

3.3 聚類實驗

3.3.1 一般的K均值聚類

3.3.2 基於馬氏距離K-means++聚類

3.3.3 基於肘部法則的K-means++聚類

4.參考資料


1. 距離測度

1.1 歐式距離

數學中歐氏距離歐幾里德度量

歐幾里得空間中兩點之間的“普通” 直線距離。通過這個距離,歐幾里德空間成為度量空間。相關的規範稱為歐幾里得範數較早的文獻將度量指為畢達哥拉斯度量廣義的歐幾里得範數項是L2範數L2距離

通常,對於n維空間來說,歐幾里得距離可以表示為:

R^{2}中的歐式距離如圖1.1-1所示:

                                                                           圖1.1-1 R^{2}

中歐幾里得距離的表達

標準的歐幾里德距離可以是平方的,以便逐漸將更大的重量放置在更遠的物體上。在這種情況下,等式變為:

                            {\ displaystyle d ^ {2}(\ mathbf {p},\ mathbf {q})=(p_ {1} -q_ {1})^ {2} +(p_ {2} -q_ {2})^ {2} + \ cdots +(p_ {i} -q_ {i})^ {2} + \ cdots +(p_ {n} -q_ {n})^ {2}。}

平方歐幾里德距離不是一個度量,因為它不滿足三角不等式 ; 然而,它經常用於僅需要比較距離的優化問題。

它在理性三角學領域也被稱為quadrance

1.2 馬氏距離

馬氏距離是對點P和分佈D的距離度量。馬氏距離對多維資料進行了歸一化,並測量了P點相對於D的平均值的標準差。如果P在D分佈中心,那麼馬氏距離為0。如果對資料進行主成分分析,如圖1.2-1所示,那麼,當P相對於主軸越遠,馬氏距離的數值也就隨之增長。當我們對主軸進行進行歸一化後,馬氏距離也就等同於在歐式空間的仿射變換。因此,馬氏距離具有“無單位”和“尺度不變性”的特性,並且考慮了資料集的相關性。

                                                                     圖1.2-1 資料的主成分分析

馬哈拉諾比斯觀察距離 :

{\ displaystyle {\ vec {x}} =(x_ {1},x_ {2},x_ {3},\ dots,x_ {N})^ {T}}

 從一組帶有均值的觀察中得出:

{\ displaystyle {\ vec {\ mu}} =(\ mu _ {1},\ mu _ {2},\ mu _ {3},\ dots,\ mu _ {N})^ {T}}

那麼,觀察值與集合的距離使用協方差矩陣S表示為:

{\ displaystyle D_ {M}ï¼{\ vec {x}}ï¼= {\ sqrt {ï¼{\ vec {x}}  -  {\ vec {\ mu}}ï¼^ {T} S ^ { -  1}ï¼ {\ vec {x}}  -  {\ vec {\ mu}}ï¼}}ã\ï¼}

集合中兩個隨機變數的距離為:

d(\ vec {x},\ vec {y})= \ sqrt {(\ vec {x}  -  \ vec {y})^ TS ^ { -  1}(\ vec {x}  -  \ vec {y} )}。\,

如果協方差矩陣是單位矩陣,則馬哈拉諾比斯距離減小到歐幾里得距離。如果協方差矩陣是對角矩陣,那麼得到的距離度量稱為標準歐式距離

d(\ vec {x},\ vec {y})= \ sqrt {\ sum_ {i = 1} ^ N {(x_i  -  y_i)^ 2 \ over s_ {i} ^ 2}},

其中,s_{i}表示變數x_{i},y_{i}的標準差。

1.2.1 利用馬氏距離進行資料歸一化

如圖1.2.1-1所示,當資料在空間中以非常不對稱的形式進行分佈時,k-means演算法總是試圖挖掘出一些與聚類相關的資訊,因為k-means聚類的核心觀點在於資料是以不均勻的方式進行聚類的。然而,“不對稱”和“不均勻”之間卻有著重要的區別。例如,當資料在某個維度上分佈很遠,而在其他維度上距離相對較小時,k-means必然不會收斂到好的結果。

         圖1.2.1-1 (a)原始資料的垂直距離比水平距離小 (b)對空間進行方差歸一化,資料之間的水平距離小於垂直距離

這種情況往往只是因為資料向量在不同的維度有不同的底層單位。例如,如果使用身高、年齡和學齡代表一個人,那麼由於單位的不同,資料在不同的維度上必然產生不同的分佈。同樣,年齡和學齡雖有著相同的單位,但是,在自然人口中也存在差異。

當我們將整個資料集作為一個整體來看待時,我們可以通過協方差矩陣來重新調整整個資料集,這種技術也稱之為資料歸一化

傳統意義上,馬氏距離是以該點在特定的方向上的方差來度量資料點與分佈之間的距離。我們使用分佈的協方差的倒數來計算馬氏距離:

\Sigma_{i,j} =E\left [ \left ( \mathbf{X}_{i}-\mu\right ) \left ( \mathbf{X}_{i}-\mu\right )^{T}\right ]=\frac{1}{N}\Sigma _{i}\left [ \left ( \mathbf{X}_{i}-\mu\right ) \left ( \mathbf{X}_{i}-\mu\right )^{T}\right ]

其中, E\left [ \cdot \right ] 表示數學期望,馬氏距離可表示為:

D_{mahananobis}(\mathbf{x},\mathbf{y})=\sqrt{\left ( \mathbf{x}-\mathbf{y}\right )^{T}\Sigma ^{-1}\left ( \mathbf{x}-\mathbf{y}\right )}

對此有以下兩種方式來解決:在K-means演算法中使用馬氏距離而不是歐式距離,或對資料進行尺度變化,然後在放縮的空間中使用歐幾里得距離。第一種方法更為直觀,但第二種方法更容易計算,因為資料轉化是線性的。

資料集的轉化D^{*}=D\Sigma ^{-1/2}

在這裡,D^{*}是即將使用的一組新的資料向量,D是原始資料。\Sigma ^{-1/2}的因子是逆協方差的平方根。

1.2.2 利用馬氏距離分類

利用K-means聚類或者任何其他的方法給定一組聚類標籤,利用這些標籤可以預測一些新的資料點最可能屬於哪個聚類。假設聚類本身服從高斯分佈,將馬氏距離的概念引進該分類問題也是非常有意義的。

分類方法:對每個聚類計算平均協方差,這樣我們可以基於協方差來計算新的資料點到每個聚類中心的馬氏距離,從而確定分類。

但是,是不是說具有最小馬氏距離的類是最優的解呢?

事實並非如此,當且僅當所有的聚類都擁有相同數量的資料點時(每個叢集內的資料點的先驗概率大小相等),此結論才成立。

貝葉斯規則簡明扼要的說了這一區別:

對於給定的命題A和命題B,給定B的情況下,A發生的概率一般來說並不等於給定A的情況下,B發生的概率

                                                                      P(A|B)\neq P(B|A))

相反,給定A的情況下B發生的概率,乘以A的概率必然等於給定B的情況下A發生的概率誠意B發生的概率

                                                              P(A|B)P(B)=P(B|A)P(A)

同樣的道理,馬氏距離表示的是一個特定樣本來自特定聚類的概率,事實上,我們想知道的是:給定變數x,它出現在某一特定類中的概率。顯然,馬氏距離告訴我們的結論恰恰是相反的。例如,對於給定的變數x出現在C類的概率:

                                                           P(C|\mathbf{x})\oe (\frac{N_{C}}{N_{D}})\left ( \left | \Sigma \right |^{-1/2} e^{-\frac{1}{2}r^{2}M}\right )

其中,N_{C},N_{D}表示聚類的資料量和整體資料裡。

也就是說,為了比較不同聚類的馬氏聚類,我們必須考慮聚類的大小。考慮到每個聚類的概率服從高斯分佈,聚類之間的似然性可以表示如下:

                                                                P(C|\mathbf{x})=\frac{P(C)}{P(\mathbf{x})}P(\mathbf{x}|C)

這就好比說,發現一隻看似恐龍的蜥蜴屬於恐龍的概率遠遠小於它屬於蜥蜴的概率,正是因為兩個叢集的先驗概率是不同的。

2. K均值的基本理論

2.1 K均值的基本原理和實現

K-means嘗試尋找資料的自然類別,使用者設定類別的個數,從而尋找到好的類別中心。

演算法的流程如下:

1.輸入資料集合和類別數K;

2.隨機分配類別中心點的位置;

3.將每個點放入離它最近的類別中心點所在的集合;

4.移動類別中心點到它所在的集合;

5.轉到第3步,直至收斂。

K均值的迭代示意圖如圖2.1-1所示:

                                                                圖2.1-1 K均值的迭代示意圖

2.2 K均值的缺點

K-means是一個極其高效的聚類演算法,但是它存在以下三個問題:

1.它不能保證定位到聚類中心的最佳方案,但能保證收斂到某個解決方案;

2.K-means無法指出應該使用多少個類別。在同一資料集中,選擇不同的類別得到的結果是不一樣的,甚至是不合理的;

3.K-means假定空間的協方差矩陣不會影響到結果。

2.3 K-means的改進

1.多進行幾次聚類,每次初始化的聚類中心不一樣,最後選擇方差最小的那個結果;

2.首先將類別設定為1,然後提高類別數(到達某個上限),一般情況下,總方差會快速下降直到達到某個拐點,這意味著再加一個新的聚類中心也不會顯著減少總體方差。在拐點處停止,儲存此時的類別數;

3.將資料乘以逆協方差矩陣進行資料的歸一化。

3.演算法實現

3.1獲取樣本

基本思路:利用隨機數RNG函式運算元對矩陣進行填充,這裡採用了正態分佈生成點集。

實現程式碼

        //定義實驗允許的最大類別數
	const int MAX_CLUSTERS = 5;
	//定義資料點的顏色
	Scalar colorTab[] = { Scalar(0,0,255),Scalar(0,255,0),Scalar(255,100,100),
                                                     Scalar(255,0,255),Scalar(0,255,255)};
	//基於均勻分佈生成隨機的K值和樣本數量
	RNG rng(12345);
	int k=rng.uniform(2, MAX_CLUSTERS);
	int sampleCount = rng.uniform(100,1000);
	Mat points(sampleCount, 1, CV_32FC2);
	Mat lables;
	Mat centers(k,1,CV_32FC1);
	Mat img(500, 500,CV_8UC3);
	//基於初始化的K值和sampleCount生成點集
	for (int i = 0; i < k; i++)
	{
		//指定ROI
		Mat pointChunk = points.rowRange(i*sampleCount/k,i==k-1?
                                                      sampleCount:(i+1)*sampleCount/k);
		//基於正態分佈填充ROI
		rng.fill(pointChunk,RNG::NORMAL,Scalar(rng.uniform(0,img.cols),
                         rng.uniform(0,img.rows)),Scalar(img.cols*0.05, img.rows*0.05));
	}
	//打亂點的生成順序
	randShuffle(points,1,&rng);

3.2 協方差逆陣平方根的計算方法

傳統意義上,馬氏距離是以該點在特定的方向上的方差來度量資料點與分佈之間的距離。我們使用分佈的協方差的倒數來計算馬氏距離:

\Sigma_{i,j} =E\left [ \left ( \mathbf{X}_{i}-\mu\right ) \left ( \mathbf{X}_{i}-\mu\right )^{T}\right ]=\frac{1}{N}\Sigma _{i}\left [ \left ( \mathbf{X}_{i}-\mu\right ) \left ( \mathbf{X}_{i}-\mu\right )^{T}\right ]

其中, E\left [ \cdot \right ] 表示數學期望,馬氏距離可表示為:

D_{mahananobis}(\mathbf{x},\mathbf{y})=\sqrt{\left ( \mathbf{x}-\mathbf{y}\right )^{T}\Sigma ^{-1}\left ( \mathbf{x}-\mathbf{y}\right )}

對此有以下兩種方式來解決:在K-means演算法中使用馬氏距離而不是歐式距離,或對資料進行尺度變化,然後在放縮的空間中使用歐幾里得距離。第一種方法更為直觀,但第二種方法更容易計算,因為資料轉化是線性的。

資料集的轉化D^{*}=D\Sigma ^{-1/2}

在這裡,D^{*}是即將使用的一組新的資料向量,D是原始資料。\Sigma ^{-1/2}的因子是逆協方差的平方根。

那麼如何求解\Sigma ^{-1/2}呢?

由於協方差矩陣是對稱矩陣,基於特殊矩陣而言採用特徵值分解的方法可以更準確的解出逆矩陣,現在我們假設\Sigma ^{-1}已經求出,且有:

                                                                   \Sigma ^{-1}=\begin{bmatrix} 0 & 1 & 1 &-1 \\ 1 & 0 &-1 &1 \\ 1&-1 &0 &1 \\ -1 & 1 & 1 &0 \end{bmatrix}

由線性代數的基本理論知道,計算矩陣的冪的基本思想是對矩陣進行對角化,即:

                                                                                 A^{n}=S\Lambda ^{n}S^{-1}

當n=2時,只需要對特徵向量計算平方,隨後反乘特徵向量S即可,計算程式碼如下:

	//矩陣對角化計算A的平方
	Mat A = (Mat_<float>(4,4)<<0,1,1,-1,1,0,-1,1,1,-1,0,1,-1,1,1,0);
	Mat value;
	Mat vector;
	bool M=eigen(A,value,vector);
	Mat eigenValueMat= Mat::zeros(Size(value.rows,value.rows),CV_32FC1);
        value.convertTo(value,CV_32FC1);
        vector.convertTo(vector,CV_32FC1);
	for (int i = 0; i < value.rows; i++)
	{
		eigenValueMat.at<float>(i, i) = pow(value.at<float>(i),2);
	}
	Mat eigenVectorMat;
	transpose(vector, eigenVectorMat);
	//對角化方法計算的解
	Mat A_1 = eigenVectorMat*eigenValueMat* eigenVectorMat.inv();
	//理論解
	Mat A_2 = A*A;
	//誤差
	Mat error = A_1 - A_2;

最終的誤差矩陣error如圖3.2-1所示:

                                                                                         圖3.2-1 誤差矩陣

顯然,這個解是完全可以接受的

3.3 聚類實驗

3.3.1 一般的K均值聚類方法

在該條件下,對資料的樣本點不加任何歸一化處理,並基於已知的K值來對3.1中的樣本進行聚類分析,實現的程式碼如下:

	//1.一般的均值聚類
	//1.1 直接呼叫聚類函式
	kmeans(points,k,lables,TermCriteria(TermCriteria::EPS | TermCriteria::EPS,10,1.0),3,KMEANS_RANDOM_CENTERS,centers);
	//1.2 繪製聚類結果
	for (int i = 0; i < sampleCount; i++)
	{
		int colorIndex = lables.at<int>(i);
		Point pt = points.at<Point2f>(i);
		circle(img,pt,2,colorTab[colorIndex],-1);
	}

聚類的結果如下圖3.3.1-1所示:

                                                                                 圖3.3.1-1 一般的均值聚類

3.3.2 基於馬氏距離k-means++的聚類方法的實現

首先,為了創造馬氏距離聚類的基本條件,在這裡為了更好的視覺化,資料採用二維向量,但是在x和y兩個維度上的服從不同的分佈,為了使資料更具可信度,我參考了《機器學習》(周志華著)原型聚類章節裡面的樣本資料,如表3.3.2-1:

                                                                        表3.3.2-1 西瓜資料集

其次,在原型聚類的例項中,書中給出的最終的迭代解,如圖3.3.2-1所示:

                                                              圖3.3.2-1 西瓜資料集均值聚類的結果

因此,為了方便進行資料實驗,不妨設定聚類數目為k=3,實驗的程式碼如下:

        //2.基於馬氏距離的k-means++聚類
	//定義實驗允許的最大類別數
	const int MAX_CLUSTERS = 5;
	//定義資料點的顏色
	Scalar colorTab[] = { Scalar(0,0,255),Scalar(0,255,0),Scalar(255,100,100),Scalar(255,0,255),Scalar(0,255,255) };
	//設定基本引數,k=3
	RNG rng(12345);
	int k = 3;
	Mat lables;
	Mat centers(k, 1, CV_32FC1);
	Mat img(500, 500, CV_8UC3, Scalar(0));
	Mat re_points = (Mat_<float>(30, 2) << 0.697,0.460,0.774,0.376,0.634,0.264,0.608,0.318,0.556,
		0.215,0.403,0.237,0.481,0.149,0.437,0.211,0.666,0.091,0.243,0.267,0.245,0.057,0.343,0.099,
		0.639,0.161,0.657,0.198,0.360,0.370,0.593,0.042,0.719,0.103,0.359,0.188,0.339,0.241,0.282,
		0.257,0.748,0.232,0.714,0.346,0.483,0.312,0.478,0.437,0.525,0.369,0.751,0.489,0.532,0.472,
		0.473,0.376,0.725,0.445,0.446,0.459);
	//求解協方差矩陣
	Mat covar;
	Mat mean;
	calcCovarMatrix(re_points, covar,mean,CV_COVAR_NORMAL | CV_COVAR_ROWS);
	//求解協方差的逆陣的方根
	Mat covar_inv;
	double inv=invert(covar,covar_inv,DECOMP_EIG);
	Mat value;
	Mat vector;
	bool M = eigen(covar_inv, value, vector);
	Mat eigenValueMat = Mat::zeros(Size(value.rows, value.rows), CV_32FC1);
	value.convertTo(value,CV_32FC1);
	vector.convertTo(vector, CV_32FC1);
	for (int i = 0; i < value.rows; i++)
	{
		float tem = sqrt(value.at<float>(i, 0));
		eigenValueMat.at<float>(i, i) = tem;
	}
	Mat eigenVectorMat;
	transpose(vector, eigenVectorMat);
	Mat result = eigenVectorMat*eigenValueMat* eigenVectorMat.inv();
	//基於馬氏距離進行歸一化,重新生成點集
	Mat New_points = re_points*result;
	//聚類
	kmeans(New_points, k, lables, TermCriteria(TermCriteria::EPS | 
                           TermCriteria::EPS, 10, 1.0), 3, KMEANS_PP_CENTERS, centers);
	for (int i = 0; i < 30; i++)
	{
		int colorIndex = lables.at<int>(i);
		int x = re_points.at<float>(i, 0) * 500;
		int y= 500-re_points.at<float>(i, 1) * 500;
		Point pt = Point(x,y);
		circle(img, pt, 2, colorTab[colorIndex], -1);
	}

最後,我想說的是k-means++與k-means方法最大的不同點在於中心的初始化方法,因為對於不同的初始點,聚類的結果必然會有差異,如果初始點選擇的不好,那麼極有可能陷入區域性最小值下,從而得不到好的分類結果,而cv::KMEANS_PP_CENTERS選項會使cv::kmeans使用文獻"Arthur, David, and S. Vassilvitskii. "k-means++: the advantages of careful seeding"中所謂的kmeans++的方法來重新分配聚類中心。這種方法謹慎的選擇了聚類的中心起點,通常能以少於預設方法的迭代得出更好的結果。

上述程式碼的結果如圖3.3.2-2所示:

                                                                圖3.3.2-2 基於馬氏距離的kmeans++聚類結果

3.3.3 基於肘部法則的K-means++聚類

        如果問題中沒有指定的值,可以通過肘部法則這一技術來估計聚類數量。肘部法則會把不同值的成本函式值畫出來。隨著值的增大,平均畸變程度會減小;每個類包含的樣本數會減少,於是樣本離其重心會更近。但是,隨著值繼續增大,平均畸變程度的改善效果會不斷減低。值增大過程中,畸變程度的改善效果下降幅度最大的位置對應的值就是肘部。

在這裡,我將畸變程度定義為所有類的樣本點距該類中心的方差和

考慮k從k=2到k=5逐步進行迭代,求解產生的畸變係數的程式如下:

        //肘部法則的聚類
	//定義實驗允許的最大類別數
	const int MAX_CLUSTERS = 5;
	//定義資料點的顏色
	Scalar colorTab[] = { Scalar(0,0,255),Scalar(0,255,0),Scalar(255,100,100),
                                   Scalar(255,0,255),Scalar(0,255,255) };
	//設定基本引數
	RNG rng(12345);
	int k;
	Mat lables;
	Mat img(500, 500, CV_8UC3, Scalar(0));
	Mat points = (Mat_<float>(30, 2) << 0.697, 0.460, 0.774, 0.376,
                     0.634, 0.264,  0.608, 0.318, 0.556,
		0.215, 0.403, 0.237, 0.481, 0.149, 0.437, 0.211, 0.666, 
                     0.091, 0.243, 0.267, 0.245, 0.057, 0.343, 0.099,
		0.639, 0.161, 0.657, 0.198, 0.360, 0.370, 0.593, 0.042,
                     0.719, 0.103, 0.359, 0.188, 0.339, 0.241, 0.282,
		0.257, 0.748, 0.232, 0.714, 0.346, 0.483, 0.312, 0.478,
                     0.437, 0.525, 0.369, 0.751, 0.489, 0.532, 0.472,
		0.473, 0.376, 0.725, 0.445, 0.446, 0.459);
	//定義方差
	Mat varMat(5, 1, CV_32FC1, Scalar(0));
	for (int k = 2; k <= MAX_CLUSTERS; k++)
	{
		Mat centers(k, 1, CV_32FC1);
		kmeans(points, k, lables, TermCriteria(TermCriteria::EPS | 
        TermCriteria::EPS, 10, 1.0), 3, KMEANS_PP_CENTERS, centers);
		vector<float> sumVar(k, 0);
		//統計每類點的數量
		Mat lableCount(k, 1, lables.type(), Scalar(0));
		for (int i = 0; i < 30; i++)
		{
			lableCount.at<int>(lables.at<int>(i)) += 1;
		}

		//計算方差
		for (int i = 0; i < 30; i++)
		{
			int index = lables.at<int>(i);
			float x = points.at<float>(i, 0);
			float y = points.at<float>(i, 1);
			Point2f center = Point2f(centers.at<float>(index, 0), 
      centers.at<float>(index, 1));
			sumVar[index] += ((x - center.x)*(x - center.x) + 
      (y - center.y)*(y - center.y)) / lableCount.at<int>(index);

		}
		for (auto m : sumVar)
		{
			varMat.at<float>(k - 2, 0) += m;
		}
	}
        //資料歸一化[0,500]
	normalize(varMat,varMat,0,400,NORM_MINMAX);
	//繪圖
	vector<Point> pts;
	for (int i = 0; i < varMat.rows; i++)
	{
		Point tem;
		tem.x = i*400/4+50;
		tem.y = 450 - varMat.at<float>(i);
		pts.push_back(tem);
	}
	for (int i = 0; i < pts.size()-1;i++)
	{
		circle(img,pts[i],3,colorTab[i],-1);
	}
	line(img,pts[0],pts[1],Scalar(255,255,255));
	line(img, pts[1], pts[2], Scalar(255, 255, 255));
	line(img, pts[2], pts[3], Scalar(255, 255, 255));

最終的結果如圖3.3.3-1所示:

                                                                圖3.3.3-1 基於肘部法則確定拐點

很顯然,拐點發生在k=3的地方,因此最佳的聚類結果為k=3

4.相關資料

1.OpenCV隨機數發生器RNG:

《學習OpenCV3》:Page 157-161

2.對小矩陣使用逗號分隔式初始化函式:

《OpenCV程式設計入門》:Page 91

3.矩陣對角化理論、矩陣對角化Matlab實現:

https://blog.csdn.net/qq_18343569/article/details/49823441

https://blog.csdn.net/compression/article/details/49180775

4.聚類演算法的一些思考:

https://blog.csdn.net/u013719780/article/details/51755124?utm_source=blogxgwz2

https://www.jianshu.com/p/95a4bcff2198

https://www.cnblogs.com/data-miner/p/6288227.html

5.原型聚類的基本理論:

《機器學習》-周志華:Page 202-211

6.維基百科:歐式距離與馬氏距離

https://en.wikipedia.org/wiki/Euclidean_distance

https://en.wikipedia.org/wiki/Mahalanobis_distance