本文已参与「新人创作礼」活动，一起开启掘金创作之路。

吴恩达机器学习作业7 - K-means和PCA主成分分析（Python实现）

Introduction

在本实验中，将实现K-means聚类算法，并将其应用于图像压缩。在第二部分实验中，将使用主成分分析来寻找人脸图像的低维表示。

1 K-means Clustering

在本节实验中，将实现K-means算法，并将其用于图像压缩。首先从一个2D数据集开始，直观地了解K-means算法是如何工作。然后使用K-means算法进行图像压缩，通过减少图像中出现的颜色的数量，最后只剩下那些在图像中最常见的颜色。

1.1 Implementing K-means

K-means算法可以自动将类似的数据实例聚类在一起。比如，给定一个训练集$\lbrace x^{(1)},...,x^{(m)}\rbrace(x^{(i)}\in R^n)$，希望将每个数据$x^{(i)}$分配给最接近的簇中心。K-means是一个迭代的、无监督的聚类算法，将类似的实例聚合成簇。 该算法通过随机选取每个簇的初始聚类中心开始，然后重复将实例数据分配给最近的簇，并重新计算该簇的聚类中心。

该算法内部重复执行两个步骤：

• 将每个训练数据$x^{(i)}$分配给其最近的簇心；
• 使用分配给簇心的点来重新计算每个簇心的平均值。

K-means算法最终收敛于簇心的最终均值，但是收敛值并不总是最优的，它取决于簇心的初始设置。因此，在实验中，K-means算法通常使用不同的随机初始化的簇心运行多次，然后从不同的随机初始化的簇心中选择代价函数值（失真）最低的一个。

1.1.1 Finding closest centroids

在K-means算法的分配簇的阶段，算法将每一个训练样本$x^{(i)}$分配给最接近的簇心。具体来讲，对于每一个样本$i$：

其中，$c^{(i)}$表示离样本$x^{(i)}$最近的簇心的下标，$u_j$表示第$j$个簇心的坐标。

编写函数findClosestCentroids，该函数可以读取数据矩阵X和簇心集内所有簇心的坐标，并输出一个一维数组idx（表示每个训练样本最近的簇心的下标）。

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

def findClosestCentroids(X,centroids):
    idx=[]
    max_dist=1000000#限制一下最大距离
    for i in range(len(X)):
        d=(X[i]-centroids)
        dist=d[:,0]**2+d[:,1]**2
        if(dist.min()<max_dist):
            index_i=np.argmin(dist)#给出dist最小值的下标
            idx.append(index_i)
    return np.array(idx)

接下来使用作业提供的例子，自定义centroids：$[[3, 3], [6, 2], [8, 5]]$，测试前$3$个样本分配的簇心应该是$[0,2,1]$（下标从$0$开始，而作业中的pdf给定的下标从$1$开始，所以值比练习中的值小一个）。

data=loadmat('./data/ex7data2.mat')
X=data['X']
centroids=np.array([[3, 3], [6, 2], [8, 5]])
idx=findClosestCentroids(X, centroids)
print(idx[0:3])

[0 2 1]

1.1.2 Computing centroid means

运行上述代码可以得到每个点分配的簇心下标。第二阶段：每个簇心重新计算分配到它的点的平均值。具体来讲，对于每个质心$k$：

其中$C_k$是分配给簇心$k$的样本集。比如，如果有两个样本$x^{(3)}$和$x^{(5)}$被分配给簇心$k=2$，那么就可以更新$u_2=\frac{x^{(3)}+x^{(5)}}{2}$

编写函数computeCentroids，完成上述的更新功能，可以考虑循环，也可以考虑向量化。（向量化效率会更高）

def computeCentroids(X,idx):
    centroids=[]
    for i in range(len(np.unique(idx))):
        u_k=X[idx==i].mean(axis=0)#求每列的平均值
        centroids.append(u_k)
    return np.array(centroids)

下面的输出符合实验中的预期值。

computeCentroids(X, idx)

array([[2.42830111, 3.15792418],
       [5.81350331, 2.63365645],
       [7.11938687, 3.6166844 ]])

1.2 K-means on example dataset

