論文研讀-ClusteringGA研讀與實現

此為課題組2024級研究生何諾飄同學近期學習內容匯報
更多內容請關注許志偉課題組官方中文主頁：https://JaywayXu.github.io/zh-cn/

1. 前言

ClusteringGA算法的整體流程與NSGAII相似，主要包括初始種群的產生、聚類、交叉、變異、非支配排序、選擇等。其中聚類這個操作是ClusteringGA的主要創新點，而交叉變異也圍繞聚類而進行，故本文將主要圍繞聚類操作進行解讀。

2. 聚類(Cluster)的介紹

聚類是一種機器學習技術，它涉及到數據點的分組。給定一組數據點，我們可以使用聚類算法將每個數據點劃分為一個特定的組。理論上，同一聚類中的數據點應該具有相似的屬性和/或特征，而不同聚類中的數據點應該具有高度不同的屬性和/或特征。聚類是一種無監督學習的方法，是許多領域中常用的統計數據分析技術。

K-means聚類：最常用的聚類算法，通過指定K個中心點，將數據點分配到最近的中心點所在的簇。
AP聚類（Affinity Propagation）：通過在數據點之間傳遞消息來識別聚類中心，不需要預先指定聚類數量，能自動確定最優的聚類數量。
層次聚類（Hierarchical Clustering）：通過構建樹狀層次結構來進行聚類，可以自底向上（凝聚）或自頂向下（分裂）進行。
DBSCAN（基于密度的空間聚類）：基于密度的聚類方法，可以發現任意形狀的簇，對噪聲數據具有較強的魯棒性。
期望最大化（EM）聚類：基于概率模型的聚類方法，通過迭代優化來估計模型參數。
譜聚類（Spectral Clustering）：利用數據點之間的相似度矩陣的特征向量來降維，然后在低維空間進行聚類。

而在ClusteringGA算法中主要采用的則是AP聚類。

3. AP聚類

3.1 AP聚類的介紹

3.1.1 AP聚類的簡介

AP聚類（Affinity Propagation）是一種基于消息傳遞的聚類算法，它將數據點視為網絡中的節點，通過節點之間的信息交換來確定聚類中心和分類。算法的核心思想是讓所有數據點之間通過傳遞消息來"協商"，從而選擇出最佳的聚類中心（exemplar）。

在AP聚類中，每個數據點都可能成為聚類中心。算法通過迭代計算兩種消息：

責任度（Responsibility）矩陣：表示數據點k作為數據點i的聚類中心的適合程度
可用度（Availability）矩陣：表示數據點i選擇數據點k作為其聚類中心的適合程度

通過這兩種消息的不斷傳遞和更新，算法最終會收斂到一組穩定的聚類中心，從而完成數據的聚類過程。每個數據點都會被分配到與其最相似的聚類中心所在的類別中。

3.1.2 AP聚類的優點

AP聚類（Affinity Propagation）算法具有以下幾個主要優點：

不需要預先指定聚類數量：算法能夠自動確定最優的聚類數量，這在實際應用中非常有價值，特別是當我們不清楚數據應該分為多少類時
適用性廣：可以處理相似度矩陣不對稱的情況，也可以處理非歐氏距離的相似度度量
聚類效果好：通過消息傳遞機制，能夠找到具有代表性的樣本點作為聚類中心，聚類結果更加合理
收斂性好：算法具有良好的收斂性，且收斂速度相對較快
并行計算：算法的消息傳遞過程可以并行化，提高計算效率

3.2 AP聚類的流程

3.2.1 聚類前期準備工作

在進行AP聚類之前，我們需要進行一些準備操作。在這個基準問題中，決策變量是一條條路徑，而AP聚類是基于適應度值進行操作的，因此我們需要對決策變量進行轉換。在種群初始化時，我們根據預定的種群大小生成整個種群，并用數組列表表示，其中列表的每一行代表一個個體。在AP聚類中，我們使用種群中個體所在的行號作為該個體的標識。在計算相似度時，需要將每個個體的路徑長度統一；對于長度不等的路徑，較短的路徑將在末尾補0以對齊長度。相似度值通過計算路徑中對應位置的路徑點之間的負歐式距離來表示。

下方圖片中表示的則是種群以及種群中的部分個體

3.2.2 相似度/參考度矩陣(Similarity Matrix)

AP聚類的第一步則是生成相似度矩陣（以下將簡稱 s 矩陣），我們通常記為 s(i,k)，是指點k作為點i的聚類中心的相似度，數據點 i 和點 k 的相似度記為 s(i,k)，從點 i 發送至候選聚類中心點 k,反映了在考慮其他潛在聚類中心后，點 k 適合作為點i的聚類中心的程度。這里采用歐氏距離來計算，將點與點的相似度值全部取為負值，因此，s(i,k) 相似度值越大說明點 i **與點 k 的距離越近，AP算法中理解為數據點 k 作為數據點i的聚類中心的能力。

為了更加方便理解，以下提供一個簡單的例子：假設點分布在實數軸上，坐標分別為：

