PCA的数学原理

1 预备知识

1.1 奇异值分解 SVD

任何 $m\times n$ 实矩阵 A, $A\in R^{m\times n}$ 的奇异值分解一定存在，表示为三个矩阵的乘积：
$$A=U\Sigma V^{\mathsf{T}}$$
其中 U 是 m 阶正交矩阵，V 是 n 阶正交矩阵; Σ 是 m×n 矩形对角矩阵，其对角线元素称为奇异值，值非负，且按降序排列，这种分解方式称为矩阵的完全奇异值分解

在实际应用中，SVD 被用来对矩阵做近似表示，也即 $U\Sigma V^{\mathsf{T}}$ 得到的结果并不需要与 A 完全相同，只要值尽量相同即可。这时就可以取 Σ 中最大的 k 个奇异值，得到矩阵的截断奇异值分解，这时有
$$A\approx U_k{\Sigma }_k{V}_k^{\mathsf{T}}$$
其中 $U_k$是 m×k 矩阵， $V_k$是 n×k 矩阵，${\Sigma }_k$ 是 k 阶对角矩阵； 矩阵 $U_k$由完全奇异值分解中 U 的前 k 列得到、矩阵$V_k$由完全奇异值分解中 V 的前 k 列得到、矩阵 ${\Sigma }_k$ 由 $\Sigma$ 的前 k 个对角元素得到

如下列矩阵 A 完全奇异值分解为
$$A={\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 4\\ 0& 3& 0& 0\\ 0& 0& 0& 0\\ 2& 0& 0& 0\end{array}\right)}_{5\times 4}={\left(\begin{array}{ccccc}0& 0& \sqrt{0.2}& 0& \sqrt{0.8}\\ 1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& \sqrt{0.8}& 0& -\sqrt{0.2}\end{array}\right)}_{5\times 5}{\left(\begin{array}{cccc}4& 0& 0& 0\\ 0& 3& 0& 0\\ 0& 0& \sqrt{5}& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)}_{5\times 4}{\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right)}_{4\times 4}$$
矩阵 A 截断奇异值分解为(取 2 个奇异值)

$$A={\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 4\\ 0& 3& 0& 0\\ 0& 0& 0& 0\\ 2& 0& 0& 0\end{array}\right)}_{5\times 4}\approx {\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 4\\ 0& 3& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)}_{5\times 4}={\left(\begin{array}{cc}1& 0\\ 0& 0\\ 0& 3\\ 0& 0\\ 2& 0\end{array}\right)}_{5\times 2}{\left(\begin{array}{cc}4& 0\\ 0& 3\end{array}\right)}_{2\times 2}{\left(\begin{array}{cc}0& 0\\ 0& 1\\ 1& 0\\ 0& 0\end{array}\right)}_{2\times 4}^{\mathsf{T}}$$
可以看到可以两者结果不完全相同，截断奇异值分解对矩阵做做了近似，精度有所损失。对于高维度矩阵截断奇异值分解可以大大提高压缩空间,如原始矩阵 A 大小为 20000=100 * 200, 完全奇异值分解后大小为 70000=200 * 200+200 * 100+100 * 100,截断奇异值分解后存储大小为(取前10个奇异值) 2100=100 * 10+10* 10+10 * 200 ，压缩到原来大小的 1/10

1.2 协方差矩阵

理解协方差矩阵的关键就在于牢记它的计算是不同维度之间的协方差，而不是不同样本之间。拿到一个样本矩阵，最先要明确的就是一行是一个样本还是一个维度

有矩阵 $X_{m\times n}$ ,m 表示数据维度， n 表示样本个数，协方差矩阵即为
$$\Sigma =\frac{1}{n}X{X}^T$$

2 PCA 原理

2.1 原理

下图为二维数据，横轴为 $x_1$,纵轴为 $x_2$,可以看出 $x_1$ 与 $x_2$ 是线性相关的。具体的，当知道其中一个变量如 $x_1$的取值时,对另一个变量$x_2$的预测不是完全随机的，反之亦然

我们希望找到一组轴，将数据投影到新轴上后各维数据均线性无关，也即知道一个轴的值后对其它维数据的猜测是随机的

对于给定的 n 条 m 维数据，{$v_1,v_2,...,v_n$}, 其中所有向量$v_i$均列向量，中心化后表示为{$x_1,x_2,...,x_n$}={$v_1-\text{u},v_2-\text{u},...,v_n-\text{u}$}，其中 $\text{u}=\frac{1}{n}{\sum }_{i=1}^n{v}_i$。