本节实验涉及实际运行K-means算法的一些迭代次数和可视化结果。为了运行算法，只需要在将样本分配给最近的簇心时重新计算簇的聚类中心。

先编写函数plotData用于绘制数据点。

def plotData(X,centroids,idx=None,iter=0):

    color = ['red', 'green', 'blue', 'darkorange', 'salmon', 'olivedrab',
              'maroon', 'navy', 'sienna', 'tomato', 'lightgray', 'gainsboro'
              'coral', 'aliceblue', 'dimgray', 'mintcream', 'mintcream']#颜色集合

    subx=[]#不同簇的样本点集合
    if(idx is not None):
        for i in range(centroids[0].shape[0]):
            x_i=X[idx==i]#最终每个点 归属的 簇心 下标
            subx.append(x_i)
    else:
        subx=[X]# X中所有元素在一个簇内
    plt.figure(figsize=(8,6))
    for i in range(len(subx)):
        point=subx[i]
        plt.scatter(point[:,0], point[:,1],c=color[i],label="Cluster %d"%i)
    plt.legend()
    plt.grid(True)
    plt.xlabel('x1')
    plt.xlabel('x2')
    plt.title("Iteration number %d"%iter)

    #绘制簇心的移动轨迹
    cx,cy=[],[]
    for centroid in centroids:
        cx.append(centroid[:,0])
        cy.append(centroid[:,1])
    plt.plot(cx,cy,'x--',color='black',markersize=7)
plotData(X, [centroids])

​
​

接下来编写函数runKmeans，了解K-means算法是如何工作的，当运行下一步时，K-means算法将生成一个可视化过程，在每次迭代中逐步完成算法操作。可视化K-means算法的每个步骤如何改变簇心和簇分配。

def runKmeans(X,centroids,max_iters=10):
    K=len(centroids)#K表示簇心个数
    centroids_all=[centroids]
    centroids_now=centroids
    for i in range(max_iters):
        idx=findClosestCentroids(X, centroids_now)
        centroids_now=computeCentroids(X, idx)
        centroids_all.append(centroids_now)
    return idx,centroids_all

idx,centroids_all=runKmeans(X, centroids, 10)
plotData(X,centroids_all,idx,10)

​
​

1.3 Random initialization（随机初始化）

在实验中，对簇心进行初始化的时候可以考虑从训练样本中随机选择样本作为簇心。

接下来编写kMeansInitCentroids函数，先随机排列这些样本的索引，然后根据结果选择最好的$K$个样本点作为簇心。

def kMeansInitCentroids(X,K):
    n,m=X.shape#n个样本，每个样本有m个特征值（坐标为m维）
    idx=np.random.choice(n,K)#从[0,n)中选取K个样本
    centroids=X[idx]
    return centroids

进行三次随机初始化，看下各自的效果：

for i in range(3):
    centroids=kMeansInitCentroids(X,3)
    idx,centroids_all=runKmeans(X,centroids,10)
    plotData(X, centroids_all,idx,10)

​
​

发现第$1$次的效果不理想，这是因为落入了局部最优，这是比较常见的问题，所以多次随机化参数进行运行，选取最优的即可。

1.4 Image compression with K-means

在本实验中，将把K-means应用于图像压缩。在一个$24$位$bit$表示颜色的图像中：每个像素点表示为三个$8$位无符号整数($0\sim255$)，指定了红、绿和蓝色的强度值，这种编码被称为RGB编码。该图像包含数千种颜色，在这部分的实验中，需要把颜色的数量减少到$16$种颜色。

这可以有效地压缩图像，只需要存储$16$种颜色的RGB值，而对于图中的每个像素点，只需要将该颜色的索引存储在该位置(只需要$4bit$就能表示$16$种颜色)。

接下来使用K-means算法选$16$种颜色，用于图像压缩，把原始图像的每个像素看作一个数据样本，然后利用K-means算法找分组最好的$16$种颜色。

1.4.1 K-means on pixels

读入图像文件bird small.png，得到一个三维矩阵A，它的前两个值表示像素位置，最后一个值表示红色、绿色或蓝色。比如，$A(50,33,3)$表示坐标为$(50,33)$的像素点的蓝色强度值，$A(50,33,1)$表示坐标为$(50,33)$的像素点的红色强度值（这里下标从$1$开始，代码中下标从$0$开始，注意区分）。接下来，利用该矩阵，创建像素颜色矩阵$m\times 3(m=16384=128\times 128)$，并调用K-means相关函数。

