【原创】PCA算法详细推导

bbbb北邮模式识别实验室

• Author：李英民 | Henry
• E-mail: li_yingmin@outlookdotcom
• $Beijing~~U~of~Posts~and~Telecom$$(BUPT~)$
• Home: https://liyingmin.wixsite.com/henry

• 快速了解我: About Me

• 转载请保留上述引用内容，谢谢配合！

---

一、线性代数基本知识

1、特征值和特征向量

定义是这样的，设A为n阶方阵，如果存在数$\lambda$,使得$A\alpha=\lambda\alpha$，称$\lambda$为特征值，$\alpha$为特征向量。

2、特征值和特征向量的意义

某方阵$A$经过其特征向量构成的矩阵的线性变化，可以成为对角阵。

即，如果有个方阵$A$，其所有特征值组成的特征向量是下面这样的：

$$U=(\alpha_1，\alpha_2，\alpha_3 \cdots，\alpha_n)^T$$
那么方阵$A$就有可能（不能保证）可以经过线性变化成一个对角阵：
$$UAU^{-1}= \left ( \begin{array}{c} \lambda_1& 0 & 0 & 0 & \cdots\\ 0 & \lambda_2 & \cdots & 0 & 0 \\ 0 & 0 & \lambda_3 & \cdots &0 \\ 0 & 0 & 0 &\lambda_4 & \cdots \end{array} \right)$$

3、线性变化的意义

线性变化的意义就是坐标轴的旋转，当然，这只在一位到三维空间比较容易想象，实际上对于高维，我们一般称之为基的选择变化

如果选择的一组基是标准正交基，那么就是标准正交变换。

标准正交基就是每组基之间两两线性无关，而且相互正交。

上边的$\alpha_1,\alpha_2$就是一组基，但是我们还不能保证其一定线性无关，更不能保证其一定正交。

4、对称矩阵的特殊性

• 实对称矩阵的$\lambda$都是实数

• 实对称矩阵对应不同特征值的特征向量是正交的。

• 实对称矩阵对应相同的特征值（重根）的特征向量一定是线性无关的，可以用方法（施密特正交化方法）使其正交

• 实对称矩阵一定可以被某个正交矩阵$U$对角化：

• $$UAU^{-1}= \left ( \begin{array}{c} \lambda_1& 0 & 0 & 0 & \cdots\\ 0 & \lambda_2 & \cdots & 0 & 0 \\ 0 & 0 & \lambda_3 & \cdots &0 \\ 0 & 0 & 0 &\lambda_4 & \cdots \end{array} \right)$$

  ​

• 正交矩阵有$UU^T=E$(单位阵)，$U^{-1}=U^T$

5、一个性质的补充

对于二维平面的一堆没有明显规律的点集$X$:

$$X= \left ( \begin{array}{c} x_1& x_2 & \cdots & x_n \\ y_1 & y_2 & \cdots & y_n \\ \end{array} \right)$$
我们需要利用一组标准正交基，将$X$表示的点集，用另一组坐标表示成$Y$:
$$Y= \left ( \begin{array}{c} a_1& a_2 & \cdots & a_n \\ b_1 & b_2 & \cdots & b_n \\ \end{array} \right)$$
根据线性代数的知识，旋转坐标轴意味着矩阵的线性变换，所以根据线性变换的原理，我们可以得到：
$$x=Uy$$
其中，
$$x=(x_1,y_1),y=(a_1,b_1)$$

$$U=\left ( \begin{array}{c} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \\ \end{array} \right)$$

表示点集中的任意一点.

如果已知：

$$\left( \begin{array}{c} a_1 \\ b_1 \end{array} \right) = \left ( \begin{array}{c} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \\ \end{array} \right) \left ( \begin{array}{c} x_1 \\ y_1 \\ \end{array} \right )$$

那么一定有：

$$\left ( \begin{array}{c} a_1-\bar{a_1} \\ b_1-\bar{b_1} \\ \end{array} \right) = \left ( \begin{array}{c} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \\ \end{array} \right) \left ( \begin{array}{c} x_1-\bar{x_1} \\ y_1-\bar{y_1} \\ \end{array} \right)$$

其中，$\bar{a_1}$和$\bar{b_1}$分别是$ab$坐标的平均值，$\bar{x_1}\,and \,\bar{x_2}$分别是$xy$坐标的平均值.

---

二、统计基本知识

1、二维观测数据如下