向量内积在几何上表示为第一个向量投影到第二个向量的长度(参考 PCA的数学原理)，因此向量 $x_i$ 在 $\text{w}$（单位方向向量）上的投影坐标可以表示为内积形式 $(x_i,\text{w})=x_i^T\text{w}$,所以目标是找到一个投影方向，使得 {$x_1,x_2,...,x_n$} 在 $w$ 的投影方差尽可能大。易知，投影之后均值为 $\text{0}$(因为 $u^{{}^{\prime }}=\frac{1}{n}{\sum }_{i=1}^n{x}_i^{T}w=(\frac{1}{n}{\sum }_{i=1}^n{x}_i^T)w=0$)

方差计算公式，方差为值将去均值平方的期望，最后可以转化为平方的期望减去均值的平方

因为均值为0，因此，投影后的方差可以表示为

上式中间部分 $(\frac{1}{n}{\sum }_{i=1}^n{w}^T{x}_i{x}_i^{T}w)$ 其实就是协方差矩阵，协方差矩阵可以写为 Σ。由于要求 $w$ 为单位向量，即$w^{T}w=1$。因此,最后即求一个最大化问题
$$max(w^T\Sigma w)$$
$$s.t.(w^{T}w=1)$$

引入拉格朗日乘子， 并对 $w$ 求导令其等于 $0$
$$L(w,\lambda )=w^{\mathsf{T}}\Sigma w-\lambda (w^{\mathsf{T}}w-1)$$
对 w 求导有：
$$w^{\mathsf{T}}(\Sigma +{\Sigma }^{\mathsf{T}})-\lambda \ast 2w^{\mathsf{T}}=2w^{\mathsf{T}}\Sigma -2\lambda w^{\mathsf{T}}$$
令其为 0,于是有 $w^{\mathsf{T}}\Sigma =\lambda w^{\mathsf{T}}=>w^{\mathsf{T}}\Sigma w=\lambda$
(注:Σ为对称矩阵，有 ${\Sigma }^{\mathsf{T}}=\Sigma$)

(注：
标量对向量求导公式

)

由此可以看到，找到的最大的方差也就是协方差矩阵最大的特征值，最佳投影方向就是最大特征值所对应的特征向量。次佳投影方向与最佳投影方向正交，是第二大特征值对应的特征向量，以此类推

2.2 算法

计算矩阵 $X_{m\times n}$ 降维的方法可以按照下面步骤计算

(1) 样本组成m行n列数据，m表示维度，n表示样本数，即样本按列排布
(2) 对样本数据进行中心化处理
(3) 求样本协方差矩阵 $C=\frac{1}{n}X{X}^T$
(4) 求协方差矩阵的特征值和特征向量，将特征值从大到小排列，特征向量按行排列
(5) 取特征值前 d 大对应的特征向量 $w_1,w_2,...,w_d$，即特征向量的前 d 行，组成矩阵 $P_{d\times m}$
(6) Y=PX 即为降维到 d 维后的数据

定义降维之后的信息占比为：

2.3 例子

import numpy as np
# stp1 数据准备 m*n  n条数据m维
x =np.array([[-1,-1,0,2,0],
             [-2,0,0,1,1]])
# 每一维均值，此例中为2维，数组最后变为 2*1
x_mean = x.mean(axis=1)[:, np.newaxis]
# step2 每一维减去均值
x_new = x-x_mean
# step3 求协方差矩阵
x_conv = np.cov(x_new)
# step4 求出协方差矩阵的特征值及对应的特征向量
eigen_vals, eigen_vecs = np.linalg.eig(x_conv)#特征值是eigen_vals，特征向量是eigen_vecs， 
                                                #特征向量是按列放的，即一列代表一个特征向量

eigen_vals = eigen_vals[:, np.newaxis]
print(eigen_vals)
# [[2.5]
# [0.5]]

# 变行向量为特征向量
eigen_vecs = eigen_vecs.T
print(eigen_vecs)
#[[ 0.70710678  0.70710678]
# [-0.70710678  0.70710678]]

val_vec = np.hstack((eigen_vals, eigen_vecs))
# 特征值由大到小排序，[::-1] 表示倒序
val_vec = val_vec[np.argsort(val_vec[:,0])[::-1]]
print(val_vec)
#[[ 2.5         0.70710678  0.70710678]
#[ 0.5        -0.70710678  0.70710678]]

# 取前 d 个特征向量
d = 1
p = val_vec[:d, 1:]
print(p)
# [[0.70710678 0.70710678]]

# 对数据进行转化, 5条2维数据转为 5条1维数据
x_new = np.dot(p, x)
print(x_new)
# [[-2.12132034 -0.70710678  0.          2.12132034  0.70710678]]

参考

1 PCA的数学原理
2 统计学习方法 第 2 版 李航
3 an introduction to Principal Component Analysis(PCA)
4 期望、方差、协方差及相关系数的基本运算
5【五万字总结PCA】【李航|PRML】统计学习方法–16.主成分分析（详细推导）
6 百面机器学习