from skimage import io
A=io.imread('./data/bird_small.png')
A=A/255 #RGB的三个值[0,255]，将它们范围设置为[0,1]
io.imshow(A)
plt.show()

​
​

A.shape

(128, 128, 3)

接下来找到最合适的$16$种颜色来表示图像，调用函数findClosestCentroids为每一个像素点寻找最近的簇心。由于减少了对像素点颜色描述的位数，原始图像对$128\times 128$个像素点中的每一个的RGB颜色描述需要$24$位（$3\times 8=24$，$2^8=256$，而RGB的范围为$[0,255]$），总大小为$128\times 128\times 24=393216$位，而现在只需要$16$种颜色，每一种需要$24$位，图像本身每个像素点需要$4$位来表示选择了哪种颜色（$2^4=16$），总大小为$16\times 24+128\times 128\times 4=65920$

X=A.reshape(-1,3)#创建像素颜色矩阵m*3
K=16#共16种颜色
centroids=kMeansInitCentroids(X,K)#随机初始化得到K个簇心
idx,centroids_all=runKmeans(X, centroids,10)#centroids_all:(11*16)

img=np.zeros(X.shape)#像素点的颜色
centroids=centroids_all[-1]#得到最合适的16种颜色
for i in range(len(centroids)):
    img[idx==i]=centroids[i]
img=img.reshape(128,128,3)

fig,ax_array=plt.subplots(1,2,figsize=(8,6))
ax_array[0].imshow(A)
ax_array[1].imshow(img)
plt.show()

​
​

可以看到虽然对图像进行了压缩，但图像的主要特征仍然存在。不过也看到了一些压缩伪影。

下面使用scikit-learn来实现K-means算法对图像压缩。

A=io.imread('./data/bird_small.png')
A=A/255 #RGB的三个值[0,255]，将它们范围设置为[0,1]
io.imshow(A)
plt.show() 

​
​

from sklearn.cluster import KMeans

参数的意义：

• n_clusters:簇的个数，即想聚成几类
• init: 初始簇中心的获取方法
• n_init: 获取初始簇中心的更迭次数，为了弥补初始质心的影响，算法默认会初始10次，实现算法，然后返回最好的结果。
• max_iter: 最大迭代次数（因为kmeans算法的实现需要迭代）
• n_jobs: 并行设置

X=A.reshape(-1,3)#创建像素颜色矩阵m*3
K=16#共16种颜色
model=KMeans(n_clusters=K,n_init=100,n_jobs=-1)

model.fit(X)