$$X= \left ( \begin{array}{c} x_1& x_2 & \cdots & x_n \\ y_1 & y_2 & \cdots & y_n \\ \end{array} \right)$$

其中 $(x_1,y_1)$ 等代表二维坐标系下的一个点.

2、均值和方差

2.1均值向量

$(\bar{x}\,，\bar{y})^{T}$表示均值向量

$$\bar{x}=\frac{1}{N}\sum_{i=1}^Nx_i$$

2.2方差

$s_{xx}$表示方差

$$s_{xx}=\frac{1}{N-1}\sum_{i=1}^N(x_i-\bar{x})^2$$

2.3协方差

$s_{xy}$表示协方差

$$s_{xy}=\frac{1}{N-1}\sum_{i=1}^N(x-\bar{x})(y-\bar{y})$$
可见$s_{xy}=s_{yx}$

3、协方差矩阵

$$S= \left ( \begin{array}{c} s_{xx} & s_{xy} \\ s_{yx} & s_{yy} \\ \end{array} \right) = \frac{1}{N-1} \left ( \begin{array}{c} x_1-\bar{x} & x_2-\bar{x} \\ y_1-\bar{y} & y_2-\bar{y} \end{array} \right) \left ( \begin{array}{c} x_1-\bar{x} & y_1-\bar{y} \\ x_2-\bar{x} & y_2-\bar{y} \end{array} \right).$$

所以协方差矩阵是一个对称矩阵

---

三、详细推导过程

三.一、坐标轴旋转

现在仍然以二维的数据作为理论分析的目标。对于二维平面的一堆没有明显规律的点集$X$:

$$X= \left ( \begin{array}{c} x_1& x_2 & \cdots & x_n \\ y_1 & y_2 & \cdots & y_n \\ \end{array} \right)$$
我们需要利用一组标准正交基，将$X$表示的点集，用另一组坐标表示成$Y$:
$$Y= \left ( \begin{array}{c} a_1& a_2 & \cdots & a_n \\ y_1 & b_2 & \cdots & b_n \\ \end{array} \right)$$
根据线性代数的知识，旋转坐标轴意味着矩阵的线性变换，所以根据线性变换的原理，我们可以得到：
$$y=Ux$$
其中，
$$x=(x,y),y=(a,b)$$
表示点集中的任意一点.

称坐标轴$a$为第一主分量，坐标轴$b$为第二主分量.

$U$是一组由标准正交基组成的正交矩阵：

$$U=(\alpha_1,\alpha_2)^T$$
然后可以根据新坐标下的数据方差的大小，适当省略方差较小的坐标轴表示的维度。比如：

新的坐标轴下，第二维数据$b_1$到$b_n$的方差很小，那么就意味着$b$这一维的数据对分类的贡献没有那么大，所以我们就把$b$这一维数据忽略掉。

下面举一个例子：

相当于投影.

三.二、为什么这样旋转

问为什么这样旋转，其实也就是要问，怎么求出来这个正交矩阵$U$。

1、求$U$的方法

对于这样一个点集$X$:

$$X= \left ( \begin{array}{c} x_1& x_2 & \cdots & x_n \\ y_1 & y_2 & \cdots & y_n \\ \end{array} \right)$$
我们可以求出它的协方差矩阵$S$:
$$S_x= \left ( \begin{array}{c} s_{xx} & s_{xy} \\ s_{yx} & s_{yy} \end{array} \right)$$
由于这个矩阵是实对称矩阵，所以它就满足以下性质：

• 实对称矩阵的特征值$\lambda$都是实数
• 实对称矩阵对应不同特征值的特征向量是正交的。
• 实对称矩阵对应相同的特征值（重根）的特征向量一定是线性无关的，可以用方法（施密特正交化方法）使其正交
• 实对称矩阵一定有$n$个正交的特征向量，这$n$个特征向量一定可以表示为一个正交矩阵

用任意一种求解特征值的算法，我们就一定可以求到矩阵$S$的两个特征值，$\lambda_1,\lambda_2$，让$\lambda_1>\lambda_2$，对应的特征向量分别为：$\alpha_1,\alpha_2$，同时还要将所有的特征向量标准正交化。

$$S_x= \left ( \begin{array}{c} s_{xx} & s_{xy} \\ s_{yx} & s_{yy} \end{array} \right)$$

$$U= (\alpha_1,\alpha_2)^T$$

