第16章 主成分分析

习题16.1

  对以下样本数据进行主成分分析：

$$X=\left[\begin{array}{llllll}2& 3& 3& 4& 5& 7\\ 2& 4& 5& 5& 6& 8\end{array}\right]$$

解答：

解答思路：

1. 给出主成分分析算法
2. 自编程基于numpy库实现主成分分析算法

解答步骤：

第1步：主成分分析算法

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

1. 相关矩阵的特征值分解算法

  根据书中第16.2.2节的相关矩阵的特征值分解算法具体步骤：

（1）对观测数据按照如下公式进行规范化处理，得到规范化数据矩阵，仍以$X$表示

$$x_{ij}^{\ast }=\frac{x_{ij}-{\overline{x}}_i}{\sqrt{s_{ii}}},i=1,2,\cdots ,m;j=1,2,\cdots ,n$$

其中

$$\begin{array}{l}{\displaystyle {\overline{x}}_i=\frac{1}{n}\sum _{j=1}^n{x}_{ij},i=1,2,\cdots ,m}\\ {\displaystyle s_{ii}=\frac{1}{n-1}\sum _{j=1}^n(x_{ij}-{\overline{x}}_i{)}^2,i=1,2,\cdots ,m}\end{array}$$

（2）依据规范化数据矩阵，计算样本相关矩阵$R$

$$R=[r_{ij}{]}_{m\times m}=\frac{1}{n-1}X{X}^T$$

其中

$$r_{ij}=\frac{1}{n-1}\sum _{l=1}^n{x}_{il}{x}_{lj},i,j=1,2,\cdots ,m$$

（3）求样本相关矩阵$R$的$k$个特征值和对应的$k$个单位特征向量。 求解$R$的特征方程

$$|R-\lambda I|=0$$

得$R$的$m$个特征值

$${\lambda }_1⩾{\lambda }_2⩾\cdots ⩾{\lambda }_m$$

求方差贡献度${\displaystyle \sum _{i=1}^k{\eta }_i}$达到预定值的主成分个数$k$。
求前$k$个特征值对应的单位特征向量

$$a_i=(a_{1i},a_{2i},\cdots ,a_{mi}{)}^T,i=1,2,\cdots ,k$$

（4）求$k$个样本主成分
以$k$个单位特征向量为系数进行线性变换，求出$k$个样本主成分

$$y_i=a_i^{T}x,i=1,2,\cdots ,k$$

（5）计算$k$个主成分$y_i$与原变量$x_i$的相关系数$\rho (x_i,y_i)$，以及$k$个主成分对原变量$x_i$的贡献率$v_i$。

（6）计算$n$个样本的$k$个主成分值
将规范化样本数据带入$k$个主成分式，得到$n$个样本的主成分值。第$j$个样本$x_j=(x_{1j},x_{2j},\cdots ,x_{mj}{)}^T$的第$i$个主成分值是

$$\begin{gathered} y_{ij}=(a_{1i},a_{2i},\cdots ,a_{mi})(x_{1j},x_{2j},\cdots ,x_{mj}{)}^T=\sum _{l=1}^m{a}_{li}{x}_{lj} \\ i=1,2,\cdots ,m,j=1,2,\cdots ,n \end{gathered}$$

  根据书中第16章的定义16.2给出了方差贡献率：

定义16.2 第$k$主成分$y_k$的方差贡献率定义为$y_k$的方差与所有方差之和的比，记作${\eta }_k$

$$\begin{array}{}\text{(16.30)}& {\eta }_k=\frac{{\lambda }_k}{{\displaystyle \sum _{i=1}^m{\lambda }_i}}\end{array}$$

$k$个主成分$y_1,y_2,\cdots ,y_k$的累计方差贡献率定义为$k$个方差之和与所有方差之和的比

$$\begin{array}{}\text{(16.31)}& \sum _{i=1}^k{\eta }_i=\frac{{\displaystyle \sum _{i=1}^k{\lambda }_i}}{{\displaystyle \sum _{i=1}^m{\lambda }_i}}\end{array}$$

2. 样本矩阵的奇异值分解算法

  根据书中第16章的算法16.1主成分分析算法

