机器学习|吴恩达机器学习之PCA

发表于

K-means Clustering

导入模块

1
2
3
4
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
from skimage import io

Implementing K-means

Finding closet centroids

  • 读入数据
1
2
data=loadmat('ex7data2.mat')
X=data['X']
  • 为每个样本找到最近的聚类中心

$c^{(i)} :=j \quad$ that minimizes $\quad\left|x^{(i)}-\mu_{j}\right|^{2}$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#========================找聚类中心
def find_closet_centroids(X,centroids):
    '''
    返回每个样本所在的cluster索引
    '''
    idx=[]
    max_dist=10000 # 给一个距离限定
    for i in range(len(X)):
        minus=X[i]-centroids # minus是3x2的矩阵，每一行代表了第i个样本到一个centroids的x1,x2距离
        dist=minus[:,0]**2+minus[:,1]**2 #求范式，即直线距离,dist是3x1的向量
        if dist.min()<max_dist:
            ci=np.argmin(dist) #返回沿axis的最小值索引
            idx.append(ci)
    return np.array(idx)

可以自己初始化一组聚类中心来测试一下

1
2
3
init_centroids = np.array([[3, 3], [6, 2], [8, 5]])
idx = findClosestCentroids(X, init_centroids)
print(idx[0:3])

结果应该是

1
[0 2 1]

Computing centroids means

1
2
3
4
5
6
7
8
9
10
#=========================移动聚类中心
def compute_centroids(X,idx):
    centroids=[]
    for i in range(len(np.unique(idx))):
        # 布尔索引，idx==i运算返回bool value，根据值为true的下标来输出X中对应元素值
        u_k=X[idx==i].mean(axis=0)
        
        centroids.append(u_k)
    
    return np.array(centroids)

K-means on example dataset

  • 画图函数
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
def plot_data(X,centroids,idx=None):

    colors = ['b','g','gold','darkorange','salmon','olivedrab', 
              'maroon', 'navy', 'sienna', 'tomato', 'lightgray', 'gainsboro'
             'coral', 'aliceblue', 'dimgray', 'mintcream', 'mintcream']
	
    # 这里centroids需要从ndarray转成list,但centroids[0]仍为ndarray.
    assert len(centroids[0])<=len(colors),'colors not enough' 
    
    subX=[]
    if idx is not None:
        for i in range(centroids[0].shape[0]): #分成几类就循环几次
            x_i=X[idx==i]
            subX.append(x_i) #把数据按cluster分成不同下标的元素，存储在subx这个list中
    else:
        subX=[X]
    
    
    fig=plt.figure(figsize=(8,5))
    for i in range(len(subX)):
        xx=subX[i]
        plt.scatter(xx[:,0],xx[:,1],c=colors[i])    
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Plot of X Points')

    xx,yy=[],[]
    for centroid in centroids: 
        xx.append(centroid[:,0])
        yy.append(centroid[:,1])
    
    plt.plot(xx,yy,'rx--')

    plt.show()
  • 画出聚类中心移动过程
1
2
3
4
5
6
7
8
9
10
11
#=============================聚类中心移动过程
def run_kmeans(X,centroids,max_iters):

    centroids_all=[]
    centroids_all.append(centroids)
    centroid_i=centroids
    for i in range(max_iters):
        idx=find_closet_centroids(X,centroid_i)
        centroid_i=compute_centroids(X,idx) 
        centroids_all.append(centroid_i) #每次移动后的聚类中心坐标都记录下来
    return idx,centroids_all

结果：

EYOrCQ.png

Random initialization

关于聚类中心的初始化，一个更好的方法是从样本集中随机选取。

1
2
3
4
5
6
#==============================随机初始化聚类中心
def random_centroids(X,K):
    m=X.shape[0]
    index=np.random.choice(m,K) #在0~m中随机生成K个样本
    
    return X[index]

最后再尝试一下聚类算法：

1
2
3
4
5
for i in range(3):
        init_centroids=random_centroids(X,3)
        
        idx,centroids_all=run_kmeans(X,init_centroids,20)
        plot_data(X,centroids_all,idx=idx) #idx为每个样本所在cluster的索引

得到最终结果：

EYXtG4.md.png

Image compression with K-means

  • 导入数据
1
2
3
4
img=io.imread('bird_small.png')
plt.imshow(img)
plt.show()
print(img.shape)

图像为

EYX5eP.png

1
(128,128,3)

可以看到图像以3维矩阵的方式存储，前面两个维度代表图像的像素点个数128x128，最后一个维度代表像素点由RGB三个通道表示。又因为一个通道占用8-bit，因此在原始的图像中一个像素点需要24-bit来储存。

在这幅图中包含了上千种种颜色，而在这个实验中，我们把颜色降为16种，也就是说，一像素点只需要4bit就足够了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
img=img/255
X=img.reshape(-1,3) #转换成128x128行，3列为RGB三通道
K=16 #16个聚类中心，就是把所有颜色压缩为16种RGB颜色，那么每个像素值需要4bit存储即可
init_centroids=random_centroids(X,K)
idx,centroids_all=run_kmeans(X,init_centroids,10)
img_2=np.zeros(X.shape)
centroids=centroids_all[-1] # 只需要记录聚类中心最后移动的位置即可