而且这两个特征向量还满足实对阵矩阵的特性：

• 实对称矩阵一定有$n$个正交的特征向量，这$n$个特征向量一定可以表示为一个正交矩阵

所以矩阵$U$一定是个正交矩阵。至此，$U$就求出来了！

但还需要补充说明一点很重要的性质

对于实对称矩阵的正交变换，还有一个性质：

有$s_{xx}+s_{yy}=\lambda_1+\lambda_2$,但是不保证$s_{xx}=\lambda_1,s_{yy}=\lambda_2$，而且一般是没有这样关系的：因为变换坐标的目的就是让新坐标下每个维度的坐标对应的数据的方差依次减小，即原先数据的方差分别为$x$坐标下$s_{xx}$，$y$坐标下$s_{yy}$；而现在的新的坐标下，比如新的点集表示为$Y$：

$$y=Ux= \left ( \begin{array}{c} a\\ b \\ \end{array} \right)$$
在下文中会证明$Y$的协方差矩阵：
$$S_y= \left ( \begin{array}{c} s_{aa} & s_{ab} \\ s_{ba} & s_{bb} \end{array} \right)$$
满足$s_{aa}=\lambda_1$,$s_{bb}=\lambda_2$，所以，点集$X$经过坐标变换为点集$Y$，实际上是数据方差的重新分布。

2、为什么这样求$U$

其实上面那个补充说明已经点到了很重要的一点，就是我的点集原来在$xy$坐标下表示为$X$，即原先数据的方差分别为$x$坐标下$s_{xx}$，$y$坐标下$s_{yy}$，现在经过$Y=UX$的正交变换，坐标轴变成了$ab$.（注意，点集还是原来那些点，只不过是坐标轴发生了变化，表示的方式不同而已）而且a坐标下方差$s_{aa}$和$b$坐标下方差$s_{bb}$，满足$s_{aa}=\lambda_1$,$s_{bb}=\lambda_2$.

同时，我们还希望对于新的坐标系$ab$,坐标系的基底之间（即特征向量之间）满足相互正交:

$$\alpha_1 \alpha_2^T=\sum_{i=1}^2m_in_i=0$$
其中$\alpha_1=(m_1,m_2)^T$，$\alpha_2=(n_1,n_2)^T$.

下面证明我们求出来的$U$确实满足了这两个性质：

由于$\alpha_1=(m_1,m_2)^T$，$\alpha_2=(n_1,n_2)^T$都是实对称矩阵$S_x$的特征向量，所以满足：

$$S_x\alpha=\lambda\alpha$$
由于实对称矩阵的性质：

• 实对称矩阵一定可以被正交矩阵$U$对角化，这个正交矩阵由按照$\lambda$的降序排列对应的特征向量的组合：

$$U= (\alpha_1,\alpha_2)^T$$

$\lambda_1$对应的特征向量为$\alpha_1$,同理$\lambda_2$对应的特征向量为$\alpha_2$.

经过推导有：

$$YY^T=S_y=UXX^TU^T=US_xU^T$$
由于$S_x\alpha=\lambda\alpha$,所以有：

$$S_xU^T=U^T\left ( \begin{array}{c} \lambda_1 & 0\\ 0 & \lambda_2 \end{array} \right)$$

$$S_y= US_xU^T= \left ( \begin{array}{c} \lambda_1& 0 \\ 0 & \lambda_2 & \\ \end{array} \right)$$

又因为$S_y$为：

$$S_y= \left ( \begin{array}{c} s_{aa} & s_{ab} \\ s_{ba} & s_{bb} \end{array} \right)$$
所以$s_{ab}=s_{ba}=0$,$s_{aa}=\lambda_1$,$s_{bb}=\lambda_2$.可见，对于新的坐标系$ab$,坐标系的基底之间（即特征向量之间）满足相互正交.

3、总结一下前面的推导

总结一下前面的推导，对于二维平面的一堆没有明显规律的点集用$xy$坐标表示为$X$:

$$X= \left ( \begin{array}{c} x_1& x_2 & \cdots & x_n \\ y_1 & y_2 & \cdots & y_n \\ \end{array} \right)$$
求其协方差矩阵$S_x$：

$$S_x= \left ( \begin{array}{c} s_{xx} & s_{xy} \\ s_{yx} & s_{yy} \end{array} \right)$$

