FCM聚類演算法詳解(Python實現iris資料集)

參考：https://blog.csdn.net/on2way/article/details/47087201

模糊C均值（Fuzzy C-means）演算法簡稱FCM演算法，是一種基於目標函式的模糊聚類演算法，主要用於資料的聚類分析。理論成熟，應用廣泛，是一種優秀的聚類演算法。本文關於FCM演算法的一些原理推導部分介紹，加上自己的理解和在課題專案中的應用以文字的形式呈現出來。

首先介紹一下模糊這個概念，所謂模糊就是不確定，確定性的東西就是資訊量很小的東西，而不確定性的東西就說很像什麼，這種資訊量很大。比如說把20歲作為年輕不年輕的標準，那麼一個人21歲按照確定性的劃分就屬於不年輕，而我們印象中的觀念是21歲也很年輕，這個時候可以模糊一下，認為21歲有0.9分像年輕，有0.1分像不年輕，這裡0.9與0.1不是概率，而是一種相似的程度，把這種一個樣本屬於結果的這種相似的程度稱為樣本的隸屬度，一般用u表示，表示一個樣本相似於不同結果的一個程度指標。

基於此，假定資料集為X，如果把這些資料劃分成c類的話，那麼對應的就有c個類中心為C，每個樣本j屬於某一類i的隸屬度為uij，那麼定義一個FCM目標函式（1）及其約束條件（2）如下所示：
$J=∑_{i=1}^{n}∑_{j=1}^{n}u_{ij}^{m}||xj−ci||^2 ...............................(1)$

$∑_{i=1}^{c}u_{ij}=1,j=1,2...,n..............................(2)$

看一下目標函式（式1）而知，由相應樣本的隸屬度與該樣本到各個類中心的距離相乘組成的，m是一個隸屬度的因子，個人理解為屬於樣本的輕緩程度，就像$x^2$與$x^3$這種一樣。式（2）為約束條件，也就是一個樣本屬於所有類的隸屬度之和要為1。觀察式（1）可以發現，其中的變數有$u_{ij}、c_i$，並且還有約束條件，那麼如何求這個目標函式的極值呢？

這裡首先採用拉格朗日乘數法將約束條件拿到目標函式中去，前面加上係數，並把式（2）的所有j展開，那麼式（1）變成下列所示：

$J=∑_{i=1}^{c}∑_{j=1}^{n}u_{ij}^{m}||x_j−c_i||^2+λ_1(∑_{i=1}^{c}u_{i1}−1)+...+λ_j(∑_{i=1}^{c}u_{ij}−1)+...+λ_n(∑_{i=n}cu_{in}−1)).....................................(3)$

現在要求該式的目標函式極值，那麼分別對其中的變數uij、ci求導數，首先對uij求導。

分析式(3)，先對第一部分的兩級求和的uij求導，對求和形式下如果直接求導不熟悉，可以把求和展開如下：

這個矩陣要對uij求導，可以看到只有uij對應的$u_{m}^{ij}||x_j−c_i||^2$保留，其他的所有項中因為不含有uijuij，所以求導都為0。那麼$u_{ij}^{m}||x_j−c_i||^2$對uij求導後就為$m||x_j−c_i||^2u_{ij}^{m-1}$。

再來看後面那個對uij求導，同樣把求和展開，再去除和uij不相關的（求導為0），那麼只剩下這一項：$λ_j(u_{ij}−1)$，它對uij求導就是λj了。

那麼最終J對uij的求導結果並讓其等於0就是：
$\frac{∂J}{∂uij}=m||x_j−c_i||^2u_{ij}^{m-1}+λ_j=0$

這個式子化簡下，將uij解出來就是：

$u_{ij}^{m-1}=\frac{−λ_j}{m||xj−ci||^2}$

進一步：

$u_{ij}=(\frac{−λ_j}{m||xj−ci||^2})^\frac{1}{m−1}=(\frac{−λ_j}{m})^\frac{1}{m−1}\frac{1}{(\frac{1}{||xj−ci||})^\frac{2}{m−1}}.......(4)$