for i in range(len(centroids)):
    img_2[idx==i]=centroids[i]


img_2=img_2.reshape((128,128,3))
fig,axes=plt.subplots(1,2,figsize=(12,6))
axes[0].imshow(img)
axes[1].imshow(img_2)
plt.show()

比对一下压缩效果：

EYjBlj.md.png

Principle Component Analysis

  • 导入模块
1
2
3
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat

Example Dataset

  • 导入数据
1
2
3
4
data = loadmat('ex7data1.mat')
X = data['X']
print(X.shape)
plt.scatter(X[:,0], X[:,1], facecolors='none', edgecolors='b')
1
(50,2)
  • 显示图像：

EYz5cQ.png

Implementing PCA

关于PCA的数学原理可以参考：浅谈SVD的数学原理及应用

PCA主要分为两个计算步骤：

  • 计算数据的协方差矩阵：$\Sigma=\frac{1}{m} \sum_{i=1}^{n}\left(x^{(i)}\right)\left(x^{(i)}\right)^{T}$
  • 用SVD函数进行奇异值分解

在开始上面的步骤之前，需要先对数据进行特征缩放和归一化：

1
2
3
4
5
6
def feature_normalize(X):
    means = X.mean(axis=0)
    stds = X.std(axis=0, ddof=1)
    X_norm = (X - means) / stds

    return X_norm, means

PCA：

1
2
3
4
5
6
#===================================PCA
def pca(X):
    sigma = (1 / len(X)) * (X.T @ X)  #求出协方差矩阵
    U, S, V = np.linalg.svd(sigma)  #奇异值分解

    return U, S, V

可视化主成成分：

1
2
3
4
5
6
7
8
9
10
11
12
13
def plot_reduce(means, U, S):
    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='r',
             linewidth=3,
             label='First Principle 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='g',
             linewidth=3,
             label='Second Principal Component')
    plt.axis('equal')
    plt.legend()

得到结果：

EtpMM4.md.png

Dimensionality Reduction with PCA

Projecting the data onto the principal components

1
2
3
4
5
6
7
'''
返回值
	- Z：投影到主成成分后的样本点
'''
def project_data(X, U, K):
    Z = X @ U[:, 0:K]
    return Z
1
2
3
K = 1
Z = project_data(X_norm, U, K)
print(Z[0])

得到结果：

1
array([1.48127391])

Reconstructing an approximation of the data

1
2
3
4
5
#===================================重建数据
def recover_data(Z, U, K):
    X_rec = Z @ U[:, 0:K].T

    return X_rec
1
2
3
4
# you will recover an approximation of the first example and you should see a value of
# about [-1.047 -1.047].
X_rec = recover_data(Z, U, 1)
X_rec[0]
1
array([-1.04741883, -1.04741883])

Visualizing the projections

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#====================================PCA样本投影可视化
def visualize_data(X_norm, X_rec):
    fig, ax = plt.subplots(figsize=(8, 5))
    plt.axis('equal')
    plt.scatter(
        X_norm[:, 0],
        X_norm[:, 1],
        s=30,
        facecolors='none',
        edgecolors='blue',
        label='Original Data Points')

    plt.scatter(
        X_rec[:, 0],
        X_rec[:, 1],
        s=30,
        facecolors='none',
        edgecolors='red',
        label='PCA Reduced Data Points')

    plt.title("Example Dataset: Reduced Dimension Points Shown",fontsize=14)
    plt.xlabel('x1 [Feature Normalized]',fontsize=14)
    plt.ylabel('x2 [Feature Normalized]',fontsize=14)

    for x in range(X_norm.shape[0]):
        plt.plot([X_norm[x,0],X_rec[x,0]],[X_norm[x,1],X_rec[x,1]],'k--')
        plt.legend(loc=0)
        # 输入第一项全是X坐标，第二项都是Y坐标
    plt.show()

结果如图：

Et9DcF.md.png

Face Image Dataset

  • 导入数据
1
2
data_2=loadmat('ex7faces.mat')
X_2=data_2['X']
  • 数据可视化
1
2
3
4
5
6
7
8
def plot_face(X,row,col):
    fig,ax=plt.subplots(row,col,figsize=(8,8))
    for i in range(row):
        for j in range(col):
            ax[i][j].imshow(X[i*col+j].reshape(32,32).T,cmap='Greys_r')
            ax[i][j].set_xticks([])
            ax[i][j].set_yticks([])
    plt.show()

结果如图：

Et9Ije.png

PCA on Faces

1
2
X_2_norm,means_2=feature_normalize(X_2)
U_2,S_2,V_2=pca(X_2_norm)

显示主成成分

1
plot_face(U[:,:36].T, 6, 6)

结果如下：

Et9O4P.png

Dimensionality Reduction

把图片降维

1
2
3
Z_2=project_data(X_2_norm,U_2,K=36)
X_2_rec=recover_data(Z_2,U_2,K=36)
plot_face(X_2_rec,10,10)

结果：

Etn09H.png

-------------End-------------
梦想总是要有的，万一有人有钱呢？