第16章 主成分分析：四个案例及课后习题

1.假设$x$为$m$ 维随机变量，其均值为$\mu$，协方差矩阵为$\Sigma$。

考虑由$m$维随机变量$x$到$m$维随机变量$y$的线性变换
$$y_i={\alpha }_i^{T}x=\sum _{k=1}^m{\alpha }_{ki}{x}_k,i=1,2,\cdots ,m$$

其中${\alpha }_i^T=({\alpha }_{1i},{\alpha }_{2i},\cdots ,{\alpha }_{mi})$。

如果该线性变换满足以下条件，则称之为总体主成分：

（1）${\alpha }_i^T{\alpha }_i=1,i=1,2,\cdots ,m$；

（2）$\mathrm{cov}(y_i,y_j)=0(i\ne j)$;

（3）变量$y_1$是$x$的所有线性变换中方差最大的；$y_2$是与$y_1$不相关的$x$的所有线性变换中方差最大的；一般地，$y_i$是与$y_1,y_2,\cdots ,y_{i-1},(i=1,2,\cdots ,m)$都不相关的$x$的所有线性变换中方差最大的；这时分别称$y_1,y_2,\cdots ,y_m$为$x$的第一主成分、第二主成分、…、第$m$主成分。

2.假设$x$是$m$维随机变量，其协方差矩阵是$\Sigma$，$\Sigma$的特征值分别是${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_m\ge 0$，特征值对应的单位特征向量分别是${\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_m$，则$x$的第2主成分可以写作

$$y_i={\alpha }_i^{T}x=\sum _{k=1}^m{\alpha }_{ki}{x}_k,i=1,2,\cdots ,m$$
并且，$x$的第$i$主成分的方差是协方差矩阵$\Sigma$的第$i$个特征值，即$$\mathrm{var}(y_i)={\alpha }_i^T\Sigma {\alpha }_i={\lambda }_i$$

3.主成分有以下性质：

主成分$y$的协方差矩阵是对角矩阵$$\mathrm{cov}(y)=\Lambda =\mathrm{diag}({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_m)$$

主成分$y$的方差之和等于随机变量$x$的方差之和
$$\sum _{i=1}^m{\lambda }_i=\sum _{i=1}^m{\sigma }_{ii}$$
其中${\sigma }_{ii}$是$x_2$的方差，即协方差矩阵$\Sigma$的对角线元素。

主成分$y_k$与变量$x_2$的相关系数$\rho (y_k,x_i)$称为因子负荷量（factor loading），它表示第$k$个主成分$y_k$与变量$x$的相关关系，即$y_k$对$x$的贡献程度。
$$\rho (y_k,x_i)=\frac{\sqrt{{\lambda }_k}{\alpha }_{ik}}{\sqrt{{\sigma }_{ii}}},k,i=1,2,\cdots ,m$$

4.样本主成分分析就是基于样本协方差矩阵的主成分分析。

给定样本矩阵
$$X=\left[\begin{array}{llll}{x}_1& x_2& \cdots & x_n\end{array}\right]=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1n}\\ x_{21}& x_{22}& \cdots & x_{2n}\\ ⋮& ⋮& & ⋮\\ x_{m1}& x_{m2}& \cdots & x_{mn}\end{array}\right]$$

其中$x_j=(x_{1j},x_{2j},\cdots ,x_{mj}{)}^T$是$x$的第$j$个独立观测样本，$j=1,2，\dots ,n$。

$X$的样本协方差矩阵
$$\begin{array}{c}S=[s_{ij}{]}_{m\times m},s_{ij}=\frac{1}{n-1}{\sum }_{k=1}^n(x_{ik}-{\overline{x}}_i)(x_{jk}-{\overline{x}}_j)\\ i=1,2,\cdots ,m,j=1,2,\cdots ,m\end{array}$$

给定样本数据矩阵$X$，考虑向量$x$到$y$的线性变换$$y=A^{T}x$$
这里
$$A=\left[\begin{array}{llll}{a}_1& a_2& \cdots & a_m\end{array}\right]=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1m}\\ a_{21}& a_{22}& \cdots & a_{2m}\\ ⋮& ⋮& & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mm}\end{array}\right]$$

如果该线性变换满足以下条件，则称之为样本主成分。样本第一主成分$y_1=a_1^{T}x$是在$a_1^T{a}_1=1$条件下，使得$a_1^T{x}_j(j=1,2,\cdots ,n)$的样本方差$a_1^{T}S{a}_1$最大的$x$的线性变换；