算法16.1（主成分分析算法）
输入：$m\times n$样本矩阵$X$，其每一行元素的均值为零；
输出：$k\times n$样本主成分矩阵$Y$。
参数：主成分个数$k$
（1）构造新的$n\times m$矩阵

$$X^{\prime }=\frac{1}{\sqrt{n-1}}{X}^T$$

$X^{\prime }$每一列的均值为零。
（2）对矩阵$X^{\prime }$进行截断奇异值分解，得到

$$X^{\prime }=U\mathrm{\Sigma }{V}^T$$

有$k$个奇异值、奇异向量。矩阵$V$的前$k$列构成$k$个样本主成分。
（3）求$k\times n$样本主成分矩阵

$$Y=V^{T}X$$

第2步：自编程基于numpy库实现主成分分析算法

import numpy as np


def pca_svd(X, k):
    """
    样本矩阵的奇异值分解的主成分分析算法
    :param X: 样本矩阵X
    :param k: 主成分个数k
    :return: 特征向量V，样本主成分矩阵Y
    """
    n_samples = X.shape[1]

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

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

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

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])

# 设置精度为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]]

习题16.2

  证明样本协方差矩阵$S$是总体协方差矩阵$\mathrm{\Sigma }$的无偏估计。

解答：

解答思路：

1. 给出总体协方差矩阵$\mathrm{\Sigma }$
2. 给出样本协方差矩阵$S$
3. 给出无偏估计的定义
4. 证明样本协方差矩阵$S$是总体协方差矩阵$\mathrm{\Sigma }$的无偏估计

解答步骤：

第1步：总体协方差矩阵$\mathrm{\Sigma }$

  根据书中第16章的定理16.1：

