基于Python实现聚类算法

一、理论知识

1.1 K-Means

给定一组数据集，聚类算法将它们分成不同的子组。我们希望类内实例高度相似，类间实例低相似。

在样本集中，随机选取K个点作为中心${\mu }_k$，计算每个样本到中心点的距离，并将样本划分到离它最近的那个点的集群中。使用变量$r_{nk}$表示数据样本${\mathbf{x}}^{(n)}$是否属于集群k：
$$r_{nk}=\{\begin{array}{c}1,k=arg{\min }_j||{\mathbf{x}}^{(n)}-{\mu }_j|{|}^2\\ 0,otherwise\end{array}$$
对于每个集群，用所有样本的平均位置更新中心点的位置：
$${\mu }_k=\frac{{\sum }_{n=1}^N{r}_{nk}{\mathbf{x}}_n}{{\sum }_{n=1}^N{r}_{nk}}$$
重复上面的样本分配和中心更新过程即可，该过程是保证收敛的。类内距离之和会随着K的增大而减小，因此选择K的时候考虑寻找使得距离基本收敛的K值或是依据下游应用需求选择K值。

K平均聚类的结果高度依赖于初始化的中心点。选择初始化中心点的一般方法有：

• 随机选择：随机选择数据样本作为初始中心。问题在于有可能选择了距离很近的多个样本作为初始中心
• 基于距离的方法：一开始随机选一个样本作为中心点，其他时候选择离各个中心最远的样本作为中心点。问题在于可能选到异常值
• 随机+距离方法：第一个中心点为随机选择的样本，从全部远离现有中心的样本中随机选择下一个中心

一个点要么属于一个集群，要么不属于，存在赋值难的问题。可以用软K平均算法解决：
$$r_{nk}=\frac{e^{-\beta }||{\mathbf{x}}^{(n)}-{\mu }_k|{|}^2}{{\sum }_{i=1}^K{e}^{-\beta }||{\mathbf{x}}^{(n)}-{\mu }_i|{|}^2}$$
这样一来，$r_{nk}$可以理解为样本属于某个集群的几率。

该算法还对异常值敏感。而且欧几里得距离决定了决策边界是球形的，当不同簇的形状不规则时，性能可能不佳。

1.2 GMM与EM算法

在监督学习中，回归和分类都可以理解为学习的条件概率分布。类似地，从概率建模的角度来看，无监督学习可以理解为对输入数据的概率分布进行学习，最简单的方法是假设该分布为高斯形式。但是一般情况下，样本的分布比高斯分布复杂得多。为了更好地建模，就需要更好的分布表示方法，即使用复杂的分布，但是训练成本又会增加。

因此，一个直观的思路是使用很多简单模型的加和表示复杂模型：
$$p(\mathbf{x})={\pi }_1{p}_1(\mathbf{x})+{\pi }_K{p}_K(\mathbf{x})$$
其中，所有$\pi$的和为1，用来使得分布合法。而每个$p_i$表示简单的分布。上面的每一项又可以看作是联合分布的边缘分布：
$$p_z(\mathbf{x}){\pi }_z=p(\mathbf{x}|z)p(z)=p(\mathbf{x},z)$$
这种模型被称为潜变量模型LVM。我们将各个分布看做是高斯分布，则该模型又称为高斯混合模型，即GMM。如果一个高斯混合模型由$K$个高斯分布组成，则总的分布表达式为：
$$p(\mathbf{x})=\sum _{k=1}^K{\pi }_{k}N(\mathbf{x};{\mu }_k,{\mathbf{\Sigma }}_k)$$
其中，$K$表示高斯分布的个数，${\pi }_k$表示第$k$个分布的权重。${\sum }_{k=1}^K{\pi }_k=1$。${\mu }_k$和${\mathbf{\Sigma }}_k$是第$k$个高斯分布的均值向量和协方差矩阵。

EM算法可以学习高斯混合模型的参数，使得该模型在特定的数据集下获得更大的对数似然函数值。算法分为E和M两步。E指评估预期，M指更新数据：

• E-step：$Q(\theta ;{\theta }^{(t)})=E_{p(\mathbf{z}|\mathbf{x};{\theta }^{(t)})}[\log p(\mathbf{x},\mathbf{z};\theta )]$
• M-step：${\theta }^{(t+1)}=\arg {\max }_{\theta }Q(\theta ;{\theta }^{(t)})$

对于上述的高斯混合模型，在E-step中，可以求得对数似然函数为：
$$\log p(\mathbf{x},\mathbf{z};\theta )=\sum _{k=1}^K{z}_k\cdot [\log N({\mathbf{x}}_n|{\mu }_k,{\mathbf{\Sigma }}_k)+\log {\pi }_k]$$
令：
$${\gamma }_k^{(t)}\triangleq \frac{N({\mathbf{x}}_n|{\mu }_k^{(t)},{\mathbf{\Sigma }}_k^{(t)}){\pi }_k}{{\sum }_{i=1}^{K}N({\mathbf{x}}_n|{\mu }_k^{(t)},{\mathbf{\Sigma }}_k^{(t)}){\pi }_i}$$
则有：
$$Q(\theta ;{\theta }^{(t)})=\sum _{k=1}^K{\gamma }_k^{(t)}[\log N({\mathbf{x}}_n|{\mu }_k,{\mathbf{\Sigma }}_k)+\log {\pi }_k]$$
然后在M-step更新各个参数：
$$\begin{gathered} {\mu }_k^{(t+1)}=\frac{1}{N_k}\sum _{n=1}^N{\gamma }_{nk}{\mathbf{x}}_n \\ {\mathbf{\Sigma }}_l^{(t+1)}=\frac{1}{N_k}\sum _{n=1}^N{\gamma }_{nk}({\mathbf{x}}_n-{\mu }_n^{(t+1)})({\mathbf{x}}_n-{\mu }_n^{(t+1)}{)}^T{\pi }_k^{(t+1)}=\frac{N_k}{N} \end{gathered}$$
其中$N_k={\sum }_{n=1}^N{\gamma }_{nk}$，表示第$k$类的样本的有效数量。反复执行E和M两步训练参数，直至模型收敛或到达最大训练次数即可。

---

二、代码实现

2.1 K-means

在K-means算法中，我们将训练数据集分成若干个聚类，对每个类打上标签后对训练样本进行分类，将属于的类的标签当做预测结果即可。

2.1.1 中心点的选择

K平均算法的核心就是中心点的选择和更新。对于本次实验使用的MNIST数据集，样本共有10类，因此聚类时选择10个不同的中心点。初始化中心点有多种方法，我会在之后的实验结果分析中分析不同的初始化方法对结果的影响。

下面的代码以随机选择为例子。在全部样本点中选择十个样本点作为初始化中心：

# 设置基本参数
k = 10  # 类别数
t = 100 # 迭代次数
# 生成初始中心
# 完全随机
center_points = np.random.choice(train_images.shape[0], k, replace = False).tolist()
for i in range(len(center_points)):
    center_points[i] = train_images[center_points[i]]
center_points = np.asarray(center_points)

同时，我还对基于距离的方法进行了尝试。即，第一个中心点为随机选择的样本点，第二个中心点为离第一个中心点最远的样本点，第三个中心点为距离第一和第二个中心点距离之和最大的样本点，以此类推：

center_points = genCenterPoints(train_images, k)

# 基于距离选择中心点
def genCenterPoints(train_images, k):
    # 第一个中心点随机选择
    center_points = [train_images[random.randint(0,train_images.shape[0]-1)].tolist()]
    while len(center_points) != k:
        print(len(center_points))
        dist = None
        point = None
        # 遍历每个样本点
        for each in train_images:
            # 不考虑已经作为中心点的样本点
            if each.tolist() in center_points:
                continue
            # 计算当前样本到各个中心点的距离之和
            tmp_dist = np.sum(np.linalg.norm(np.asarray(center_points) - each,axis = 1))
            # 选择距离最短的样本点作为下一个中心点
            if dist==None or tmp_dist<dist:
                point = each
                dist = tmp_dist
        center_points.append(point.tolist())
    return np.asarray(center_points)

2.1.2 训练过程

之后就是对各个样本点进行分类和中心点的迭代更新：

# 进行迭代
for i in range(t):
    # 对训练样本进行分类
    r = classification(center_points,train_images, k)
    # 更新中心点位置
    center_points = updateCenterPoints(r, train_images, k)

对样本点进行分类时，使用数组r[i][j]表示第i个样本是否属于第j类。是则为1，否则为0。计算出当前样本和各个中心点之间的距离，选择最近的中心点作为该样本属于的类别即可。这里我调用了np.linalg.norm来计算距离，得到的是欧式距离，也就是各个特征值的平方和开根号作为距离。样本距离可以采用不同的标准，会对分类的结果产生不同的影响。

def classification(center_points,train_images, k):
    # r[i][j]表示第i个样本是(1)否(0)属于第j类
    r = np.zeros((train_images.shape[0],k)).astype(np.int16)
    for i in range(train_images.shape[0]):
        # 计算样本距离各个中心点的距离
        distance = np.linalg.norm(center_points - train_images[i], axis=1)
        # 将相距最近的中心点所属的类作为该样本的类
        r[i][np.argwhere(distance==np.min(distance))[0][0]] = 1
    return r

将每个样本成功归类后就要更新中心点的位置。 这里的中心点位置需要到该类中各个样本点的距离最短，因此是各个样本点的平均位置。同样，这里采用不同的距离度量方式也可以产生不同的结果。

# 更新中心点位置
def updateCenterPoints(r, train_images, k):
    center_points = np.empty((k,train_images.shape[1]))
    for i in range(k):
        # 类内样本特征值之和除以样本数，即平均特征值
        center_points[i] = np.sum(train_images*(np.repeat(r[:,i].reshape(1,-1),\
train_images.shape[1],axis=0).T),axis=0) / np.sum(r[:,i])
    return center_points

2.1.3 聚类标签的获取

现在我们得到了十个聚类的中心点。经过训练，我们期望这些点能够对数据进行分类。然而我们不知道每个中心点对应的是哪一个标签。因为训练集中，每个类中有多个样本，可以让这些样本进行投票，每个样本为自身表示的标签进行投票，票数最多的标签就作为该中心点所代表的类的标签：

labels = (genLabels(r, train_labels, k)).astype(np.int16)

# 对每个聚类打上标签
def genLabels(r, train_labels, k):
    labels = np.empty(10)
    # 第一层循环遍历每个聚类
    for i in range(k):
        # 第二层循环统计该类中样本属于某一类的个数
        max_count = 0
        max_label = -1
        # 分别对每个标签进行投票
        for j in range(k):
            tmp = np.sum(train_labels[r[:,i]==1]==j)
            # 记录当前票数最多的标签与对应的票数
            if tmp > max_count:
                max_count = tmp
                max_label = j
        labels[i] = max_label
    return labels

2.1.4 验证结果

在验证时，只需要比较验证样本与各个中心点的距离，选择距离最近的中心点所代表的类的标签作为预测结果，与实际的验证样本标签相比较即可。

right = 0
for i in range(test_images.shape[0]):
    distance = np.linalg.norm(center_points - test_images[i], axis=1)
    predict = labels[np.argwhere(distance==np.min(distance))[0][0]]
    if predict == test_labels[i]:
        right+=1
print(right,right/test_images.shape[0])

2.2 GMM

在GMM方法中，我们将不同的训练样本看做一个高斯混合模型的分布。不同的标签对应不同的分布，而每个分布由多个高斯分布组合而成。因此，对于每一类的分布都需要单独训练。

2.2.1 数据预处理

在GMM算法中，需要计算每个样本属于某个分布中的高斯分布的概率。而当特征向量的维数很高时，高斯分布的概率密度函数在特定点的大小就会特别的小。本次的实验中每个样本的维度为784，这么高的维度是不能够应用GMM的，因此首先要考虑对数据进行降维。这里我采用了sklearn.decomposition中的PCA模型，使用主成成分分析的方法对数据进行降维：

pca_model = PCA(n_components=50)
pca_model.fit(train_images.T)
train_images = pca_model.components_.T
pca_model.fit(test_images.T)
test_images = pca_model.components_.T

每个GMM应该用于拟合一类的样本，因此我们需要给样本进行分类。按照标签将训练样本分为十类，每一类单独用来训练一个GMM：

# 将训练数据按照类别分组
type_num = 10   # 不同标签的个数
tmp = [[] for i in range(type_num)]
for i in range(train_num):
    tmp[train_labels[i]].append(train_images[i])
# train_images[i][j]表示类型为i的第j个样本
train_images = tmp[:][:]

2.2.2 训练过程

用每一类单独训练一个GMM：

# 设置基本参数
k = 10      # 每个GMM中高斯模型的个数
t = 100     # 单个GMM迭代次数
GMM = []    # 记录每个类的GMM
for i in range(type_num):
    GMM.append(genGMM(k,t,np.asarray(train_images[i])))

训练单个GMM时，首先设置GMM的基本参数和初始值。需要混合高斯模型中，每个高斯分布的均值gauss_means、每个高斯分布的协方差矩阵gauss_var，以及这些高斯分布在模型中所占的比例pi。同时还需要样本属于某个分布的概率值gamma。分布的均值由数据集中随机抽取k个样本得到，协方差矩阵初始化为单位矩阵。比重和概率一开始都均等。

def genGMM(k,t,train_images):
    # 每个高斯分布的均值，初始化为随机的k个样本点的位置
    train_num = train_images.shape[0]
    gauss_means = np.asarray(train_images[random.sample([i for i in range(train_num)],k)])
    # 每个高斯分布的协方差矩阵
    gauss_var = np.asarray([(1*np.eye(train_images.shape[1])) for j in range(k)])# 对角且相等
    # gamma[i][j] 表示样本i属于第j个高斯分布的概率
    gamma = np.ones((train_num, k)) / k
    # 每个分布的比重
    pi = gamma.sum(axis=0) / gamma.sum()

下面开始迭代更新参数。

for i in range(t):

本次实验采用EM算法来训练GMM。EM算法的执行过程简单来说，就是在给定的参数下得到对数似然函数的期望，然后进一步更新参数，使得原来的参数下得到的对数似然函数取得最大值。

E步骤的过程如下：

# E-step
# 计算对数似然函数
pdf = np.zeros((train_num, k))
for j in range(k):
    pdf[:, j] = pi[j] * multivariate_normal(gauss_means[j], gauss_var[j],allow_singular=True).pdf(train_images)
    log_likehood = np.mean(np.log(pdf.sum(axis=1)))
    print(log_likehood)

这里，我使用了scipy库的multivariate_normal来生成单个多元高斯分布模型。给定该分布的均值和协方差矩阵，可以得到其概率密度分布。调用pdf功能，给出样本，就可以得到该样本在该高斯分布下出现的概率。用pdf[i][j]表示在该高斯混合模型下，第i个样本属于第j个高斯分布的概率。

然后，进行M步骤，按照先前算出来的公式更新参数，使得当前的概率密度期望达到最大值：

# M-step
# 更新各个分布的参数
gamma = pdf / pdf.sum(axis=1).reshape(-1,1)
pi = gamma.sum(axis=0) / gamma.sum()
for j in range(k):
    gauss_means[j] = np.average(train_images, axis=0, weights=gamma[:, j])
    cov = [np.dot((train_images[t]-gauss_means[j]).reshape(-1,1),(train_images[t]-gauss_means[j]).reshape(1,-1)) for t in range(train_num)]
    cov_sum = np.zeros((train_images.shape[1], train_images.shape[1]))
    for l in range(train_num):
        cov_sum += gamma[l][j] * cov[j]
        gauss_var[j] = cov_sum / np.sum(gamma[:,j])

每个高斯混合模型需要保留各个高斯分布的均值、协方差，以及他们的权重，这三个参数决定了一个完整的高斯混合模型：

return [gauss_means, gauss_var, pi]

2.2.3 验证结果

在得到了所有的高斯混合模型后，需要测试时，将测试样本在每个高斯混合模型下算出其对数似然函数值，用于表示它符合各个高斯混合模型的概率。选择概率最高的那一个模型对应的标签作为预测结果：

right = 0
for l in range(test_images.shape[0]):
    prob = []
    # 依次使用每个GMM进行预测
    for i in range(type_num):
        pdf = np.zeros(k)
        # 用每个模型中的每个高斯分布得到样本分布概率
        for j in range(k):
            pdf[j] = GMM[i][2][j] * multivariate_normal(GMM[i][0][j], GMM[i][1][j],allow_singular=True).pdf(test_images[l])
            # 计算对数似然
        prob.append(pdf.sum())
    prob = np.asarray(prob)
    # 选择对数似然最高的GMM对应的标签作为预测结果
    ans = np.argwhere(prob==np.max(prob))[0][0]
    if ans == test_labels[l]:
        right+=1
print(right, right/test_images.shape[0])

---

三、实验结果与分析

3.1 K-means

使用MNIST数据集，同样进行100次迭代，K-means算法中，如果使用随机初始化中心点，准确率为58.9%；如果是基于距离选择的中心点，准确率为59.45%。二者差别不大。

K平均聚类的结果高度依赖于初始化的中心点。选择中心点的策略不同的话可能会得到不同的结果。我们希望选择的初始点有各个类别的代表性，能够很好地区分各个类别。这就要求选择的初始点相距较大。初始时如果是随机选择中心点，即随机选择数据样本作为初始中心，有可能选择了距离很近的多个样本作为初始中心从而导致结果不好。当然也有可能随机选择的中心点较优。总的来说，存在选择较差点的风险。为了避免选到太近的点作为中心点，可以使用基于距离的方法，即一开始随机选一个样本作为中心点，其他时候选择离各个中心最远的样本作为中心点。但是如果这样选择的话，问题在于可能选到异常值，导致聚类结果不准确。更好的方法是结合二者的思想，设定一个距离阈值，在所有超过距离阈值的点中随机选择，这样既避免了选择的点太近，又可以减少选到异常点的概率。

MNIST数据集为手写数字图片的集合。可以想见，相同的位置的像素点，不同的数字更有可能不同，而相同数字更有可能相同。K平均算法处理这样的数据是合理且有效的。

3.2 GMM

使用EM训练的GMM算法，初始化数据降维到50，每个GMM含有50个高斯分布，迭代10次，如果协方差矩阵初始化为对角且元素相等，准确率为37.86%；如果对角且不相等则为35.23%。下面的分析以协方差矩阵初始化为对角且元素相等的情况为例。

下图为第一个GMM训练时，对数似然函数值关于训练次数的变化：

模型很快收敛了，但最终的对数似然函数的数值仍然较小，且几乎收敛，说明模型已经尽可能地拟合了数据，但拟合效果仍然不理想，模型不能很好地表示数据的分布。通过输出结果，我发现0、6等标签出现的概率高于其他标签，这说明这几个标签对应的GMM能更好地表示整个数据集，因此在这几个模型下会有更高的对数似然值，甚至比实际标签对应的模型算出来的对数似然值更高。这也侧面体现了模型拟合性能的不足。

可以看到，GMM算法的准确率非常低，而且实际训练过程中，训练时间相较K-means也更长。这应该是数据维度导致的。本次实验处理的数据维度高达784，而单个多变量高斯函数随着维度的上升，概率密度函数会越来越扁平。也就是说，在高维空间上，由高斯分布拟合的数据的概率会非常非常小，即便是在概率最大的均值处也会相当小。高达784维的数据根本不能简单地使用GMM算法运算。使用PCA方法降维又会产生样本信息的损失，且样本本身也不一定符合高斯分布，由此导致了结果性能较差。GMM更适合样本特征向量维度较小的情形。

总的来说，在MNIST数据集下，相比EM训练的GMM算法，K-means算法有着更好的性能。