A = 1，B = 2，C = 3，D = 5，E = 6

A B C D E

0—1—2—3—4—5—6—7—8—9—>

用兩個點之間的距離的負數作為兩個點之間的相似度，則s矩陣生成為：

s(i/k)	A	B	C	D	E
A	?	-1.0	-2.0	-4.0	-5.0
B	-1.0	?	-1.0	-3.0	-4.0
C	-2.0	-1.0	?	-2.0	-3.0
D	-4.0	-3.0	-2.0	?	-1.0
E	-5.0	-4.0	-3.0	-1.0	?

對角線 s(i,i) 表示自己與自己的相似度，AP聚類算法通常選擇上述矩陣中元素的最小值或者中位數來填充對角線元素。接下來我們選擇中位數得到完整的s矩陣：

s(i/k)	A	B	C	D	E
A	-3.0	-1.0	-2.0	-4.0	-5.0
B	-1.0	-3.0	-1.0	-3.0	-4.0
C	-2.0	-1.0	-3.0	-2.0	-3.0
D	-4.0	-3.0	-2.0	-3.0	-1.0
E	-5.0	-4.0	-3.0	-1.0	-3.0

以下是源碼中 s 矩陣的實現片段：

// 生成相似度矩陣Similarity
double[][] s = InitialArray(clustering,paths,maxLength); // 生成s矩陣：每個個體與其他個體的相似度（歐式距離）
double[][] S = new double[paths.size()][paths.size()]; // 填充s的對角線元素生成完整S矩陣
int l = 0;
for(int i=0;i<paths.size();i++){
    for(int j=0;j<paths.size();j++){
        if(i==j){S[i][j]=p;} // 對角線存相似度中位數
        else {S[i][j] = s[l][2];l++;} // 除對角線外元素存s矩陣第三列相似度值
    }
}
        
 // 創建初始化s數組
private static double[][] InitialArray(int[][] clustering, List<Path> paths, int maxLength) {
	int N = paths.size();
	int M = N * (N - 1); // 排除對角線
	double[][] s = new double[M][3]; // 創建一個 M 行 3 列的矩陣，用于存放個體之間的相似度
	int j = 0;
	
	for (int i = 0; i < N; i++) {
	    for (int k = 0; k < N; k++) {
	        if(i == k) continue; // 跳過自己
	        double sum = 0; // 每對 (i, k) 重新計算 sum
	        for (int l = 0; l < maxLength; l++) {
	            double first = clustering[i][l];
	            double second = clustering[k][l];
	            if(first == 0 && second == 0) {continue;} // 若x、y都為0，則沒必要進行操作直接跳出
	            // 將路徑點拆分為x、y坐標方式進行之后的歐式距離計算
	            double first_x = (double) first / 1000;
	            double first_y = (double) first % 1000;
	            double second_x = (double) second / 1000;
	            double second_y = (double) second % 1000;
	            // 計算歐式距離
	            double aum = Math.sqrt(Math.pow(first_x - second_x, 2) + Math.pow(first_y - second_y, 2));
	            sum += aum * aum; // 歐式距離平方累加
	        }
	        s[j][0] = i+1;
	        s[j][1] = k+1;
	        s[j][2] = -sum;// 負歐式距離
	        j++;
	    }
	}
return s;
}

3.2.3 責任度/吸引度矩陣(Responsibility Matrix)

在完成好相似度矩陣后，便進入了主循環的迭代，所以AP聚類的第二步便是在主循環迭代中更新責任度矩陣（以下將簡稱 r 矩陣），r 用來描述點 k 適合作為數據點 i 的聚類中心的程度，上面的 s 矩陣雖然在一定程度上反映了程度，但是，不同行之間數據是不能進行比較的。例如，第1行第2列最大值是-1，意味著個體A會選個體B作為聚類中心。但是第2行的最大值-1有兩個，這樣不利于對聚類中心的選擇，所以我們需要再創建一個r矩陣，來將這些數據就行優化利于聚類中心的選擇。

r 矩陣按照以下數學公式來進行生成（在初始化時，a(i,k’) = 0）：

