机器学习：主成分分析 （Principal Component Analysis）

介绍（Introduction）

在本次实验中，将实现主成分分析方法，并使用它获得人脸图像的低维表示。

本次实验需要用到的数据集包括：

• exdata1.mat -2D 仿真数据集
• exdata2.mat -LFW人脸数据集

评分标准如下：

• 要点1：数据标准化--------------------（10分）
• 要点2：实现PCA算法-----------------（30分）
• 要点3：降维仿真数据-----------------（10分）
• 要点4：重构仿真数据-----------------（10分）
• 要点5：降维人脸数据-----------------（20分）
• 要点6：重构人脸数据-----------------（20分）

# 引入所需要的库文件
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sb
from scipy.io import loadmat
import scipy

%matplotlib inline

1 实现主成分分析

在本部分实验中，将实现主成分分析算法。 算法步骤如下：

Step 1: 对所有样本进行标准化使得样本的均值为0，标准差为1

Step 2: 计算样本的协方差矩阵: ，其中N为样本个数

Step 3: 对协方差矩阵   做奇异值分解（或特征值分解）

Step 4: 取最大的 d� 个特征值所对应的特征向量 ,⋯,

输出: 投影矩阵 W=[,⋯,]

**要点 1：** 在下方cell中，请**实现''数据标准化''的代码**。

# ====================== 在这里填入代码 =======================
def Normalize(X):
    """   
    输入
    ----------
    X : 尺寸为 (D, N)的矩阵，第i列为第i个样本，D为样本的维数，N为样本的个数。 
    
    输出
    -------
    X_norm : 尺寸为 (D, N)的矩阵。 
    """  
    m = np.mean(X, axis=1)
    sigma = np.std(X, axis=1)
    X_norm = (X - m[:, None]) / sigma[:, None]
    return X_norm 
# ============================================================= 

**要点 2：** 在下方cell中，请**实现''主成分分析''的代码**。

# ====================== 在这里填入代码 =======================
def PCA(X,d):
    """   
    输入
    ----------
    X : 尺寸为 (D, N)的矩阵，第i列为第i个样本，D为样本的维数，N为样本的个数。 
        
    d:  期望的维数
    
    输出
    -------
    W : 尺寸为 (D, d)的投影矩阵，第i列为协方差矩阵的第i个特征向量。 
    """
    
    # Step 1: 标准化（调用要点1中的代码）
    X_norm = Normalize(X)
    
    # Step 2: 计算协方差矩阵
    cov_matrix = np.dot(X_norm, X_norm.T) / X.shape[1]
    
    # Step 3: 奇异值分解（查询numpy中进行奇异值分解的函数）
    U, S, V = np.linalg.svd(cov_matrix)
    
    # Step 4: 取最大的𝑑个特征值所对应的特征向量并组成矩阵W
    W = U[:, :d] 
    
    return W 
# ============================================================= 

如果完成了上述函数 pca，以下代码可用于测试。如果结果为[-0.707107, -0.707107]，则计算通过。

# Load the dataset into the variable X 
data = loadmat('ex4data1.mat')
X = data['X']
X=X.T
#  Run PCA
W= PCA(X,2)

print('第一个特征向量: W[:, 0] = [{:.6f}, {:.6f}]'.format(W[0, 0], W[1, 0]))

2 将PCA应用于仿真数据

在本部分实验中，将已实现的PCA算法应用于数据集1，该数据集中的样本维数为2，因此降维结束后，可通过可视化观察降维前后的结果。

#  可视化原始数据集
plt.plot(X[0, :], X[1, :], 'bo', ms=10, mec='k', mew=1)
plt.axis([0.5, 6.5, 2, 8])
plt.gca().set_aspect('equal')
plt.grid(False)

#  可视化标准化之后的数据集
X_norm=Normalize(X)
plt.plot(X_norm[0, :], X_norm[1, :], 'bo', ms=10, mec='k', mew=1)
plt.axis([-3, 2.75, -3, 2.75])
plt.gca().set_aspect('equal')
plt.grid(False)

2.1 利用PCA降维

利用上述PCA代码可将仿真数据降至1维。

**要点 3：** 在下方cell中，请**实现''ProjectData''的代码**，输出降维后的数据。

# ====================== 在这里填入代码 =======================
def ProjectData(X, W):
    """   
    输入
    ----------
    X : 尺寸为 (D, N)的矩阵，第i列为第i个样本，D为样本的维数，N为样本的个数。 
    W : 尺寸为 (D, d)的投影矩阵，第i列为协方差矩阵的第i个特征向量。   
    
    输出
    -------
    Z : 尺寸为 (d, N)的矩阵，第i列为第i个降维后的数据。 
    """

    # Step 1: 标准化（调用要点1中的代码）
    X_norm = Normalize(X)
    
    # Step 2: 计算降维后的数据
    Z = np.dot(W.T,X_norm)

    return Z
