K均值算法及利用其图片压缩

算法实现

一、K均值算法
K均值算法是典型的原型聚类算法，此类算法假设聚类结构能够通过一组原型刻画。通常情况下，算法先对原型进行初始化，然后对原型进行迭代更新求解。
K均值算法主要过程概括为：
1>初始化均值向量{${\mu }_1,{\mu }_2,..,{\mu }_k$},$k$值由用户给出，为聚类簇数。
2>对样本进行归类：对给定数据集$X$中每个样本$x$利用距离度量$d=||x-\mu |{|}_2$计算每个样本到均值向量的距离，把样本$x$归为其距离最近的均值向量所对应的簇。比如：对第$i$个样本$x_i$计算其到每个均值向量${\mu }_1,{\mu }_2,...，{\mu }_k$的距离，如果他到第$j$簇$C_j$,就把$x_i$划分到$C_j$中。
3>更新均值向量：现在我们每个均值向量${\mu }_j$对应的簇$C_j$下都有我们第二步通过归类所得到的样本。我们对每个簇下的样本进行求平均得到：${\mu }_j^{{}^{\prime }}=\frac{1}{|C_j|}{\sum }_{x\in C_i}x$。这样我们就可以得到把${\mu }_j^{{}^{\prime }}$当做新的均值向量${\mu }_j$。
4>跳至第二步，进行迭代，直至满足迭代终止条件为止。一般迭代终止条件为第二步中所有样本的类所属簇不在改变，或者由用户设定迭代次数。

二、K均值算法的实现
首先，根据我们上面的算法可以知道主要是对样本进行归类、更新均值向量，这样我们就需要编写两个函数。然后还需要将以上两个函数合并，给出迭代终止条件。
1、导入需要使用的模块

import numpy as np
from scipy.io import loadmat
import pandas as pd
import matplotlib.pyplot as plt

2、初始化均值向量
这里我们还需要编写一个能够初始化均值向量的函数。聚类的最终结果一定程度上取决于我们初始均值向量的选取，我们这里采用随机方法，产生初始均值向量。

def initial_centriod(X,k):
    m,d = X.shape
    initial_centriod=np.zeros((k,d))
    index=np.random.randint(0,m,k)
    for i in range(k):
        initial_centriod[i,:]=X[index[i],:]
    return initial_centriod

代码第三行创建一个$k\ast d$零矩阵用于存放均值向量，其中$k$是我们设定的聚类簇数，$d$为样本的维度或者说是我们研究的特征空间的维度，矩阵每一行就是用来存放初始均值向量，因为我们要产生$k$个均值向量，所以矩阵为k行。
第四行我们产生一个大小为$k$，范围为$1—m$的数组，$m$为样本的大小。最后对零向量每一行进行遍历，使其每一行为$m$个样本中随机产生的一个。例如，我们把$k$设置为3，一共有20个样本，那么index就有可能是数组$[1,3,4]$,到底下对零矩阵每一行进行遍历时，初始均值向量的第一行就会是第1个样本，第二行为第3个样本，第三行为第4个样本。也就是说我们把第1，3，4个样本作为初始均值向量。

3、对每个样本进行归类

def find_closest_centriod(X,centriod):
    m=X.shape[0]
    k = centriod.shape[0]
    index_centriod=np.zeros(m)
    for i in range(m):
        minidist = 1000000000000
        for j in range(k):
             dist = (centriod[j]-X[i]).T @ (centriod[j]-X[i])
             if dist < minidist:
                minidist = dist
                index_centriod[i]=j
    return index_centriod

这里代码第4行的index_centriod它的维度为训练集大小$m$，用来存放着$m$个样本的所属簇的编号。然后我们队每一个样本进行遍历，在对每一个均值向量进行遍历，找出每一个样本最近的均值向量，第8行即是在计算当前遍历样本与当前遍历均值向量的距离平方大小。找出$x_i$的最近均值向量${\mu }_j$后，即把index_centriod的第i个位置设置为均值向量${\mu }_j$对应的簇的编号j。

4、计算每个簇下样本的平均值，得到新的均值向量

def compute_centroid(X,index,k):
    m,d = X.shape
    centriod_mat = np.zeros((k,d))
    for j in range(k):
        l=[]
        for i in range(m):
            if index[i]==j:
                l.append(X[i])
        centriod_mat[j,:]=(sum(l))/len(l)
    return centriod_mat

这里我们遍历每个簇，对于每个簇创建一个列表，用于每个簇对应的样本，遍历每个上面分类得到的关于每个样本所属类别的数组，如果数组第$i$个位置对应的数为第$j$个簇对应的编号，那么就相当于我们的第i个样本$x_i$属于第$j$簇，就将其存放在列表中，最后对$j$簇对应的列表中的所有样本求平均，即为新的均值向量${\mu }_j$，将其存放在已经创建的矩阵的第$j$行，这个矩阵的第$i$行存放着第$i$簇的均值向量。

5、将两个函数整合在一起，即主函数

def run_k_means(X,centriod,max_iter):
    old_centriod=centriod
    for i in range(max_iter):
        old_index=find_closest_centriod(X,old_centriod)
        k = old_centriod.shape[0]
        new_centriod=compute_centroid(X,old_index,k)
        new_index=find_closest_centriod(X,new_centriod)
        if list(new_index)==list(old_index):
            break
        else:
            old_centriod=new_centriod
    return new_index,new_centriod

参数centriod为初始均值向量，max_iter为用户设置的迭代最大次数。如果我们发现迭代前后两次对样本分类结果一样，就立即停止迭代，否则执行到最大迭代次数max_iter。

对数据进行聚类

一、处理数据并查看

先查看数据的维度：