定理16.1 设$x$是$m$维随机变量，$\mathrm{\Sigma }$是$x$的协方差矩阵，$\mathrm{\Sigma }$的特征值分别是${\lambda }_1⩾{\lambda }_2⩾\cdots ⩾{\lambda }_m⩾0$，特征值对应的单位特征向量分别是${\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_m$，则$x$的第$k$个主成分是

$$y_k={\alpha }_k^{T}x={\alpha }_{1k}{x}_1+{\alpha }_{2k}{x}_2+\cdots +{\alpha }_{mk}{x}_m,k=1,2,\cdots ,m$$

$x$的第$k$主成分的方差是

$$\text{var}(y_k)={\alpha }_k^T\mathrm{\Sigma }{\alpha }_k={\lambda }_k,k=1,2,\cdots ,m$$

即协方差矩阵$\mathrm{\Sigma }$的第$k$个特征值。

  根据书中第16章的本章概要的总体主成分的性质：

主成分$y$的协方差矩阵是对角矩阵

$$\text{cov}(y)=\mathrm{\Lambda }=\text{diag}({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_m)$$

第2步：样本协方差矩阵$S$

  根据书中第16.2.1的样本协方差矩阵定义：

  给定样本矩阵$X$，可以估计样本均值，以及样本协方差。样本均值向量$\overline{x}$为

$$\begin{array}{}\text{(16.40)}& \overline{x}=\frac{1}{n}\sum _{j=1}^n{x}_j\end{array}$$

样本协方差矩阵$S$为

$$\begin{array}{}\text{(16.41)}& \begin{array}{l}S=[s_{ij}{]}_{m\times m}\\ {\displaystyle s_{ij}=\frac{1}{n-1}\sum _{k=1}^n(x_{ik}-{\overline{x}}_i)(x_{jk}-{\overline{x}}_j),i,j=1,2,\cdots ,m}\end{array}\end{array}$$

第3步：无偏估计的定义

  根据《实用多元统计分析》（方开泰，1989年9月第一版）中第88~89页的无偏性：

如果参数$\theta$的某个估计$\hat{\theta }$满足$E(\hat{\theta })=\theta$，则称$\hat{\theta }$是无偏的。······引理3.4还指出：

$$\hat{\mathrm{\Sigma }}=\frac{1}{n}{A}^d=\frac{1}{n}\sum _{\alpha =1}^{n-1}{z}_{\alpha }{z}_{\alpha }^{\prime }$$

$z_1,\cdots ,z_{n-1}$独立同分布于$N_p(0,\mathrm{\Sigma })$，于是

$$E(\hat{\mathrm{\Sigma }})=\frac{1}{n}\sum _{\alpha =1}^{n-1}{z}_{\alpha }{z}_{\alpha }^{\prime }=\frac{n-1}{n}\mathrm{\Sigma }$$

$\hat{\mathrm{\Sigma }}$不是$\mathrm{\Sigma }$的无偏估计，但可修正一下使之为无偏估计。令

$$S=\frac{1}{n-1}A$$

则$S$是$\mathrm{\Sigma }$的无偏估计。在实用中普遍采用$S$来估计$\mathrm{\Sigma }$。

  无偏估计是用样本统计量来估计总体参数的一种无偏推断。估计量的期望等于被估计参数的真实值，则称此估计量为被估计参数的无偏估计，即具有无偏性。表示在多次重复下，它们的平均数接近所估计的参数真值。

第4步：证明样本协方差矩阵$S$是总体协方差矩阵$\mathrm{\Sigma }$的无偏估计

  由样本协方差矩阵$S$定义：

$$S=\frac{1}{n-1}\sum _{i=1}^n(x_i-\overline{x})((x_i-\overline{x}){)}^{\prime }$$

  令：

$$A=\sum _{i=1}^n(x_i-\overline{x})((x_i-\overline{x}){)}^{\prime }$$

  根据无偏估计的定义，证明样本协方差矩阵$S$是总体协方差矩阵$\mathrm{\Sigma }$的无偏估计，即证明样本协方差矩阵$S$的期望等于总体协方差矩阵$\mathrm{\Sigma }$：

$$E(S)=E(\frac{1}{n-1}A)=\frac{1}{n-1}E(A)=\mathrm{\Sigma }$$

即证明：

$$E(A)=(n-1)\mathrm{\Sigma }$$

  参考《应用多元统计分析 第五版》（王学民，2017年8月第5版）书中第51页的无偏性一节的推导

依(2.2.14)式

$$V(\overline{x})=\frac{1}{n^2}\sum _{i=1}^{n}V(x_i)=\frac{1}{n^2}\sum _{i=1}^n\mathrm{\Sigma }=\frac{1}{n}\mathrm{\Sigma }$$

于是

$$\begin{array}{rl}E(\hat{\mathrm{\Sigma }})& =\frac{1}{n}E[\sum _{i=1}^n(x_i-\overline{x})(x_i-\overline{x}{)}^{\prime }]\\ & =\frac{1}{n}E\{\sum _{i=1}^n[(x_i-\mu )-(\overline{x}-\mu )][(x_i-\mu )-(\overline{x}-\mu ){]}^{\prime }\}\\ & =\frac{1}{n}E[\sum _{i=1}^n(x_i-\mu )(x_i-\mu {)}^{\prime }-n(\overline{x}-u)(\overline{x}-u{)}^{\prime }]\\ & =\frac{1}{n}[\sum _{i=1}^{n}V(x_i)-nV(\overline{x})]\\ & =\frac{1}{n}(n\mathrm{\Sigma }-n\cdot \frac{1}{n}\mathrm{\Sigma })\\ & =\frac{n-1}{n}\mathrm{\Sigma }\end{array}$$

$$\begin{array}{rl}E(A)& =E[\sum _{i=1}^n(x_i-\overline{x})(x_i-\overline{x}{)}^{\prime }]\\ & =nE(\hat{\mathrm{\Sigma }})\\ & =n\cdot \frac{n-1}{n}\mathrm{\Sigma }\\ & =(n-1)\mathrm{\Sigma }\end{array}$$

  因此，样本协方差矩阵$S$是总体协方差矩阵$\mathrm{\Sigma }$的无偏估计。

习题16.3

  设$X$为数据规范化样本矩阵，则主成分等价于求解以下最优化问题：

$$\begin{array}{cl}\underset{L}{min}& \parallel X-L{\parallel }_F\\ \text{s.t.}& \text{rank}(L)⩽k\end{array}$$

其中，$F$是弗罗贝尼乌斯范数，$k$是主成分个数。试问为什么？

解答：

解答思路：

1. 给出弗罗贝尼乌斯范数的定义
2. 给出矩阵的最优近似
3. 证明主成分是该最优化问题的最优解，该最优化问题的最优解是$X$的主成分

解答步骤：

第1步：弗罗贝尼乌斯范数

  根据书中第15章的定义15.4的弗罗贝尼乌斯范数：

  定义15.4（弗罗贝尼乌斯范数） 设矩阵$A\in R^{m\times m}$，$A=[a_{ij}{]}_{m\times m}$，定义矩阵$A$的弗罗贝尼乌斯范数为

$$\parallel A{\parallel }_F={(\sum _{i=1}^m\sum _{j=1}^n(a_{ij}{)}^2)}^{\frac{1}{2}}$$

第2步：矩阵的最优近似

  根据书中第15.3.2节的定理15.3：

  定理15.3 设矩阵$A\in R^{m\times n}$，矩阵的秩$\text{rank}(A)=r$，有奇异值分解$A=U\mathrm{\Sigma }{V}^T$，并设$\mathcal{M}$为$R^{m\times n}$中所有秩不超过$k$的矩阵的集合，$0<k<r$，若秩为$k$的矩阵$X\in \mathcal{M}$满足

$$\Vert A-X{\Vert }_F=\underset{S\in \mathcal{M}}{min}\Vert A-S{\Vert }_F$$

则

$$\Vert A-X{\Vert }_F=({\sigma }_{k+1}^2+{\sigma }_{k+2}^2+\cdots +{\sigma }_n^2{)}^{\frac{1}{2}}$$

特别地，若$A^{\prime }=U{\mathrm{\Sigma }}^{\prime }{V}^T$，其中

$${\mathrm{\Sigma }}^{\prime }=\left[\begin{array}{cccccc}{\sigma }_1& & & & & \\ & \ddots & & & 0& \\ & & {\sigma }_k& & & \\ & & & 0& & \\ & 0& & & \ddots & \\ & & & & & 0\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\Sigma }}_k& 0\\ 0& 0\end{array}\right]$$