# ============================================================= 

如果完成了上述函数 ProjectData，以下代码可用于测试。如果结果为 1.496313，则计算通过。

#  将数据降至一维 
W = PCA(X,1)
Z = ProjectData(X, W)
print('第一个样本降维后的数据: {:.6f}'.format(Z[0, 0]))

2.2 通过降维后的数据重构原始数据

利用降维后的数据矩阵Z�，可近似重构原始数据。

**要点 4：** 在下方cell中，请**实现''RecoverData''的代码**，输出重构后的数据。

# ====================== 在这里填入代码 =======================
def RecoverData(Z, W):
    """   
    输入
    ----------
    Z : 尺寸为 (d, N)的矩阵，第i列为第i个降维后的数据。 
    W : 尺寸为 (D, d)的投影矩阵，第i列为协方差矩阵的第i个特征向量。   
    
    输出
    -------
    X_rec : 尺寸为 (D, N)的矩阵，第i列为第i个重构后的数据。 
    """
    X_rec = np.dot(W,Z)
    return X_rec
# ============================================================= 

如果完成了上述函数 RecoverData，以下代码可用于测试。如果结果为[-1.058053, -1.058053]，则计算通过。

#  重构原始数据
X_rec = RecoverData(Z, W)
# print('第一个重构后的数据: {:.6f}'.format(X_rec[:, 0]))
print('第一个重构后的数据: X_rec[:, 0] = [{:.6f}, {:.6f}]'.format(X_rec[0, 0], X_rec[1, 0]))



#可视化重构前后的结果
X_norm=Normalize(X) 

#  Plot the normalized dataset (returned from featureNormalize)
fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(X_norm[0, :], X_norm[1, :], 'bo', ms=8, mec='b', mew=0.5)
ax.set_aspect('equal')
ax.grid(False)
plt.axis([-3, 2.75, -3, 2.75])

# Draw lines connecting the projected points to the original points
ax.plot(X_rec[0, :], X_rec[1, :], 'ro', mec='r', mew=2, mfc='none')
for xnorm, xrec in zip(X_norm.T, X_rec.T):
    ax.plot([xnorm[0], xrec[0]], [xnorm[1], xrec[1]], '--k', lw=1)

3 将PCA应用于人脸重构

在本部分实验中，将已实现的PCA算法应用于数据集2，为LFW人脸数据集合的子集。

#  Load Face dataset
data = loadmat('ex4data2.mat')
X = data['X']
X=X.T
X = X[:,:100]


#定义显示人脸图像函数
def displayImage(X):
    """   
    输入
    ----------
    X : 尺寸为 (D, N)的矩阵，第i列为第i个样本，D为样本的维数，N为样本的个数。   
    """
    # Compute number of items to display
    N=X.shape[1]
    display_rows = int(np.floor(np.sqrt(N)))
    display_cols = int(np.ceil(N / display_rows))
    
    figsize=(N,N)
    fig, ax_array = plt.subplots(display_rows, display_cols)

    ax_array = [ax_array] if N == 1 else ax_array.ravel()
    for i, ax in enumerate(ax_array):
        ax.imshow(X[:,i].reshape(32, 32, order='F'), cmap='gray')
        ax.axis('off')

3.1 利用PCA对人脸图像降维

**要点 5：** 利用PCA代码将人脸数据降至64维，并显示前16个特征向量。

# ====================== 在这里填入代码 =======================

#Step 1: 利用PCA对原始数据降维
W=PCA(X,64)

#Step 2: 调用前文中的displayImage函数显示前16个特征向量所对应的图像
displayImage(W[:, :16])
plt.gcf().suptitle('Top 16 eigenvectors')
plt.show()
# ============================================================= 

3.2 通过主成分重构原始人脸数据

**要点 6：** 利用降维后的数据重构原始人脸数据，并显示前16个重构前后的人脸数据。

# ====================== 在这里填入代码 =======================

#Step 1: 利用RecoverData函数得到重构数据 
Z = ProjectData(X,W)
X_rec = RecoverData(Z,W)

#Step 2: 利用displayImage函数显示前16个原始图像  

displayImage(X[:,:16])
plt.gcf().suptitle('Original faces')

#Step 3: 利用displayImage函数显示前16个重构图像 

displayImage(X_rec[:,:16])
plt.gcf().suptitle('Recovered faces')
plt.show()
# =============================================================