data = loadmat('F:\\MachineLearning\data\ex7data2.mat')
X=data['X']
print(X.shape)

得到

(300,2)

即数据为二维，可以对其进行可视化：

#查看数据分布
data2 = pd.DataFrame(X,columns=['x1','x2'])
fig1,ax1=plt.subplots(figsize=(9,6))
ax1.scatter(data2['x1'],data2['x2'],s=20,color='blue')
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()

得到数据分布图：

很明显，我们可以把数据分为三类。

二、对数据进行聚类

1、首先，我们产生初始均值向量，用上面的initial_centriod函数：

centriod = initial_centriod(X,3)

2、接着，运行主函数进行聚类，并处理数据：

index,centriod=run_k_means(X,centriod,10)
data2['index']=index
data_index1=data2[data2['index']==0]
data_index2=data2[data2['index']==1]
data_index3=data2[data2['index']==2]

，这里我们设置最大迭代次数为10，然后得到了聚类原型，也就是最后一次迭代的均值向量centriod，以及每个样本聚类的结果index。然后我们给数据添加一列，即为index，给出每个样本对应的聚类编号。然后分别把属于簇1、2、3的样本分别存放data_index1、data_index2、data_index3。
3、查看聚类结果

fig2,ax2 = plt.subplots(figsize=(9,6))
ax2.scatter(data_index1['x1'],data_index1['x2'],s=20,color='blue',label='cluster1')
ax2.scatter(data_index2['x1'],data_index2['x2'],s=20,color='red',label='cluster2')
ax2.scatter(data_index3['x1'],data_index3['x2'],s=20,color='green',label='cluster3')
plt.legend()
plt.show()

这里我们将簇1用蓝色表示，簇2用红色表示，簇3用绿色表示：

发现聚类结果良好。
当然也会发现例如一下这样的情况：

这样的聚类结果显然是不能令人满意的，这主要是因为我们的初始均值向量的选择会影响聚类结果，然而我们的均值向量是随机产生的，会有一定可能性（可能性较小）选择了“不好”的初始均值向量。如果遇到这种情况我们可以多次运行程序，当然也可以根据数据分布自行选择。

利用聚类算法进行图像压缩

一、查看原图，并产看数据
首先，查看我们需要处理的原图：

from IPython.display import Image,display
display(Image(filename='F:/MachineLearning/data/bird_small.png'))

得到原图：

接着，处理数据：
先查看数据维度

image_data=loadmat('F:/MachineLearning/data/bird_small.mat')
A=image_data['A']
print(A.shape)

(128, 128, 3)

即数据为三维数组。前两个维度即128和128表示我们的平面图由128128个像素点，第三个维度3表示我们每个像素点的红、蓝、绿的强度。所以我们有128128个样本，样本的维度为3。而我们利用聚类算法进行压缩的原理主要是将我们128*128个样本处理成k个簇，再把每个样本代表的簇的原型代替原先样本，利用这些点组成新的图。

二、利用聚类算法进行压缩图片

A=A/255
X=np.reshape(A,(A.shape[0]*A.shape[1],A.shape[2]))
initial_image_centriod=initial_centriod(X,16)
index,centriod=run_k_means(X,initial_image_centriod,10)
X_recovered = centriod[index.astype(int),:]
X_recovered=np.reshape(X_recovered,(A.shape[0],A.shape[1],A.shape[2]))
plt.imshow(X_recovered)
plt.show()

第一行：归一化数据，让imshow函数可以处理。第二行：我们把样本展开，形成128*128行，3列的样本矩阵，可以让程序处理。第三行：初始化聚类均值向量，这里我们聚类簇数设置为16，簇数决定了压缩图片的还原原图的质量，一般这个值设置的越大，质量越好，同时图片内存也越大。第四行：运行程序。第五行：将原始数据点用聚类原型代替。第六行：再次改变维度使变为图片维度。第七行：形成压缩图片。最后，显示图片。得到图：

这是我们聚类簇数为16的图。如果我们将第三行聚类簇数设置为32、8和4，会看到如下三图：



会发现参数为32好于参数16，参数为16明显好于参数为8，参数为8明显好于参数为4。这不难理解，聚类簇数相当压缩图片所使用的色彩数，色彩越少，图片质量肯定越差。

三、用sklearn中K-Means库处理此问题

from skimage import io
pic=io.imread('F:/MachineLearning/data/bird_small.png')
pic=pic/255
print(pic.shape)
io.imshow(pic)
plt.show

这里我们直接从原图得到数据，第一行：导入需要的模块。第二行：io.imread将我们的原图处理为数组形式。第三行：归一化。第四行：查看数据维度，为(128, 128, 3)。第五、六行再利用得到数据重构原图。我们会得到：

这实际上就是原图。

from sklearn.cluster import KMeans
#系数解释：n_cluster:聚类数。n_init：使用初始点次数。n_jobs：是否使用多线程，-1使用。
Cluster = KMeans(n_clusters=4,n_init=100,n_jobs=-1)
Cluster.fit(X)
centriod_pic=Cluster.cluster_centers_
print(centriod_pic)
C=Cluster.predict(data_pic)
pic_recovered=centriod_pic[C].reshape((128,128,3))

上面是调用K-Means库的情况。第一行：导入模块。第三行：设置参数。第四行：拟合数据。第五行：读取聚类原型矩阵。第七行：读取聚类结果C，即每个样本属于簇组成的数组。第八行：再把数据处理成可以图片数组。
这里我们可以得到centriod_pic，聚类原型矩阵：

[[0.48994379 0.40168338 0.32354671]
 [0.1287319  0.13038871 0.12085462]
 [0.91127734 0.8588342  0.74082312]
 [0.79460449 0.63997646 0.41903748]]

最后查看原图和压缩图片的差别：