样本第二主成分$y_2=a_2^{T}x$是在$a_2^T{a}_2=1$和$a_2^T{x}_j$与$a_1^T{x}_j(j=1,2,\cdots ,n)$的样本协方差$a_1^{T}S{a}_2=0$条件下，使得$a_2^T{x}_j(j=1,2,\cdots ,n)$的样本方差$a_2^{T}S{a}_2$最大的$x$的线性变换；

一般地，样本第$i$主成分$y_i=a_i^{T}x$是在$a_i^T{a}_i=1$和$a_i^T{x}_j$与$a_k^T{x}_j(k<i,j=1,2,\cdots ,n)$的样本协方差$a_k^{T}S{a}_i=0$条件下，使得$a_i^T{x}_j(j=1,2,\cdots ,n)$的样本方差$a_k^{T}S{a}_i$最大的$x$的线性变换。

5.主成分分析方法主要有两种，可以通过相关矩阵的特征值分解或样本矩阵的奇异值分解进行。

（1）相关矩阵的特征值分解算法。针对$m\times n$样本矩阵$X$，求样本相关矩阵
$$R=\frac{1}{n-1}X{X}^T$$
再求样本相关矩阵的$k$个特征值和对应的单位特征向量，构造正交矩阵
$$V=(v_1,v_2,\cdots ,v_k)$$

$V$的每一列对应一个主成分，得到$k\times n$样本主成分矩阵
$$Y=V^{T}X$$

（2）矩阵$X$的奇异值分解算法。针对$m\times n$样本矩阵$X$
$$X^{\prime }=\frac{1}{\sqrt{n-1}}{X}^T$$
对矩阵$X^{\prime }$进行截断奇异值分解，保留$k$个奇异值、奇异向量，得到
$$X^{\prime }=US{V}^T$$
$V$的每一列对应一个主成分，得到$k\times n$样本主成分矩阵$Y$
$$Y=V^{T}X$$

案例1：高维数据的可视化

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
from matplotlib_inline import backend_inline

backend_inline.set_matplotlib_formats('svg')

iris = load_iris()
X = iris.data  # 形状(150,4)
y = iris.target

pca = PCA(n_components=2)
X_dr = pca.fit_transform(X)
# print(X_dr) # 形状(150,2)

colors = ['red', 'black', 'orange']
plt.figure()
for i in range(3):
    plt.scatter(X_dr[y == i, 0],
                X_dr[y == i, 1],
                c=colors[i],
                alpha=0.7,
                label=iris.target_names[i])
plt.legend()
plt.title('PCA of IRIS dataset')
plt.show()

# 一些关键信息
print('各主成分的方差（即协方差矩阵的特征值）：', pca.explained_variance_, '\n特征向量:',
      pca.components_, '\n可解释性方差贡献率：', pca.explained_variance_ratio_,
      '\n保留了原数据多少信息（累计方差贡献率）：', pca.explained_variance_ratio_.sum())

各主成分的方差（即协方差矩阵的特征值）： [4.22824171 0.24267075] 
特征向量: [[ 0.36138659 -0.08452251  0.85667061  0.3582892 ]
 [ 0.65658877  0.73016143 -0.17337266 -0.07548102]] 
可解释性方差贡献率： [0.92461872 0.05306648] 
保留了原数据多少信息（累计方差贡献率）： 0.977685206318795

接下来说明选择主成分个数（n_components）的一些方法

# 法一：选择累计方差贡献率曲线的突然变平坦的点对应的主成分个数
pca_ = PCA(n_components=4).fit(X)
plt.plot(range(1, 5), np.cumsum(pca_.explained_variance_ratio_))
plt.xticks([1, 2, 3, 4])
plt.show()  # 故选择主成分个数为2较合适

# 法二：令n_components='mle'使用极大似然估计 自动选择主成分个数
pca_mle = PCA(n_components='mle').fit(X)
X_mle = pca_mle.fit_transform(X)
print(X_mle.shape, pca_mle.explained_variance_ratio_.sum())
# 可以看出这种方法选择主成分的个数为3

(150, 3) 0.9947878161267247

# 法三：输出浮点数n_components=0.9按信息量占比选择主成分个数
# 确保选择的主成分能够解释至少 97% 的数据方差，并使用 完全SVD方法 进行计算
pca_3 = PCA(n_components=0.97, svd_solver="full").fit(X)
X_3 = pca_3.fit_transform(X)
print(X_3.shape, pca_3.explained_variance_ratio_.sum()) # 选择两个主成分个数就能包含原数据信息97%以上了

(150, 2) 0.977685206318795

svd_solver 参数决定了计算主成分的奇异值分解 (SVD) 方法。以下是 svd_solver 参数的可选值及其含义：

