14、降维，PCA算法

使用PCA将数据从2D减少到1D :

"""
这部分，你将运用PCA来实现降维。
您将首先通过一个2D数据集进行实验，以获得关于PCA如何工作的直观感受，然后在一个更大的图像数据集上使用它。
"""
from scipy.io import loadmat
import matplotlib.pyplot as plt
import numpy as np

# 加载二维数据
mat = loadmat('data\ex7data1.mat')
X = mat['X']
print(X.shape)  # (50, 2)
plt.scatter(X[:, 0], X[:, 1], facecolors='none', edgecolors='b')  # 蓝色空心散点图
# plt.show()

计算PCA主成分,并画出主成分：

def featureNormalize(X):
    '''
    标准化数据
    Args:
        X:数据集

    Returns:
        X_norm：标准化后的数据集
        means：均值from scipy.io import loadmat
import matplotlib.pyplot as plt
import numpy as np
        stds：方差
    '''
    means = X.mean(axis=0)  # 对各列求均值，返回 1* n 矩阵
    stds = X.std(axis=0, ddof=1)  # 对各列求无偏样本标准差（除以n-1）
    X_norm = (X - means) / stds  # 数据标准化
    return X_norm, means, stds


def pca(X):
    '''
    计算主成分
    Args:
        X: 训练集

    Returns:
        U,S,V: 主成分
    '''
    sigma = (X.T @ X) / len(X)  # 协方差矩阵
    U, S, V = np.linalg.svd(sigma)  # 用SVD计算主成分
    return U, S, V

X_norm, means, stds = featureNormalize(X) # 数据标准化
U, S, V = pca(X_norm) # 计算主成分
print(U)
'''
[[-0.70710678 -0.70710678]
 [-0.70710678  0.70710678]]
'''
print(U[:, 0])  # 从2D降到1D，所以只取U的第一列  [-0.70710678 -0.70710678]
print(S)  # [1.70081977 0.25918023]
# 画出主成分
plt.figure(figsize=(7, 5))
plt.scatter(X[:, 0], X[:, 1], facecolors='none', edgecolors='b')  # 数据散点图
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 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='g', linewidth=3, label='Second Principal Component')  # 第二个主成分
plt.grid()  # 显示网格
plt.axis('equal')  # 将横轴纵轴的定标系数设成相同值 ,即单位长度相同。不然两条线是不垂直的。
plt.legend()
plt.show()

将数据投影到PCA主成分上：

def projectData(X, U, K):
    '''
    将数据投影到主成分上
    Args:
        X:数据集
        U:投影矩阵
        K:降维后的维数

    Returns:
        Z:降维后的数据集
    '''
    Z = X @ U[:, :K]  # U[:, :K] 取U的前K列
    return Z


Z = projectData(X_norm, U, 1)  # 数据投影降维
# print(Z)
'''
[[ 1.48127391]
 ...
 [ 0.36423084]
 [-1.42814204]]
'''

数据恢复(降维后的坐标)：

def recoverData(Z, U, K):
    '''
    数据恢复(降维后的坐标)
    Args:
        Z:降维后的数据集
        U:投影矩阵
        K:降维的维数

    Returns:
        X_rec：恢复后的数据
    '''
    X_rec = Z @ U[:, :K].T
    return X_rec


X_rec = recoverData(Z, U, 1)  # 恢复后的数据
# print(X_rec)
'''
[[-1.04741883 -1.04741883]
 ...
 [-0.2575501  -0.2575501 ]
 [ 1.00984892  1.00984892]]
'''

可视化降维投影：

# 可视化投影
plt.figure(figsize=(7, 5))
plt.axis("equal")
plot = plt.scatter(X_norm[:, 0], X_norm[:, 1], s=30, facecolors='none',
                   edgecolors='b', label='Original Data Points')  # 可视化源数据
plot = plt.scatter(X_rec[:, 0], X_rec[:, 1], s=30, facecolors='none',
                   edgecolors='r', 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)
plt.grid(True)
# 画出每一个样本降维前后的直线距离长度，黑色虚线连接。
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--')
    # 输入第一项全是X坐标，第二项都是Y坐标
plt.legend()
plt.show()

将PCA运用到人脸图像：

# 载入图像数据集
mat = loadmat('data/ex7faces.mat')
X = mat['X']
print(X.shape)  # (5000, 1024)


# 将图片展示出来
def displayData(X, row, col):
    fig, axs = plt.subplots(row, col, figsize=(8, 8))  # 画一个row行col列的子图
    for r in range(row):
        for c in range(col):
            axs[r][c].imshow(X[r * col + c].reshape(32, 32).T, cmap='Greys_r')  # 将数据恢复成单通道灰度图片矩阵(32*32=1024)
            # 关闭坐标刻度
            axs[r][c].set_xticks([])
            axs[r][c].set_yticks([])


displayData(X, 10, 10)  # 先展示100个图片
# plt.show()

# 找出图像中的主成分，并可视化
X_norm, means, stds = featureNormalize(X)  # 数据集标准化，均值，方差
U, S, V = pca(X_norm)  # 计算主成分
print(U.shape, S.shape)  # (1024, 1024) (1024,)
displayData(U[:, :36].T, 6, 6)  # 画出主成分(36维=6*6)
# plt.show()

# PCA降维
z = projectData(X_norm, U, K=36)  # 将图片降到36维
X_rec = recoverData(z, U, K=36)  # 将降维后的数据恢复为图像矩阵
displayData(X_rec, 10, 10)#显示降维后的图像
plt.show()