\[r(i,k) = s(i,k) - \max_{k' \neq k}\{s(i,k') + a(i,k')\} \]

讓我們詳細解讀公式 r(i,k) = s(i,k) - max_{k' ≠ k}{s(i,k') + a(i,k')}：

r(i,k)：表示點i選擇點k作為其聚類中心的程度（責任度）
s(i,k)：表示點i與點k之間的相似度
max_{k' ≠ k}{s(i,k') + a(i,k')}：表示除了k以外，其他所有點作為i的聚類中心的競爭力最大值
- k'表示除k以外的其他點
- s(i,k')表示i與其他點k'的相似度
- a(i,k')表示i選擇k'作為聚類中心的適合程度（在初始時為0）

該公式的計算過程：將 s 矩陣中每一個元素與當前這一行中最大的元素做差。

這個公式的本質是：通過將點i和k的相似度，減去i與其他所有候選中心點的最大競爭力，來確定k作為i的聚類中心的程度。如果結果為正值，說明k比其他所有候選點更適合作為i的聚類中心。

根據以上描述，舉例表格中的元素更新為：

r(i/k)	A	B	C	D	E
A	-2.0	1.0	-1.0	-3.0	-4.0
B	0.0	-2.0	0.0	-2.0	-3.0
C	-1.0	0.0	-2.0	-1.0	-2.0
D	-3.0	-2.0	-1.0	-2.0	1.0
E	-4.0	-3.0	-2.0	2.0	-2.0

為防止 r 矩陣的更新幅度過大造成結果不收斂后我們會生成 rold 矩陣來記錄剛更新的 r 矩陣，并引入阻尼系數 λ，來控制更新幅度。一般取 λ = 0.5 ，更新公式如下：

\[r = (1-\lambda)r + \lambda \cdot r_{old} \]

對于舉例表格而言，由于沒有 rold，故我們直接將表格中的數據乘以0.5:

r(i/k)	A	B	C	D	E
A	-1.0	0.5	-0.5	-1.5	-2.0
B	0.0	-1.0	0.0	-1.0	-1.5
C	-0.5	0.0	-1.0	-0.5	-1.0
D	-1.5	-1.0	-0.5	-1.0	0.5
E	-2.0	-1.5	-1.0	1.0	-1.0

以下是源碼中 r 矩陣的實現片段：

// 更新R責任度矩陣（將S中的每一個元素與該元素所處行最大值相減）
double[][] Xpro = new double[paths.size()][paths.size()]; // 將每行最大值填充到一個新矩陣
for (int i = 0; i < paths.size(); i++) {
    for (int j = 0; j < paths.size(); j++) {
        Xpro[i][j] = Y[i];
    }
}
R = calMetrices(S,Xpro,2); // 按公式計算R矩陣（calMertices實現了矩陣加減操作，后面數字2表示是減法操作）
// 采取手段使本來優秀的個體更加優秀
for(int k = 0;k<paths.size();k++){
    int num = I[k]; // 用來記錄S中每一行最大值的位置（列索引）
    R[k][num] = S[k][num] - Y2[k]; // 將R中原來每一行最大值替換為最大值 - 第二大值
}

以下是源碼中阻尼系數實現片段：

R = addLam(R,Rold,0.5); // 加入阻尼系數(R=(1-lam)*R+lam*Rold)
//加入阻尼系數
private static double[][] addLam(double[][] metrix1,double[][] matrix2,double lam) {
    int rows = metrix1.length;
    int cols = metrix1[0].length;
    double[][] result = new double[rows][cols]; // 創建存儲結果的矩陣
    double[][] temp1 = new double[rows][cols]; // 新矩陣（更新后的矩陣）
    double[][] temp2 = new double[rows][cols]; // 舊矩陣（更新前的矩陣）
    for (int i = 0; i < metrix1.length; i++) { // 遍歷行
        for (int j = 0; j < metrix1[i].length; j++) { // 遍歷列
            temp1[i][j] = metrix1[i][j] * (1-lam); // 新矩陣元素乘以系數
        }
    }
    for (int i = 0; i < matrix2.length; i++) {
        for (int j = 0; j < matrix2[i].length; j++) {
            temp2[i][j] = matrix2[i][j] * lam; // 舊矩陣元素乘以系數
        }
    }
    result = calMetrices(temp1,temp2,1); // 將乘以系數后的兩個新矩陣相加生成最終結果矩陣
    return result;
}

3.2.4 可用度/歸屬度矩陣(Availability Matrix)

AP聚類的第三步是在主循環迭代中更新可用度矩陣（以下將簡稱 a 矩陣）a(i,k) 用來描述點 i 選擇點 k 作為其聚類中心的適合程度。從候選聚類中心點 k 發送至點 i，反映了在考慮其他點對點 k 成為聚類中心的支持后，點 i 選擇點 k 作為聚類中心的合適程度，從候選聚類中心點 k 發送至點 i，反映了在考慮其他點對點 k 成為聚類中心的支持后，點 i 選擇點 k 作為聚類中心的合適程度。簡單來說就是根據目前 r 矩陣提供的支持度結果，做進一步的決策調整。

a 矩陣按照以下數學公式來進行生成：

\[a(i,k) = \min\{0, r(k,k) + \sum_{i' \neq i,k}\max\{0,r(i',k)\}\} \quad (i \neq k) \]

\[a(k,k) = \sum_{i' \neq k}\max\{0,r(i',k)\} \]

其中，第一個公式用于計算非對角線元素，第二個公式用于計算對角線元素。

讓我們來詳細解讀這兩個計算公式：

\(a(i,k) = \min\{0, r(k,k) + \sum_{i' \neq i,k}\max\{0,r(i',k)\}\} \quad (i \neq k)\)

r(k,k)：表示候選聚類中心 k 對自身的責任度，即 k 作為聚類中心的可能性
∑max{0,r(i',k)}：對除了 i 和 k 之外的所有點 i'，計算它們對 k 作為聚類中心的正向支持度的總和
min{0, ...}：將結果限制在負值范圍內，這樣可以防止可用度值過大

該公式的計算過程如下：

首先計算r(k,k)值，即 k 點對自身的責任度
對當前元素所在這一列中所有除了i和k之外的點i'，找出它們對 k 的責任度r(i',k)中的正值并求和
將r(k,k)與上一步得到的和相加
最后取這個結果與0比的較小值，確保結果不會是正數

這個公式的核心思想是：當計算點 i 選擇 k 作為聚類中心的適合程度時，不僅要考慮 k 自身作為中心的傾向，還要考慮其他點對 k 的支持程度。如果其他點都強烈支持 k 作為中心（即有較多正的 r 值），那么 i 選擇 k 作為中心也是較為合適的。本質上是在評估：如果其他點都認為 k 適合作為聚類中心（正的 r 值），那么 i 選擇 k 作為其聚類中心也是合理的。

\(a(k,k) = \sum_{i' \neq k}\max\{0,r(i',k)\}\)

∑：表示需要對所有除 k 以外的點進行累加運算
i' ≠ k：表示在累加過程中，排除當前考慮的點 k 自身
max{0,r(i',k)}：只考慮正的責任度值，如果責任度為負，則取0。

該公式的計算過程如下：

遍歷責任度矩陣 r 中第 k 列的所有元素（除了第 k 行）
對每個元素，如果值為正數則保留，如果為負數則記為0
將所有處理后的值相加，得到的結果就是a(k,k)

這個公式的核心思想是：它直接反映了其他點對 k 作為聚類中心的集體支持度，通過只考慮正值，消除了反對票的影響，使得聚類中心的選擇更加穩健，對角線元素的值越大，表示該點越有可能成為最終的聚類中心。

根據非對角線元素更新公式的描述，舉例表格中的元素更新為：

a(i/k)	A	B	C	D	E
A	-1.0	-1.0	-1.0	0.0	-0.5
B	-1.0	-1.0	-1.0	0.0	-0.5
C	-0.5	-0.5	-1.0	0.0	-0.5
D	-0.5	-0.5	-1.0	-1.0	-1.0
E	-0.5	-0.5	-1.0	-1.0	-1.0

進一步再根據對角線元素更新公式，舉例表格中的元素更新為：

a(i/k)	A	B	C	D	E
A	0.0	-1.0	-1.0	0.0	-0.5
B	-1.0	0.5	-1.0	0.0	-0.5
C	-0.5	-0.5	0.0	0.0	-0.5
D	-0.5	-0.5	-1.0	1.0	-1.0
E	-0.5	-0.5	-1.0	-1.0	0.5

最后還是為防止 a 矩陣的更新幅度過大造成結果不收斂，生成 aold 矩陣來記錄剛更新的 a 矩陣，并引入阻尼系數 λ，來控制更新幅度。取 λ = 0.5 ，舉例表格中的元素更新為(aold取0)：

a(i/k)	A	B	C	D	E
A	0.0	-0.5	-0.5	0.0	-0.25
B	-0.5	0.25	-0.5	0.0	-0.25
C	-0.25	-0.25	0.0	0.0	-0.25
D	-0.25	-0.25	-0.5	0.5	-0.5
E	-0.25	-0.25	-0.5	-0.5	0.25

以下是源碼中 a 矩陣的實現片段：

// 更新A可用度矩陣
Aold = A; // 將更新前的A提前記錄下來(為后續加入阻尼系數提供方便）
Rp = R; // Rp作為過渡矩陣
for(int i=0;i<paths.size();i++){
    for(int j=0;j<paths.size();j++){
        if(Rp[i][j] < 0 && i != j){
            Rp[i][j] = 0; // 除R(k,k)外，將R中的負數變為0，忽略不適合的點的不適合程度信息
        }
    }
}
for(int i=0;i<paths.size();i++){
    Rp[i][i] = R[i][i]; // Rp中對角線的值替換為之前R中的對角線的值
}
double[][] Ypro_add = new double[1][paths.size()]; // 按列求和
double[][] Ypro = new double[paths.size()][paths.size()]; // 將求和結果擴充成一個矩陣
for (int j = 0; j < paths.size(); j++) {
    double a = 0;
    for (int i = 0; i < paths.size(); i++) {
        a += Rp[i][j];
        Ypro_add[0][j] = a; // 求和
    }
}
for (int i = 0; i < paths.size(); i++) {
    for (int j = 0; j < paths.size(); j++) {
        Ypro[i][j] = Ypro_add[0][j]; // 擴充
    }
}
A = calMetrices(Ypro,Rp,2); // A=repmat(sum(Rp,1),[N,1])-Rp;按公式計算A矩陣
for(int i=0;i<paths.size();i++){
    dA[i][i] = A[i][i]; // 將A矩陣對角線元素保存
}
for(int i=0;i<paths.size();i++){
    for(int j=0;j<paths.size();j++){
        if(i != j && A[i][j] > 0){ A[i][j] = 0;} // 除A(k,k)外其他大于0的值置為0
    }
}
for(int i=0;i<paths.size();i++){
    A[i][i] = dA[i][i]; // 將之前保留的A對角線元素換到新的A中
}
A = addLam(A,Aold,0.5); // 加入阻尼系數A=(1-lam)*A+lam*Aold;
for(int i=0;i<paths.size();i++){
    dA[i][i] = A[i][i]; // 將A矩陣對角線元素保存
}

3.2.5 最終判斷

在完成了上述 s 矩陣、r 矩陣與a 矩陣的生成后，我們設定 E 為一個邏輯向量，記錄了哪些點被判定為簇中心也就是聚類中心，判定規則是 r(i,i)與a(i,i)依次相加(也就是 r 矩陣與 a 矩陣對角線元素之和)，若相加之和大于0，則令E(i) = 1；記錄完后，將每一次的循環結果，依次循環存入矩陣 e 中，即第一次循環得出的 E 放到 e 矩陣的第一列，第51次的結果又放在第一列。再然后將 e 中元素按行求和放進 se 中，se 是一個列向量，用來記錄每個點在過去 convits 次迭代中的簇中心狀態之和。如果一個點在所有 convits 次迭代中都為 1 或都為 0，說明這個點已經收斂。通過sum 統計已經收斂的點的數量，如果所有點都收斂（即 sum == 初始生成種群大小），則說明算法整體已經收斂。最終如果所有點收斂且 K > 0（存在簇中心），或者迭代次數達到最大值 maxits，算法結束。很顯然只有在所有點收斂時表示聚類成功，而如果僅僅只是迭代次數達到最大值則代表著聚類失敗。

上述操作反應到舉例表格中則是：

s(i/i)+a(i/i)	A	B	C	D	E
A	-1.0
B		-0.75
C			-1.0
D				-0.5
E					-0.75

很顯然上述舉例表格中的對角線元素都小于0，那么他們都無法成為簇1中心點，這時則說明沒有聚類中心，那么將繼續下輪的迭代。上述的迭代計算結果舉例如下，供大家參考：

s(i/i)+a(i/i)	A	B	C	D	E
A	-2.75	0.5	-2.375	-2.75	-4.5
B	-1.625	-0.75	-1.625	-1.5	-3.25
C	-2.375	0.25	-2.75	-0.75	-2.5
D	-4.125	-1.25	-2.125	-1.0	-0.5
E	-5.125	-2.25	-3.125	0.5	-1.75

s(i/i)+a(i/i)	A	B	C	D	E
A	-1.71875	1.46875	-1.28125	-1.875	-3.4375
B	-0.9375	0.96875	-0.75	-0.1875	-1.75
C	1.59375	1.0625	-1.59375	0.0	-1.5625
D	-2.84375	-0.25	-0.90625	0.375	0.0
E	-3.84375	-1.25	-1.65625	0.875	-0.4375

可以看到在最后一張表中，已經選出了最終的聚類中心。

以下是源碼中最終判斷的實現片段：

int[] E = new int[paths.size()]; // E是一個邏輯向量，記錄了哪些點被判定為簇中心
for(int i=0;i<paths.size();i++){
    if(dA[i][i] + dR[i][i] > 0){
        E[i] = 1;
    }
}
// 將循環計算結果列向量E放入矩陣e中，注意是循環存放結果
// 即第一次循環得出的E放到N*50的e矩陣的第一列，第51次的結果又放到第一列
double[][] e = new double[paths.size()][convits];
// 計算列索引
int colIndex = (counter - 1) % convits; // Java 的 % 操作符等效于 MATLAB 的 mod
if (colIndex < 0) colIndex += convits; // 處理負數情況，確保結果在 [0, convits-1]
colIndex++; // MATLAB 的索引從 1 開始，Java 的索引從 0 開始，所以加 1
// 將 E 的內容賦值到 e 的第 colIndex 列
for (int row = 0; row < E.length; row++) {
    e[row][colIndex - 1] = E[row]; // 注意 Java 的列索引從 0 開始
}

int K = 0;
for (int i = 0; i < paths.size(); i++) {
    if(E[i] == 1){K++;} // 計算 E 中值為 1 的數量，即當前被判定為簇中心的點的個數K
}

// 最終判斷
if (counter >= convits || counter >= maxits) {
    // 將e中元素按行求和放進se中
    double[] se = new double[paths.size()]; // se是一個列向量，E的convits次迭代和
    for (int i = 0; i < paths.size(); i++) {
        double sumse = 0;
        for (int j = 0; j < convits; j++) {
            sumse += e[i][j];
        }
        se[i] = sumse;
    }
    boolean unconverged; // 用來判斷收斂點的數量是否等于總數
    double sum = 0;
    for (int i = 0; i < paths.size(); i++) {
        if(se[i] == convits || se[i] == 0){
            sum++;
        }
    }
    unconverged = (sum != paths.size()); // 如果收斂的點數量不等于總點數N,說明有點還未收斂,unconverged 為真
    // 收斂
    if ((!unconverged && (K > 0)) || (counter == maxits)) {
        dn = true; // 如果條件成立，將 dn 設置為 true
    }
}

3.2.6 聚類后續完善

經過以上的整個操作，我們選出了聚類的中心點，最后我們需要按照我們所計算出的簇中心點，來將簇中心點以及每個所屬的子個體放在同一個簇中，來方便后續的交叉變異操作，并且也需要將每一個個體還原成原來的路徑表示。

該階段的大致流程為：

找出潛在的簇中心點
分配點到最近的簇中心
細化簇中心點
計算聚類適合度
去重和排序簇中心點
構建簇并分配路徑數據
輸出聚類結果

以下是源碼中后續完善的實現片段：

// 經過上面的的循環，便確定好了哪些點可以作為簇中心點，找出簇1中心點
// int[] I = new int[paths.size()];
ArrayList<Integer> I = new ArrayList<>();
int K = 0; // 用來記錄I中存放多少數據，及多少個簇中心點
int sum_k = 0; // 用來方便在I中按行記錄簇中心點索引
for (int i = 0; i < paths.size(); i++) {
    // 將R與A矩陣對角線元素相加，若為正數則說明該點為簇中心點，且i反應了簇中心點的行列（位置）
    if(R[i][i] + A[i][i] > 0){
        //I[sum_k] = i;
        I.add(i);
        sum_k++;
        K++;
    }
}

double[] tmpidx = new double[paths.size()]; // 用來存放最后每個點的聚類中心
double tmpnetsim = 0; // 存放各個點到聚類中心的負歐式距離的和（衡量這次聚類的適合度）
double tmpexpref = 0; // 存放所有被選為簇中心的點的適合度

if(K > 0){
    // [~,c]=max(S(:,I),[],2);
    // 將S中對應I的所有列取出來，例如I中是3，4，則將S的第3和第4列的所有元素取出來
    // （判斷除此簇中心點外，其他的點是否以此簇中心點為中心）
    // 對取出的這幾列，按行查找最大值并記錄這個最大值的列索引
    double[][] Sp = new double[paths.size()][paths.size()]; // 用來存S中對應I的列
    // 將Sp內部填充負無窮大，方便后續查找最大值的列索引
    for (int i = 0; i < paths.size(); i++) {
        for (int j = 0; j < paths.size(); j++) {
            Sp[i][j] = Double.NEGATIVE_INFINITY;
        }
    }
    int[] c = new int[paths.size()]; // 存放列索引（存放的是每個點的簇中心點是哪個）
    // 取出S對應I中的所有列
    for (int k = 0; k < I.size(); k++) {
        int sum = 0;
        //sum = I[k];
        sum = I.get(k);
        for (int i = 0; i < paths.size(); i++) {
            Sp[i][sum] = S[i][sum];
        }
    }
    // 按行查找最大值,并記錄最大值列索引
    for(int i=0;i<paths.size();i++){
        double max = Double.NEGATIVE_INFINITY; // 定義成負無窮大方便查找最大值
        int maxIndex = -1; // 記錄列索引
        for(int j=0;j<paths.size();j++){
            if(Sp[i][j] > max){
                max = Sp[i][j];
                //maxIndex = j+1;
                maxIndex = j;
            }
        }
        c[i] = maxIndex;
    }
    // 將c數組內的值換成I中的序號(索引)
    for(int i=0;i<c.length;i++){
        int numc = c[i];
        for(int j=0;j<I.size();j++){
            if(I.get(j) == numc){
                c[i] = j;
            }
        }
    }
    // 將 1:K 的值賦給向量 c 的索引為 I 的位置（c(I)=1:K;）
    for (int i = 0; i < I.size(); i++) {
        c[I.get(i)] = i;
    }

    // 細化聚類中心
    // 當k=2時，在c中找到個體1、2、4都以2為聚類中心，提取S中1、2、4行和1、2、4列組成的3*3矩陣
    // 分別算出3列之和取最大值，y記錄最大值，j記錄最大值所在的行
    for(int k = 0; k < K; k++){
        ArrayList<Integer> ii = new ArrayList<>();
        int sum=0;
        // 首先找到同一聚類中心的每個個體
        for (int m = 0; m < c.length; m++) {
            if(c[m] == k){
                ii.add(m);
                sum++;
            }
        }
        // 提取這些個體在S中的行列，生成子矩陣
        double[][] jj = new double[ii.size()][ii.size()];
        for(int m = 0; m < ii.size(); m++){
            for(int n = 0; n < ii.size(); n++){
                jj[m][n] = S[ii.get(m)][ii.get(n)];
            }
        }
        // 并計算三列之和，記錄最大值與最大值的列
        int maxCols = jj.length; // 由于不知道jj為多大，首先找到jj的大小
        int[] caltemp = new int[maxCols]; // 用來記錄每列最大值，且最大值所在的列數為索引
        int[] yCols = new int[1]; // 記錄最大值所在的列
        double y = Double.NEGATIVE_INFINITY; // 記錄最大值
        for(int m=0;m<maxCols;m++){
            for(int n=0;n<maxCols;n++){
                caltemp[m] += jj[m][n]; // 將jj一列的值相加存入caltemp數組
            }
        }
        for(int m=0;m<caltemp.length;m++){
            if(caltemp[m] > y){
                y = caltemp[m]; // 最大值
                yCols[0] = m; // 最大值所在的列
            }
        }
        I.set(k, ii.get(yCols[0])); // 將子矩陣中列和最大的列對應的點在原矩陣中的索引賦值給I
    }
    // [tmp,c]=max(S(:,I),[],2);
    // 在矩陣S的指定列I中，逐行尋找最大值，并返回這些最大值及其對應的列索引
    double[] tmp = new double[paths.size()];
    for (int row = 0; row < S.length; row++) {
        double maxVal = Double.NEGATIVE_INFINITY;
        int maxIndex = -1;
        // 遍歷 I 中的列索引
        for (int colIndex = 0; colIndex < I.size(); colIndex++) {
            int actualCol = I.get(colIndex); // 轉換為實際列索引
            if (S[row][actualCol] > maxVal) {
                maxVal = S[row][actualCol];
                maxIndex = colIndex; // 局部索引
            }
        }
        // 保存最大值和對應列索引
        tmp[row] = maxVal;
        c[row] = maxIndex;
    }
    // c(I)=1:K;
    for (int i = 0; i < I.size(); i++) {
        c[I.get(i)] = i;
    }
    for(int i = 0; i < paths.size(); i++){
        tmpidx[i] = I.get(c[i]); // 每個點的聚類中心是誰
    }
    // 將各點到簇中心的歐式距離負值的和來衡量這次聚類的適合度
    for (int i = 0; i < tmpidx.length; i++) {
        //int row = (int)tmpidx[i] - 1; // 轉換為 0-based 索引
        int row = (int)tmpidx[i];
        tmpnetsim += S[row][i]; // 累加第 row 行的第 i 列的元素
    }
    // 表示所有被選為簇中心的點的適合度之和
    for (int i = 0; i < I.size(); i++) {
        tmpexpref += dS[I.get(i)][I.get(i)];
    }
    //int qwe = 0; // 測試
}
    else{
        Arrays.fill(tmpidx, Double.NaN); // 用 Double.NaN 填充數組
        tmpnetsim = Double.NaN;
        tmpexpref = Double.NaN;
    }

    double netsim = tmpnetsim; // 反應這次聚類的適合度
    double expref = tmpexpref;
    double dpsim = tmpnetsim - tmpexpref;
    double[] idx = new double[tmpidx.length];
    idx = Arrays.copyOf(tmpidx, tmpidx.length);
    // 返回數組 idx 中的唯一元素（去重后的元素），并按升序排列
    // 使用 HashSet 去重
    Set<Double> uniqueSet = new HashSet<>();
    for (double element : idx) {
        uniqueSet.add(element);
    }
    // 將去重后的元素轉換為數組
    double[] ans = uniqueSet.stream().mapToDouble(Double::doubleValue).toArray();
    // 對數組進行排序（升序）
    Arrays.sort(ans);

    // 根據聚類結果，將數據分配到各個簇中，并提取每個簇的中心點信息
    ArrayList<ArrayList<Path>> cluster = new ArrayList<>(); // 用來儲存最后的每個簇中的路徑(包含所有屬性)
    int tt = 0; // 計數器，記錄當前是哪個簇

    // 開始遍歷所有的簇（idx中是所有個體歸屬于哪個聚類中心，ans是整個種群選出來的聚類中心）
    for(int k = 0; k < ans.length; k++){
        ArrayList<Integer> ii = new ArrayList<>();
        // 按照ans的順序找到當前簇的所有點存入ii
        for(int i = 0; i < idx.length; i++){
            if(idx[i] == ans[k]){
                ii.add(i);
            }
        }
        // 將當前簇的簇中心點表示的路徑存入center
        int aum = 0;
        aum = (int) ans[tt];
        paths.get(aum).center = 1; // 若該條路徑是聚類中心點，則修改其center屬性為1
        // 將當前簇的每個點所對應的路徑存入cluster中
        ArrayList<Path> interRow = new ArrayList<>(); // 儲存當前簇內部的路徑(包含屬性)，并確定該聚類中的聚類中心點路徑
        for(int j = 0; j < ii.size(); j++){
            int num = ii.get(j);
            interRow.add(paths.get(num));
        }
        cluster.add(interRow);
        tt++;
    }

    //return result;//測試
    return cluster;

最后所生成的聚類為下圖所示：

至此整個聚類操作完成。

4. 交叉變異操作

4.1 交叉操作

ClusteringGA的交叉操作，針對于每一個簇中的個體，分為三種情況

1、當簇中只有1個個體時：

不進行交叉操作

2、當簇中有且只有2個個體時：

父代個體1的前部與父代個體2的后部交叉生成子代個體1

父代個體2的前部與父代個體1的后部交叉生成子代個體2

3、當簇中個體大于等于3個個體時:

父代個體1的前部與父代個體2的后部交叉生成子代個體1

父代個體2的前部與父代個體3的后部交叉生成子代個體2

……

最后一個父代個體的前部與父代個體1的后部交叉生成最后一個子代個體

// 交叉
switch (ss) {
    //  case 1:若當前聚類只有一個個體則不進行交叉操作
    case 2: // 若當前聚類只有兩個個體
        List<Integer> Ppath1 = APCluster.get(i).get(0).path;
        List<Integer> Ppath2 = APCluster.get(i).get(1).path;
        // 由第一個個體前部與第二個個體后部交叉形成第一個個體的子代
        APCluster.get(i).get(0).path = ClusteringGACrossover.Crossover(Ppath1, Ppath2);
        // 由第二個個體前部與第一個個體后部交叉形成第二個個體的子代
        APCluster.get(i).get(1).path = ClusteringGACrossover.Crossover(Ppath2, Ppath1);
        Non_repeating_sets.add(APCluster.get(i).get(0));
        Non_repeating_sets.add(APCluster.get(i).get(1));
        break;
    default: // 若當前聚類有三個及以上個體
        for (int j = 0; j < APCluster.get(i).size(); j++) {
            if (j != APCluster.get(i).size() - 1) { // 最后一個元素之前的元素，兩兩配對交叉
                List<Integer> Ppath_1 = APCluster.get(i).get(j).path;
                List<Integer> Ppath_2 = APCluster.get(i).get(j + 1).path;
                //第一個個體與第二個個體交叉生成第一個個體的子代
                APCluster.get(i).get(j).path = ClusteringGACrossover.Crossover(Ppath_1, Ppath_2);
                Non_repeating_sets.add(APCluster.get(i).get(j));
            } else { // 最后一個元素與第一個元素交叉
                List<Integer> Ppath_1 = APCluster.get(i).get(j).path; // 最后一個元素
                List<Integer> Ppath_2 = APCluster.get(i).get(0).path; // 第一個元素
                APCluster.get(i).get(j).path = ClusteringGACrossover.Crossover(Ppath_1, Ppath_2);
                Non_repeating_sets.add(APCluster.get(i).get(j));
            }
        }
        break;
}

4.2變異操作

變異操作按照正常GA的變異操作進行，在此不做過多描述。

5 快速非支配排序

5.1 拆解聚類

在進行快速非支配排序之前，首先拆除聚類，將種群換回最開始的矩陣列表

// 將聚類拆開形成原來的Path對象
public static List<Path> Discluster(ArrayList<ArrayList<Path>> APcluster){
    List<Path> paths = new ArrayList<>();
    int n = 0;
    for (int i = 0; i < APcluster.size(); i++) {
        for (int j = 0; j < APcluster.get(i).size(); j++) {
            paths.add(APcluster.get(i).get(j));
        }
    }
    return paths;
}

5.2 快速非支配排序

快速非支配排序按照正常GA的快速非支配排序進行，在此不做過多描述。

6 總結

至此，整個ClusteringGA便介紹完了，對于AP聚類算法，不可否認有很多的優點，但是或許對于多模態問題來說，并不是一個很有效的辦法，存在許多的缺點：

AP算法需要事先計算每對數據對象之間的相似度，如果數據對象太多的話，AP算法的時間復雜度較高，一次迭代大概O(N^3)
聚類的好壞受到參考度和阻尼系數的影響。
對于ClusteringGA算法，普通的GA算法每個個體都可以進行交叉，而ClusteringGA的交叉僅僅只是針對每一個聚類中的個體，破壞了多樣性，以至于其實看一些測試問題的結果ClusteringGA并不能找到所有的解。

7 參考文獻

https://www.science.org/doi/10.1126/science.1136800

https://link.springer.com/chapter/10.1007/978-3-030-95391-1_36

https://blog.csdn.net/qq_38195197/article/details/78136669

https://blog.csdn.net/quicmous/article/details/129489166

posted @ 2024-12-30 20:10 WUST許志偉閱讀(53) 評論(0) 收藏舉報

刷新頁面返回頂部

cloud&ken

論文研讀-ClusteringGA研讀與實現

論文研讀-ClusteringGA研讀與實現

1. 前言

2. 聚類(Cluster)的介紹

3. AP聚類

3.1 AP聚類的介紹

3.1.1 AP聚類的簡介

3.1.2 AP聚類的優點

3.2 AP聚類的流程

3.2.1 聚類前期準備工作

3.2.2 相似度/參考度矩陣(Similarity Matrix)

3.2.3 責任度/吸引度矩陣(Responsibility Matrix)

3.2.4 可用度/歸屬度矩陣(Availability Matrix)

3.2.5 最終判斷

3.2.6 聚類后續完善

4. 交叉變異操作

4.1 交叉操作

4.2變異操作

5 快速非支配排序

5.1 拆解聚類

5.2 快速非支配排序

6 總結

7 參考文獻

公告