• auto：默认值。根据输入数据的大小和所需的主成分数量，自动选择合适的 SVD 方法。
• full：使用标准的完全 SVD 方法。这种方法比较慢，但适用于所有情况，特别是当数据矩阵的形状为 n_samples > n_features 或 n_features > n_samples 时。
• arpack：使用 scipy.sparse.linalg.svds 实现的截断 SVD。这种方法适合计算少数几个主成分的情况，特别是当所需的主成分数远小于数据矩阵的形状时。
• randomized：使用随机 SVD 方法，适用于大数据集，且需要计算的主成分数不多的情况。这种方法通常比完全 SVD 方法快很多。

案例2：人脸识别

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import fetch_lfw_people
from matplotlib_inline import backend_inline

# 设置 matplotlib 的显示格式为 SVG 矢量图
backend_inline.set_matplotlib_formats('svg')

# 获取 LFW (Labeled Faces in the Wild) 数据集
faces = fetch_lfw_people(min_faces_per_person=60)
# 数据集包含每个人至少 60 张脸部图像

# 查看数据集的形状
# print(faces.images.shape) # (1348, 62, 47)
# 62 是每个图像的行数，47 是每个图像的列数，因此每个图像有 62*47 = 2914 个像素
# 数据集中共有 1348 张图像

# print(faces.data.shape) # (1348, 2914)
# 数据集的形状为 (1348, 2914)，表示有 1348 个样本，每个样本有 2914 个特征（像素值）

# 创建图像网格，显示前 24 张图像，每个子图不显示坐标轴刻度
fig, axes = plt.subplots(3, 8, figsize=(8, 4), subplot_kw={"xticks": [], "yticks": []})
for i, ax in enumerate(axes.flat): # axes 是一个 3x8 的数组，包含了 24 个 AxesSubplot 对象。axes.flat 将其展开为一个一维数组
    ax.imshow(faces.images[i, :, :], cmap='gray') # 获取第 i 张图像，这是一张 2D 数组（灰度图像）

plt.show()

# 原本有2914维，现在试着降到200维

# 将数据矩阵 X 定义为 faces 数据集的 data 属性
X = faces.data # # (1348, 2914)
pca_face = PCA(n_components=200).fit(X)
X_face=pca_face.fit_transform(X) # (1384,200)
V = pca_face.components_
print('奇异值分解中的V的形状（即特征矩阵）:', V.shape, '\n200个主成分保留了原数据多少信息（累计方差贡献率）：',
      pca_face.explained_variance_ratio_.sum())

fig, axes = plt.subplots(3, 8, figsize=(8, 4), subplot_kw={"xticks": [], "yticks": []})
for i, ax in enumerate(axes.flat): 
    ax.imshow(V[i, :].reshape(62,47), cmap='gray') # 可见第一主成分（第一张图）获取的是人的轮廓信息，其他主成分不好解释……

# 使用inverse_transform可以逆转回去（但这种逆转不是完全可逆的，毕竟降维时只使用了200个主成分）
X_inv=pca_face.inverse_transform(X_face)

# 可视化逆转回的数据
fig, ax = plt.subplots(2, 8, figsize=(8, 2.5), subplot_kw={"xticks": [], "yticks": []})
# 现在的ax中是2行8列，作为对比：第一行是原数据，第二行是inverse_transform后返回的数据
# 一次循环画同一列上的两张图，而不是把ax拉平
for i in range(8):
    ax[0, i].imshow(faces.images[i, :, :], cmap="binary_r")
    ax[1, i].imshow(X_inv[i].reshape(62, 47), cmap="binary_r")
plt.show()

奇异值分解中的V的形状（即特征矩阵）: (200, 2914) 
200个主成分保留了原数据多少信息（累计方差贡献率）： 0.9420423

案例3：用PCA做噪音过滤（手写数字识别）

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
import numpy as np

digits = load_digits()
# digits.data.shape  # (1797,64)


# 定义一个画图的函数
def plot_digits(data):
    fig, axes = plt.subplots(4, 10, figsize=(10, 4), subplot_kw={"xticks": [], "yticks": []})
    for i, ax in enumerate(axes.flat):
        ax.imshow(data[i].reshape(8, 8), cmap="binary")
    plt.show()

plot_digits(digits.data)

# 将均值为0、标准差为2的正态分布随机数添加到 digits.data 中，从而人为创建一个带有噪声的新数据集 noisy
np.random.RandomState(42)
noisy = np.random.normal(digits.data, 2)
plot_digits(noisy)

# 利用PCA的inverse_transform函数去噪
pca_digits = PCA(n_components=0.7, svd_solver="full").fit(noisy)
noisy_ = pca_digits.fit_transform(noisy)
no_noisy_data = pca_digits.inverse_transform(noisy_)
plot_digits(no_noisy_data)

案例4：手写数字识别