D:\Python\lib\site-packages\sklearn\cluster\_kmeans.py:792: FutureWarning: 'n_jobs' was deprecated in version 0.23 and will be removed in 1.0 (renaming of 0.25).
  warnings.warn("'n_jobs' was deprecated in version 0.23 and will be"





KMeans(n_clusters=16, n_init=100, n_jobs=-1)

centroids=model.cluster_centers_
centroids.shape

(16, 3)

C=model.predict(X)
C.shape

(16384,)

img=centroids[C]#对C种16384个像素点选取它对应于centroids中的颜色
img.shape

(16384, 3)

img=img.reshape(128,128,3)
img.shape

(128, 128, 3)

fig,ax_array=plt.subplots(1,2,figsize=(8,6))
ax_array[0].imshow(A)
ax_array[1].imshow(img)
plt.show()

​
​

2 Principal Component Analysis（主成分分析）

在本节实验中，将使用主成分分析(PCA)来执行降维操作。先对一个2D数据集进行实验，直观地了解PCA的工作过程，然后在$5000$个人脸图像数据集上使用它。

2.1 Example Dataset

为了理解PCA是如何工作的，首先从一个二维数据集开始，它有一个变化较大的方向和一个变化较小的方向。在这节实验中，将看到使用PCA将数据从2D减少到1D时会发生什么。

data=loadmat("./data/ex7data1.mat")
X=data["X"]
X.shape

(50, 2)

plt.figure(figsize=(8,6))
plt.scatter(X[:,0], X[:,1],facecolors='none',edgecolors='blue')
plt.title("Example Dataset 1")
plt.show()

​
​

2.2 Implementing PCA

在本部分实验中，将实现PCA，它主要两个步骤组成：

• 计算数据的协方差矩阵
• 用SVD计算特征向量$(U_1,U_2,...,U_n)$

在使用PCA之前，首先将数据集中的每个特征减去它们的平均值来对数据进行归一化，并对每个维度进行缩放，使它们在相同的范围内。（标准化数据）

数据标准化后，运行PCA来计算主成分。编写函数pca：首先，计算数据的协方差矩阵$\sum=\frac{1}{m}X^TX$，其中$X$是数据矩阵（每个样本数据是以行的形式表示），$m$时数据样本数量，$\sum$是一个$n\times n$的矩阵，而不是求和。然后，对$\sum$执行SVD来计算主成分，$[U,S,V]=svd(Sigma)$，其中$U$包含了主成分，每一列就是数据要映射的向量，$S$为对角矩阵，为奇异值。

先编写函数featureNormalize用于标准化数据。

def featureNormalize(X):
    means=X.mean(axis=0)
    stds=X.std(axis=0,ddof=1)
    X_norm=(X-means)/stds
    return X_norm,means,stds

协方差矩阵$\sum=\frac{1}{m}X^TX$，而$X$中每行表示一个样本，而此时需要对特征（列）进行压缩，并且对协方差矩阵做SVD（左奇异矩阵用于行数的压缩；右奇异矩阵可以用于列数即特征维度的压缩，也就是我们的PCA降维），所以出口基为$V$。

def pca(X):
    sigma=(X.T@X)/len(X)
    U,S,V=np.linalg.svd(sigma)
    return U,S,V

先对数据进行标准化，然后对协方差矩阵进行$SVD$操作，然后绘制图形查看结果。

X_norm,means,stds=featureNormalize(X)
U,S,V=pca(X_norm)
U,S,V

(array([[-0.70710678, -0.70710678],
        [-0.70710678,  0.70710678]]),
 array([1.70081977, 0.25918023]),
 array([[-0.70710678, -0.70710678],
        [-0.70710678,  0.70710678]]))

plt.figure(figsize=(8,6))
plt.scatter(X[:,0], X[:,1],facecolors='none',edgecolors='blue')#绘制原数据点
plt.plot([means[0],means[0]+1.5*S[0]*U[0][0]],[means[1],means[1]+1.5*S[0]*U[0][1]],c='red',linewidth=1,label="First Principal Component")
plt.plot([means[0],means[0]+1.5*S[1]*U[1][0]],[means[1],means[1]+1.5*S[1]*U[1][1]],c='green',linewidth=1,label="econd Principal Component")
plt.grid()
plt.axis("equal")#将单个方格的长宽设置为一样，看起来垂直
plt.legend()
plt.show()

​
​

2.3 Dimensionality Reduction with PCA（使用PCA进行降维）

有了主成分（矩阵U），可以用它将原始数据投影到一个较低维的空间中。对于这个任务，可以实现一个计算投影并且仅选择顶部$K$个分量的函数，这样就有效地减少了维数，比如$x^{(i)}\rightarrow z^{(i)}$（将$2D$数据降为$1D$）。在本部分实验中，将使用PCA返回的特征向量，将数据集投影到一个一维空间中。

2.3.1 Projecting the data onto the principal components（将数据投射到主成分上）

编写函数projectData：传入数据集X，主成分U和最终的期望维数K，将X中的每个样本投影到主成分矩阵U的最主要的$K$个主成分上（实际上就是主成分矩阵U的前$K$列）。将第一个样本投影到第一个维度上，期望看到的结果约为$1.481$或$-1.481$。

def projectData(X,U,K):
    return X@U[:,:K]

Z=projectData(X_norm, U, 1)

Z[0]

array([1.48127391])

2.3.2 Reconstructing an approximation of the data（重建近似数据）

在将数据投影到低维空间后，也可以将数据投影回原始的高维空间来近似地恢复数据。接下来编写函数recoverData：将Z中的每个样本投影回原始空间，并返回X_rec中恢复的近似值。将第一个样本进行恢复，得到近似值约为$[-1.047 -1.047]$。

def recoverData(Z,U,K):
    return Z@U[:,:K].T

X_rec=recoverData(Z,U,1)

X_rec[0]

array([-1.04741883, -1.04741883])

2.3.3 Visualizing the projections（可视化投影）

接下来同时执行投影和近似重建，来观察投影如何影响数据。在下图中，用蓝色的圆圈表示原始的数据点，用红色的圆圈表示投影的数据点。投影有效地保留了U1方向上的信息。

plt.figure(figsize=(8,6))
plt.scatter(X_norm[:,0], X_norm[:,1],facecolors='none',edgecolors='blue',label="Original Data Points")
plt.scatter(X_rec[:,0], X_rec[:,1],facecolors='none',edgecolors='red',label="Projected Data Points")
plt.xlabel("x1 (Feature Normalized)")
plt.ylabel("x2 (Feature Normalized)")
plt.title("The normalized and projected data after PCA")
plt.grid()
plt.legend()
plt.axis("equal")

#绘制投影线
for x in range(len(X_norm)):
    plt.plot([X_norm[x][0],X_rec[x][0]],[X_norm[x][1],X_rec[x][1]],'--',color='black')

plt.show()

​
​

2.4 Face Image Dataset

在本部分实验中，将对人脸图像运用PCA并查看它如何在实验中进行降维。数据集ex7faces.mat中包含了变量X作为人脸图像，每一个都是$32\times 32$的灰色像素格，X的每一行表示一个人脸图像的特征（长度为$1024$）。接下来实现函数ex7_pca：加载并可视化前$100$个人脸图像。

data=loadmat("./data/ex7faces.mat")
X=data['X']
X.shape

(5000, 1024)

def ex7_pca(X,row,col,plt_title=None):
    fig,ax_array=plt.subplots(nrows=row,ncols=col,figsize=(8,8))
    for r in range(row):
        for c in range(col):
            ax_array[r][c].imshow(X[r*col+c].reshape(32,32).T,cmap='Greys_r')#X[r*col+c].reshape(32,32)运行出来是倒着的，所以转置一下
            ax_array[r][c].set_xticks([])
            ax_array[r][c].set_yticks([])
    plt.suptitle(plt_title,color='white',x=0.5, y=0.95, fontsize=16)

ex7_pca(X, 10, 10,plt_title="Faces dataset")

​
​

2.4.1 PCA on Faces

为了在人脸图像数据集上运行PCA，首先需要从矩阵X中减去每个特征的平均值来对数据集进行标准化。然后运行PCA，获得数据的主成分。U（每一行）中的每个主成分都是一个长度为$n$的向量(其中对于人脸数据集，$n=1024$)。可以通过将这些主成分重塑成一个$32\times 32$的矩阵，对应于原始数据集中的像素。

X_norm,means,stds=featureNormalize(X)
U,S,V=pca(X_norm)

U.shape,S.shape

((1024, 1024), (1024,))

接下来用最合适的$36$个主成分来绘制变化最大的$36$个人脸图像：

ex7_pca(U[:,:36].T,6,6,plt_title="Principal components on the face dataset")

​
​

2.4.2 Dimensionality Reduction

现在已经得到了人脸图像数据集的主成分，可以使用它来减少人脸图像数据集的维度。这样大大减小了维度，可以加快算法速度。

接下来，将每个人脸图像都用一个向量来描述$z^{(i)}\in R^{100}$（用前$100$个主成分进行表示）。要了解在降维过程中丢失了什么，可以仅使用投影的数据集来恢复数据。对数据进行近似恢复，然后使用ex7_pca函数将原始的和投影的人脸图像并排显示，这样可以观察到人脸图像的一般结构和外观被保留了下来，而细节却丢失了。

如果正在训练一个神经网络来执行识别人的操作（识别一个人脸图像，然后预测这个人的身份），可以使用仅为100维的特征，而不是原始像素。

Z=projectData(X_norm, U, K=36)
X_rec=recoverData(Z, U, K=36)

ex7_pca(X_rec, 10, 10)
ex7_pca(X, 10, 10)

​