主成分分析 （Principal Component Analysis）

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

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

ex4data1.mat -2D 仿真数据集
ex4data2.mat -LFW人脸数据集
评分标准如下：

要点1：实现PCA算法-----------------（20分）
要点2：降维仿真数据-----------------（20分）
要点3：重构仿真数据-----------------（20分）
要点4：降维人脸数据-----------------（20分）
要点5：重构人脸数据-----------------（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

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

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

Step 2: 计算样本的协方差矩阵: X X T mathbf{X}mathbf{X}^T XXT

Step 3: 对协方差矩阵 X X T mathbf{X}mathbf{X}^T XXT 做特征值分解（或奇异值分解）

Step 4: 取最大的 d d d 个特征值所对应的特征向量 w 1 , ⋯   , w d mathbf{w}_1,cdots,mathbf{w}_{d} w1​,⋯,wd​

输出: 投影矩阵 W = [ w 1 , ⋯   , w d ] mathbf{W}=[mathbf{w}_1,cdots,mathbf{w}_{d}] W=[w1​,⋯,wd​]

# ====================== 在这里填入代码 =======================
def PCA(X,d):
    """   
    输入
    ----------
    X : 尺寸为 (D, m)的矩阵，第i列为第i个样本，D为样本的维数，m为样本的个数。 
        
    d:  期望的维数
    
    输出
    -------
    W : 尺寸为 (D, d)的投影矩阵，第i列为协方差矩阵的第i个特征向量。 
    """
    
    X=[(x-np.mean(x))/np.std(x,axis=0,ddof=1) for x in X]
    matrix=np.cov(X[0],X[1])
    X=np.array(X)
    Z=np.matmul(X,X.T)
    U,S,V=np.linalg.svd(Z)
    W=U[:,:d]
    return W 
# ============================================================= 

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

# Load the dataset into the variable X 
data = loadmat(os.path.join('ex4data1.mat'))
X = data['X']

#  Run PCA
W= PCA(X.T,2)


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

第一个特征向量: W[:, 0] = [-0.707107, -0.707107]

在本部分实验中，将已实现的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)

#  可视化标准化之后的数据集
mu = np.mean(X,axis=0)
sigma = np.std(X,axis=0)
X_norm = (X - mu)/sigma
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)

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

# ====================== 在这里填入代码 =======================
def ProjectData(X, W):
    """   
    输入
    ----------
    X : 尺寸为 (D, m)的矩阵，第i列为第i个样本，D为样本的维数，m为样本的个数。 
    W : 尺寸为 (D, d)的投影矩阵，第i列为协方差矩阵的第i个特征向量。   
    
    输出
    -------
    Z : 尺寸为 (d, m)的矩阵，第i列为第i个降维后的数据。 
    """
    X=[(x-np.mean(x))/np.std(x,axis=0,ddof=1) for x in X]
    Z=np.matmul(W.T,X)
    

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

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

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

第一个样本降维后的数据: 1.481274

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

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


    return X_rec
# ============================================================= 

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

#  重构原始数据
X_rec = RecoverData(Z, W).T

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

第一个重构后的数据: X_rec[0, :] = [-1.047419, -1.047419]

mu = np.mean(X,axis=0)
sigma = np.std(X,axis=0, ddof=1)
X_norm = (X - mu)/sigma

#  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, X_rec):
    ax.plot([xnorm[0], xrec[0]], [xnorm[1], xrec[1]], '--k', lw=1)

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

#  Load Face dataset
data = loadmat(os.path.join( 'ex4data2.mat'))
X = data['X']
#定义显示人脸图像函数
def displayData(X, example_width=None, figsize=(8, 8)):

    """   
    输入
    ----------
    X : 尺寸为 (m, D)的矩阵，第i列为第i个样本，D为样本的维数，m为样本的个数。   
    """
    
    # Compute rows, cols
    if X.ndim == 2:
        m, n = X.shape
    elif X.ndim == 1:
        n = X.size
        m = 1
        X = X[None]  # Promote to a 2 dimensional array
    else:
        raise IndexError('Input X should be 1 or 2 dimensional.')

    example_width = example_width or int(np.round(np.sqrt(n)))
    example_height = int(n / example_width)

    # Compute number of items to display
    display_rows = int(np.floor(np.sqrt(m)))
    display_cols = int(np.ceil(m / display_rows))

    fig, ax_array = plt.subplots(display_rows, display_cols, figsize=figsize)
    fig.subplots_adjust(wspace=0.025, hspace=0.025)

    ax_array = [ax_array] if m == 1 else ax_array.ravel()

    for i, ax in enumerate(ax_array):
        ax.imshow(X[i].reshape(example_height, example_width, order='F'), cmap='gray')
        ax.axis('off')
 #  Display the first 100 faces in the dataset
displayData(X[:64, :], figsize=(8, 8))

# ====================== 在这里填入代码 =======================
W1=PCA(X.T,64)
print(W1[:36])
Z1=ProjectData(X.T,W1)



# ============================================================= 
[[-0.01425307 -0.03606596 -0.04561884 ...  0.0061877   0.02367001
   0.05076207]
 [-0.01474233 -0.03809858 -0.04756249 ...  0.00101133  0.01959809
   0.05250727]
 [-0.01501482 -0.03988747 -0.05082085 ... -0.00338655  0.00776122
   0.05543446]
 ...
 [-0.01711169 -0.04114578 -0.04517263 ...  0.0064823   0.0409063
   0.02625426]
 [-0.01721429 -0.04300162 -0.04873793 ...  0.00549935  0.02761684
   0.03023218]
 [-0.01742302 -0.04464136 -0.05105597 ...  0.01009193  0.00885712
   0.03474425]]

# ====================== 在这里填入代码 =======================
X1_rec=RecoverData(Z1,W1).T

displayData(X1_rec[:64,:],figsize=(8,8))
# =============================================================