对于某些算法，如果数据维度过大那么运行时间将会非常久，例如KNN。故利用PCA把原数据进行降维后再运用其他算法进行后续工作是常见的做法。本处案例进行手写数字识别。

from sklearn.decomposition import PCA
from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.neighbors import KNeighborsClassifier as KNN
from sklearn.model_selection import cross_val_score
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

plt.rcParams["font.sans-serif"]=["SimHei"] #设置字体
plt.rcParams["axes.unicode_minus"]=False #该语句解决图像中的负号的乱码问题

data = pd.read_csv("D:/PycharmProjects/ML20230605/机器学习/3、数据预处理与特征工程/digit recognizor.csv")
X = data.iloc[:, 1:] # (42000,784)
y = data.iloc[:, 0]

# 选择主成分个数，先设置200观察一下图像
pca_line = PCA(n_components=200).fit(X)
plt.figure(figsize=[20, 5])
plt.plot(np.cumsum(pca_line.explained_variance_ratio_))
plt.xlabel("主成分个数")
plt.ylabel("累计方差贡献率")
plt.show()

# 看图，缩小主成分范围。使用随机森林和KNN建模
# score_r = []
score_knn = []
for i in range(25, 30):
    X_dr = PCA(n_components = i).fit_transform(X)
#     once_r = cross_val_score(RFC(n_estimators=10, random_state=0), X_dr, y, cv=5).mean()
    once_knn = cross_val_score(KNN(), X_dr, y, cv=5).mean() 
#     score_r.append(once_r)
    score_knn.append(once_knn)
plt.figure(figsize=[20, 5])
# plt.plot(range(25, 30), score_r)
plt.plot(range(25, 30), score_knn)
plt.show() # 于是选择29作为主成分个数

X_dr=PCA(n_components=29).fit_transform(X)
print(cross_val_score(KNN(), X_dr, y, cv=5).mean()) # KNN未进行任何调参时

# 对KNN模型进行调参
score = []
for i in range(3,10):
    once = cross_val_score(KNN(i), X_dr, y, cv=5).mean()
    score.append(once)
plt.figure(figsize=[20, 5])
plt.plot(range(3,10), score)
plt.show() # 发现K=5时准确率最高，其实这就是KNN的默认参数

0.9720952380952381

习题

16.1

import numpy as np
from sklearn.decomposition import PCA

# 定义数据矩阵
X = np.array([[2, 3, 3, 4, 5, 7], [2, 4, 5, 5, 6, 8]])
X = X.astype("float64")

# 规范化变量
avg = np.average(X, axis=1)
var = np.var(X, axis=1)
for i in range(X.shape[0]):
    X[i] = (X[i, :] - avg[i]) / np.sqrt(var[i])

# 使用PCA类
pca = PCA(n_components=2)

# 拟合数据并转换
X_pca = pca.fit_transform(X.T)

# 打印结果
print("特征向量 V：\n", pca.components_.T) # pca.components_ 属性返回的是形状为 (n_components, n_features) 的矩阵
print("样本主成分矩阵 Y：\n", X_pca.T)

# 使用自定义的pca_svd函数
def pca_svd(X, k):
    """
    样本矩阵的奇异值分解的主成分分析算法
    :param X: 样本矩阵X
    :param k: 主成分个数k
    :return: 特征向量V，样本主成分矩阵Y
    """
    n_samples = X.shape[1]  # 这里的X矩阵的一列表示一个观测

    # 构造新的n×m矩阵
    T = X.T / np.sqrt(n_samples - 1)

    # 对矩阵T进行截断奇异值分解
    U, S, Vt = np.linalg.svd(T)
    V = Vt.T[:, :k]

    # 求k×n的样本主成分矩阵
    return V, np.dot(V.T, X)

# 设置精度为3
np.set_printoptions(precision=3, suppress=True)
V, vnk = pca_svd(X, 2)

print("正交矩阵V：")
print(V)
print("样本主成分矩阵Y：")
print(vnk)

特征向量 V：
 [[ 0.707 -0.707]
 [ 0.707  0.707]]
样本主成分矩阵 Y：
 [[-2.028 -0.82  -0.433 -0.     0.82   2.461]
 [-0.296  0.046  0.433  0.    -0.046 -0.137]]
正交矩阵V：
[[ 0.707  0.707]
 [ 0.707 -0.707]]
样本主成分矩阵Y：
[[-2.028 -0.82  -0.433  0.     0.82   2.461]
 [ 0.296 -0.046 -0.433  0.     0.046  0.137]]

16.2、16.3

见该链接

使用书籍：李航《机器学习方法》
习题解答：https://datawhalechina.github.io/statistical-learning-method-solutions-manual/#/