要解出uij則需要把λj去掉才行。這裡重新使用公式（2）的約束條件，並把算出來的uij代入式（2）中有：

這樣就有（其中把符號i換成k）：
把這個重新代入到式（4）中有：

好了，式子（5）就是最終的uij迭代公式。

下面在來求J對ci的導數。由公式（2）可以看到只有$∑_{i=1}^{c}∑_{j=1}^{n}u_{ij}^{m}||xj−ci||^2$這一部分裡面含有ci，對其二級求和展開如前面所示的，那麼它對ci的導數就是：

即：


好了，公式（6）就是類中心的迭代公式。

我們發現uij與ci是相互關聯的，彼此包含對方，有一個問題就是在fcm演算法開始的時候既沒有uij也沒有ci，那要怎麼求解呢？很簡單，程式開始的時候我們隨便賦值給uij或者ci其中的一個，只要數值滿足條件即可。然後就開始迭代，比如一般的都賦值給uij，那麼有了uij就可以計算ci，然後有了ci又可以計算uij，反反覆覆，在這個過程中還有一個目標函式J一直在變化，逐漸趨向穩定值。那麼當J不在變化的時候就認為演算法收斂到一個比較好的解了。可以看到uij和ci在目標函式J下似乎構成了一個正反饋一樣，這一點很像EM演算法，先E在M，有了M在E，在M直至達到最優。

公式（5），（6）是演算法的關鍵。現在來重新從巨集觀的角度來整體看看這兩個公式，先看（5），在寫一遍

假設看樣本集中的樣本1到各個類中心的隸屬度，那麼此時j=1，i從1到c類，此時上述式中分母裡面求和中，分子就是這個點相對於某一類的類中心距離，而分母是這個點相對於所有類的類中心的距離求和，那麼它們兩相除表示什麼，是不是表示這個點到某個類中心在這個點到所有類中心的距離和的比重。當求和裡面的分子越小，是不是說就越接近於這個類，那麼整體這個分數就越大，也就是對應的uij就越大，表示越屬於這個類，形象的圖如下：

再來巨集觀看看公式（6），考慮當類i確定後，式（6）的分母求和其實是一個常數，那麼式（6）可以寫成：

這是類中心的更新法則。說這之前，首先讓我們想想kmeans的類中心是怎麼更新的，一般最簡單的就是找到屬於某一類的所有樣本點，然後這一類的類中心就是這些樣本點的平均值。那麼FCM類中心怎麼樣了？看式子可以發現也是一個加權平均，類i確定後，首先將所有點到該類的隸屬度u求和，然後對每個點，隸屬度除以這個和就是所佔的比重，乘以xj就是這個點對於這個類i的貢獻值了。畫個形象的圖如下：

由上述的巨集觀分析可知，這兩個公式的迭代關係式是這樣的也是可以理解的。

資料集用的是iris資料，有需要資料的朋友可以呼叫sklearn.load_iris。或者下載下來

from pylab import *
from numpy import *
import pandas as pd
import numpy as np
import operator
import math
import matplotlib.pyplot as plt
import random

# 資料儲存在.csv檔案中
df_full = pd.read_csv("iris.csv")
columns = list(df_full.columns)
features = columns[:len(columns) - 1]
# class_labels = list(df_full[columns[-1]])
df = df_full[features]
# 維度
num_attr = len(df.columns) - 1
# 分類數
k = 3
# 最大迭代數
MAX_ITER = 100
# 樣本數
n = len(df)  # the number of row
# 模糊引數
m = 2.00


# 初始化模糊矩陣U
def initializeMembershipMatrix():
    membership_mat = list()
    for i in range(n):
        random_num_list = [random.random() for i in range(k)]
        summation = sum(random_num_list)
        temp_list = [x / summation for x in random_num_list]  # 首先歸一化
        membership_mat.append(temp_list)
    return membership_mat