然后由于$S_x$是实对称矩阵，由于实对称矩阵有很多可爱的性质，让我们可以利用求$S_x$的特征值和特征向量，进而求出一个正交矩阵$U$,使得

$$S_y=US_xU^T=\left ( \begin{array}{c} s_{aa} & s_{ab} \\ s_{ba} & s_{bb} \end{array} \right) =\left ( \begin{array}{c} \lambda_1& 0 \\ 0 & \lambda_2 \\ \end{array} \right)$$

a坐标下方差$s_{aa}$,和$b$坐标下方差$s_{bb}$，满足$s_{aa}=\lambda_1$,$s_{bb}=\lambda_2$，且$ab$坐标相互正交.

这一切正是我们需要的，我们可以根据比较前几个坐标的方差大小，忽略方差不是那么大的坐标.就如最初的那个例子一样.

4、补充一点

我们说，求出的对角矩阵$S_y$:

$$S_y =\left ( \begin{array}{c} \lambda_1& 0 \\ 0 & \lambda_2 \\ \end{array} \right)$$
根据$\lambda_1和\lambda_2$的大小，忽略方差小的那一个，那么我们会想：$\lambda$一定是正数吗？

查阅线性代数的知识，有下面这个定义：

对于实对称矩阵，如果其特征值全部为正数，那么这个矩阵是正定的；如果其特征值非负，那么这个矩阵是半正定的。满足以上两种的任意一种的实对称矩阵，我们称之为非负定的。

显然，我们希望$S_y$是非负定的.

同时，线性代数理论中，有这样判断一个实对称矩阵是否为非负定的方法：

• $一个实对称矩阵是非负定的充要条件是该矩阵的顺序主子式\ge0$

对于二维情况：

$$S_x=\left ( \begin{array}{c} s_{xx} & s_{xy} \\ s_{yx} & s_{yy} \end{array} \right)$$
根据柯西不等式，有
$$s_{xy}^2\le s_{xx}s_{yy}$$
所以有二维情况下，$s_{xx}\ge0$,$s_{xx}s_{yy}-s_{xy}s_{yx}\ge0$,所以，满足顺序主子式都$\ge0$.

对于$n$维情况？

---

四、算法代码实现

1、数据集准备

集中在第一象限的点集

2、python代码

#!/usr/bin/env python
# coding=utf-8
from numpy import *
import matplotlib
import matplotlib.pyplot as plt
def loadDataSet(fileName, delim = "\t"):
    fr = open(fileName)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    datArr    = [map(float, line) for line in stringArr]#对每行数据的两个数字都变为float型
    return mat(datArr)
def pca(dataMat, topNfeat=9999999):                      #dataMat为1000×2的矩阵，
    meanVals    = mean(dataMat, axis=0)                  #numpy中的mat类型有算均值mean的方法,并且对每列进行计算，得到2*1的矩阵
    meanRemoved = dataMat - meanVals                     #计算协方差
    covMat      = cov(meanRemoved, rowvar = 0)           #计算协方差，维度为2,得到2×2的协方差矩阵
    eigVals, eigVects = linalg.eig(mat(covMat))          #计算特征值及特征向量，2×1,2×2,得到两个特征值，两列特征向量，向量长度为2
    eigValInd   = argsort(eigVals)                       #默认快排，从小到大进行排序
    eigValInd   = eigValInd[:-(topNfeat+1):-1]           #得到topNfeat个特征值，大的特征值排在前面，若topNfeat为1。，则取1个特征值
    redEigVects = eigVects[:, eigValInd]                 #将topNfeat个最大的特征值对应的那列向量取出来，若topNfeat为1,
    lowDDataMat = meanRemoved * redEigVects              #1000×2 * 2×topNfeat得到1000×topNfeat的低维空间的数据
    reconMat    = (lowDDataMat * redEigVects.T) + meanVals#1000×topNfeat * topNfeat*2 +1000×2重构数据
    return lowDDataMat, reconMat
dataMat = loadDataSet("testSet.txt")
lowDMat, reconMat = pca(dataMat, 1)
print shape(dataMat)
print shape(lowDMat)
print shape(reconMat)
fig = plt.figure()
ax  = fig.add_subplot(111)
ax.scatter(dataMat[:,0].flatten().A[0], dataMat[:,1].flatten().A[0], marker="^", s=90)
ax.scatter(reconMat[:,0].flatten().A[0], reconMat[:,1].flatten().A[0], marker="o", s=50, c="red")
plt.show()

3、可视化