则

$$\Vert A-A^{\prime }\Vert =({\sigma }_{k+1}^2+{\sigma }_{k+2}^2+\cdots +{\sigma }_n^2{)}^{\frac{1}{2}}=\underset{S\in \mathcal{M}}{min}\Vert A-S{\Vert }_F$$

第3步：证明命题

  根据样本矩阵的奇异值分解算法（主成分分析算法），针对$m\times n$样本矩阵$X$，需要表示为

$$X^{\prime }=\frac{1}{\sqrt{n-1}}{X}^T$$

对优化问题的目标函数进行如下变换：

$$\begin{array}{cl}\underset{L}{min}& {\displaystyle {\Vert \frac{1}{\sqrt{n-1}}{X}^T-\frac{1}{\sqrt{n-1}}{L}^T\Vert }_F}\\ \text{s.t.}& \text{rank}(L)⩽k\end{array}$$

  令

$$X^{\prime }=\frac{1}{\sqrt{n-1}}{X}^T,L^{\prime }=\frac{1}{\sqrt{n-1}}{L}^T$$

  则原最优化问题等价于

$$\begin{array}{cl}\underset{L}{min}& \parallel X^{\prime }-L^{\prime }{\parallel }_F\\ \text{s.t.}& \text{rank}(L)⩽k\\ & {\displaystyle X^{\prime }=\frac{1}{\sqrt{n-1}}{X}^T}\\ & {\displaystyle L^{\prime }=\frac{1}{\sqrt{n-1}}{L}^T}\end{array}$$

  对$X^{\prime }$进行奇异值分解，根据已知主成分个数为$k$，$X^{\prime }$的主成分个数也为$k$个，则

$$X^{\prime }=U^{\prime }{\mathrm{\Sigma }}^{\prime }{V}^{\prime T}$$

  结合第2步，可知变换后的最优化问题和定理15.3类似，形如：

$$\begin{gathered} A\sim X^{\prime } \\ S\sim L^{\prime } \end{gathered}$$

  可找到 ${\mathrm{\Sigma }}^{\prime }=\text{diag}({\lambda }_1^{\prime },{\lambda }_2^{\prime },\cdots ,{\lambda }_k^{\prime },0,0,\cdots ,0)$ 使得变换后的目标函数取最小值，即存在矩阵满足主成分条件。

参考文献

【1】方开泰. 实用多元统计分析[M]. 1989-09:88-89
【2】王学民. 应用多元统计分析 第五版[M]. 2017-08:51