# 計算類中心點
def calculateClusterCenter(membership_mat):
    cluster_mem_val = zip(*membership_mat)
    cluster_centers = list()
    cluster_mem_val_list = list(cluster_mem_val)
    for j in range(k):
        x = cluster_mem_val_list[j]
        xraised = [e ** m for e in x]
        denominator = sum(xraised)
        temp_num = list()
        for i in range(n):
            data_point = list(df.iloc[i])
            prod = [xraised[i] * val for val in data_point]
            temp_num.append(prod)
        numerator = map(sum, zip(*temp_num))
        center = [z / denominator for z in numerator]  # 每一維都要計算。
        cluster_centers.append(center)
    return cluster_centers


# 更新隸屬度
def updateMembershipValue(membership_mat, cluster_centers):
    #    p = float(2/(m-1))
    data = []
    for i in range(n):
        x = list(df.iloc[i])  # 取出檔案中的每一行資料
        data.append(x)
        distances = [np.linalg.norm(list(map(operator.sub, x, cluster_centers[j]))) for j in range(k)]
        for j in range(k):
            den = sum([math.pow(float(distances[j] / distances[c]), 2) for c in range(k)])
            membership_mat[i][j] = float(1 / den)
    return membership_mat, data


# 得到聚類結果
def getClusters(membership_mat):
    cluster_labels = list()
    for i in range(n):
        max_val, idx = max((val, idx) for (idx, val) in enumerate(membership_mat[i]))
        cluster_labels.append(idx)
    return cluster_labels


def fuzzyCMeansClustering():
    # 主程式
    membership_mat = initializeMembershipMatrix()
    curr = 0
    while curr <= MAX_ITER:  # 最大迭代次數
        cluster_centers = calculateClusterCenter(membership_mat)
        membership_mat, data = updateMembershipValue(membership_mat, cluster_centers)
        cluster_labels = getClusters(membership_mat)
        curr += 1
    print(membership_mat)
    return cluster_labels, cluster_centers, data, membership_mat


def xie_beni(membership_mat, center, data):
    sum_cluster_distance = 0
    min_cluster_center_distance = inf
    for i in range(k):
        for j in range(n):
            sum_cluster_distance = sum_cluster_distance + membership_mat[j][i] ** 2 * sum(
                power(data[j, :] - center[i, :], 2))  # 計算類一致性
    for i in range(k - 1):
        for j in range(i + 1, k):
            cluster_center_distance = sum(power(center[i, :] - center[j, :], 2))  # 計算類間距離
            if cluster_center_distance < min_cluster_center_distance:
                min_cluster_center_distance = cluster_center_distance
    return sum_cluster_distance / (n * min_cluster_center_distance)


labels, centers, data, membership = fuzzyCMeansClustering()
print(labels)
print(centers)
center_array = array(centers)
label = array(labels)
datas = array(data)

# Xie-Beni聚類有效性
print("聚類有效性：", xie_beni(membership, center_array, datas))
xlim(0, 10)
ylim(0, 10)
# 做散點圖
fig = plt.gcf()
fig.set_size_inches(16.5, 12.5)
f1 = plt.figure(1)
plt.scatter(datas[nonzero(label == 0), 0], datas[nonzero(label == 0), 1], marker='o', color='r', label='0', s=10)
plt.scatter(datas[nonzero(label == 1), 0], datas[nonzero(label == 1), 1], marker='+', color='b', label='1', s=10)
plt.scatter(datas[nonzero(label == 2), 0], datas[nonzero(label == 2), 1], marker='*', color='g', label='2', s=10)
plt.scatter(center_array[:, 0], center_array[:, 1], marker='x', color='m', s=30)
plt.show()

效果一般啊。。。
竟然聚類的結果感覺和剛開始的好像差不多啊。不會被聚成一類？但是回過頭來看這個演算法，本身要得到的結果就是隸屬度u（也就是各樣本的權重）以及Ci聚類中心，有幾個聚類中心就有幾類，可以看到，還是挺準的。