深入理解PCA(主成分分析法)算法

问题引入

现有某商店分析影响销售额的因素有哪些，做了以下调查。

日期${\text{x}}_1$	2	6	7	9
空气质量${\text{x}}_2$	10	60	30	90
上班时间${\text{x}}_3$	8	8	9	7

假装目前这三个影响因素太多，需要降维处理。
原始思想
如果要降到2维，最直接的做法就是将其中一个因素${\text{x}}_i$去除，然后分析剩余的两个，但是这种做法过于粗暴，是否存在一种更加合理的方法去处理。

---

符号声明

${\text{x}}_i$代表表格中第i行，${\text{x}}_j$代表表格中第j列，$x_{ij}$代表矩阵中第(i,j)个元素，其中$i=1,2,\cdots ,m;j=1,2,\cdots ,n$。

基本思想

方案一
求出每个因素${\text{x}}_i$的方差
方差公式$$s_i^2=\frac{1}{n-1}\sum _{j=1}^n(x_{ij}-{\overline{\text{x}}}_i{)}^2$$
从概率上来说方差可以反应数值的分散程度，对于分散程度较小的因素而言(极端情况为一个不变的数值)说明分析的价值不大，所以可以舍弃。(ps:有人从熵的角度去解释，但是本人目前没有接触该方面的知识，以后会补上)
方案二
但是由于各类因素是充满联系的，比如说空气质量可能会影响上班时间，一些方差小的原因可能不是该因素本身的原因，而是其它因素影响造成的，需要将这些因素之间的影响去除掉，在除掉之前先引进协方差的概念。
向量${\text{x}}_i$与向量${\text{y}}_i$之间的协方差为
$$cov({\text{x}}_i,{\text{y}}_i)=\frac{1}{n-1}\sum _{j=1}^n(x_{ij}-{\overline{x}}_i)(y_{ij}-{\overline{y}}_i)$$

为方便求解，将${\text{x}}_i$与${\text{y}}_i$分别进行去中心化（将向量中的每个元素减去平均值）处理
$${\text{x}}_i={\text{x}}_i-{\overline{\text{x}}}_i$$
$${\text{y}}_i={\text{y}}_i-{\overline{\text{y}}}_i$$
此时，协方差公式简化为
$$cov({\text{x}}_i,{\text{y}}_i)=\frac{1}{n-1}\sum _{j=1}^n{\text{x}}_{ij}{\text{y}}_{ij}=\frac{1}{n-1}{\text{x}}_i^T{\text{y}}_i$$
协方差可以反应两个向量之间的相关性：
当cov(X,Y)<0时，两者负相关。
当cov(X,Y)=0时，两者不相关。
当cov(X,Y)>0时，两者正相关。
ps:相关性与独立性不是一个概念。
不同因素之间的相关性可以利用协方差矩阵$\text{C}$来表示，具体公式如下：
$$\text{C}=\left[\begin{array}{cccc}cov(x_1,x_1)& cov(x_1,x_2)& \cdots & cov(x_1,x_n)\\ cov(x_2,x_1)& cov(x_2,x_2)& \cdots & cov(x_2,x_n)\\ ⋮& ⋮& \ddots & ⋮\\ cov(x_m,x_1)& cov(x_m,x_1)& \cdots & cov(x_m,x_m)\end{array}\right]=\frac{1}{n-1}\sum _{j=1}^n{\text{x}}_j{\text{x}}_j^T$$
对于问题中的三个变量之间的关系，可以利用协方差矩C来表达(等式右边很关键，可以自行推导或验证)，如下：
$$C=\left[\begin{array}{ccc}cov(x_1,x_1)& cov(x_1,x_2)& cov(x_1,x_3)\\ cov(x_2,x_1)& cov(x_2,x_2)& cov(x_2,x_3)\\ cov(x_3,x_1)& cov(x_3,x_2)& cov(x_3,x_3)\end{array}\right]=\frac{1}{4-1}\sum _{j=1}^4{\text{x}}_j{\text{x}}_j^T$$
协方差矩阵对角线位置上的数值表示方差，非对角线位置上的数值表示各因素之间的相关性，假设存在一个协方差矩阵D为对角矩阵，形式如下
$$D=\left[\begin{array}{ccc}cov(y_1,y_1)& 0& 0\\ 0& cov(y_2,y_2)& 0\\ 0& 0& cov(y_3,y_3)\end{array}\right]=\frac{1}{4-1}\sum _{j=1}^4{\text{y}}_j{\text{y}}_j^T$$
此时就可以按照改进一中的方案利用方差的大小进行降维。
为使各个因素之间不相关，需要一种方法对各个因素进行坐标转换使协方差矩阵转换为对角矩阵 ***（记住这一求解目标）***。可以看到协方差矩阵C是一个对称矩阵，线性代数的知识表明，可以将C进行坐标变换，使各个元素在新的坐标系中不相关，下面进行这一方面的推导，这个推导过程实际上就是PCA的产生过程。

---

投影

向量X与向量Y的点乘公式如下
$$\text{X}\cdot \text{Y}=||\text{X}||||\text{Y}||cos<\text{X},\text{Y}>$$
以二维平面为例

向量X在向量Y上投影的长度可以表示为
$$||\text{X}||cos<\text{X},\text{Y}>$$
cos<X,Y>的大小仅仅与X与Y的夹角有关，所以X在Y上的投影可以表达为
$$\text{X}\cdot {\text{e}}_y={\text{e}}_y^T\text{X}$$
其中${\text{e}}_y$表示与Y同方向的单位向量。

坐标变换

坐标系上点的坐标是该点在各个坐标中上的投影，利用这一概念对坐标系中的点进行坐标变换。
以二维坐标为例，$e_1$与$e_2$是平面坐标系中的两个正交向量(也称为正交基)，点a在坐标系中的位置为$(1,1{)}^T$，求点a在$e_1$与$e_2$组成的坐标系中的位置。

首先对$e_1$与$e_2$单位化处理
$$e_1=(\frac{2}{\sqrt{5}},\frac{1}{\sqrt{5}}{)}^T$$
$$e_2=(-\frac{1}{\sqrt{5}},\frac{2}{\sqrt{5}}{)}^T$$
此时$e_1与e_2$是单位向量且正交，像这样一组为单位向量且两两正交的向量被称为标准正交基。
利用投影可以求出点P在新坐标中的位置。
a在$e_1$上的投影距离为
$$e_1^{T}a=\frac{3}{\sqrt{5}}$$
a在$e_2$上的投影距离为
$$e_2^{T}a=\frac{1}{\sqrt{5}}$$
写成矩阵表达式为
$$\left[\begin{array}{c}{e}_1^T\\ e_2^T\end{array}\right]a$$
a在新的坐标下的坐标为$(\frac{3}{\sqrt{5}},\frac{1}{\sqrt{5}})$。
在m维空间中，令$\text{P}=\left[\begin{array}{cccc}{e}_1& e_2& \cdots & e_m\end{array}\right]$，求$\text{x}$在$\text{P}$组成的m维标准正交基中的坐标公式如下
$${\left[\begin{array}{cccc}{e}_1& e_2& \cdots & e_m\end{array}\right]}^T\text{x}={\text{P}}^T\text{x}$$

---

证明过程

我们的目标是通过坐标变换的方式找到一个为对角矩阵协方差矩阵，为证明这一目标，不妨假设存在左边转换${\text{y}}_j={\text{P}}^T{\text{x}}_j$使协方差矩阵C转换为对角矩阵D，推导过程如下
$$\begin{array}{c}\begin{array}{rlr}D& =& \frac{1}{n-1}\sum _{j=1}^n{\text{y}}_j{\text{y}}_j^T\\ & =& \frac{1}{n-1}\sum _{j=1}^n{\text{P}}^T{\text{x}}_j{\text{x}}_j^T\text{P}\\ & =& \frac{1}{n-1}{\text{P}}^T(\sum _{j=1}^n{\text{x}}_j{\text{x}}_j^T)\text{P}\\ & =& \frac{1}{n-1}{\text{P}}^{T}C\text{P}\end{array}\end{array}\text{\hspace{1em}(推导1)}$$
ps:不知道为什么\begin{align}\end{align}报错，无法左对齐，排版就先凑合看吧。

接下来回顾一下对称矩阵的性质:
1.设$\text{C}$为是n阶实对称矩阵，则必有正交矩阵$\text{P}$，使${\text{P}}^{-1}\text{C}\text{P}={\text{P}}^T\text{C}\text{P}$为对角矩阵$\text{D}$，且D中主对角线上的元素为$\text{C}$的特征值，它们的排列次序与对应于他们的特征向量在$\text{P}$中的排列次序一致。
2.实对称矩阵的k重特征值所对应的线性无关的特征向量恰好有k个。
3.实对称矩阵的相异特征值所对应的特征向量是正交的。
由于C为对称矩阵，根据对称矩阵的性质则一定存在矩阵P使${\text{P}}^T\text{C}\text{P}$变为对角矩阵，改进二中提出的目标得以证明。
接下来只需要保留方差较大的${\text{y}}_i$即可，由上面的推导(1)以及性质1可知，坐标转换之后的协方差矩阵$\text{D}$是一个对角矩阵，并按照方差(根据性质3以及推导(1)可以看出方差就是$\text{C}$的特征值)从大向小的顺序排列如下:
$${\text{D}}_{m\times m}=\left[\begin{array}{cccc}{\sum }_{j=1}^n{\text{y}}_{1j}^2& 0& \cdots & 0\\ 0& {\sum }_{j=1}^n{\text{y}}_{2j}^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ \cdots & \cdots & {\sum }_{j=1}^n{\text{y}}_{kj}^2& \cdots \\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sum }_{j=1}^n{\text{y}}_{mj}^2\end{array}\right]$$
为方便描述定义转换坐标后的${\text{y}}_i$分量所对应的原来的基向量(也叫特征向量)为${\text{e}}_i$，那么$${\text{P}}^T={\left[\begin{array}{c}{\text{e}}_1,{\text{e}}_2,\cdots ,{\text{e}}_m\end{array}\right]}^T$$
为保留前k个方差较大的$\text{y}$分量，由$\text{D}$知前k行对应前k个${\text{y}}_i$分量，只需要保留分量${\text{y}}_1,\cdots ,{\text{y}}_k$即可。
由坐标转换公式知
$${\text{y}}_j={\text{P}}^T{\text{x}}_j={\left[\begin{array}{c}{\text{e}}_1,{\text{e}}_2,\cdots ,{\text{e}}_k,\cdots ,{\text{e}}_m\end{array}\right]}^T{\text{x}}_j$$
由于
$${\text{y}}_j=\left[\begin{array}{c}{y}_{1j},y_{2j},\cdots ,y_{kj},\cdots ,y_{mj}\end{array}\right]$$
现需要保留分量${\text{y}}_1,\cdots ,{\text{y}}_k$，所以直接取前k个元素$y_{1j},y_{2j},\cdots ,y_{kj}$即可。
为了方便计算通常将较小的方差(特征值)对应的特征向量${\text{e}}_{k+1},\cdots ,{\text{e}}_m$剔除，计算公式如下
$${\text{y}}_j={\text{P}}^T{\text{x}}_j={\left[\begin{array}{c}{\text{e}}_1,{\text{e}}_2,\cdots ,{\text{e}}_k\end{array}\right]}^T{\text{x}}_j$$

---

算法步骤

最后总结将PCA降维步骤如下:
1).将${\text{x}}_i$去中心化
$${\text{x}}_i={\text{x}}_i-{\overline{\text{x}}}_i$$
2).计算协方差矩阵C。
$$\text{C}=\frac{1}{n-1}\text{x}{\text{x}}^T$$
实际上系数$\frac{1}{n-1}$仅仅涉及到最终计算出来的方差的倍数，并不影响方差之间的比例以及大小，所以去掉也不会对结果有任何影响，为了减少计算量常常写作
$$\text{C}=\text{x}{\text{x}}^T$$
3).求出C的全部相异特征值${\lambda }_1,\cdots ,{\lambda }_{\dot{m}}$;（有些特征值可能相等，$\dot{m}\le m$）
4).对于每一个重特征值${\lambda }_i$,求出对应的$r_i$个线性无关的特征向量${\alpha }_{i1},{\alpha }_{i2},\cdots ,{\alpha }_{i{r}_i}(i=1,2,\cdots ,\dot{m})$,由性质2知${\sum }_{i=1}^{\dot{m}}{r}_i=m$;(说明正交矩阵P是mxm的)
5).利用施密特正交化方法，把对应于每一个${\lambda }_i$的线性无关特征向量先正交化再单位化，得到一组等价的两两正交的单位向量组${\eta }_{i1},{\eta }_{i2},\cdots ,{\eta }_{i{r}_i}(i=1,2,\cdots ,\dot{m})$,他们仍为矩阵$\text{P}$对应于${\lambda }_i$的特征向量。(很多博客没有正交化步骤，通过性质3可以知道如果出现相等的特征值，所对应的特征向量不一定正交)
6).将上面求得的正交单位向量作为列向量，排成一个n阶方阵$\text{P}$，则$\text{P}$即为所求的正交矩阵，此时,$${\text{P}}^T\text{C}\text{P}=\text{D}$$
为对角矩阵，D中主对角线上的元素为$\text{C}$的特征值，且它们的排列次序与对应于他们的特征向量在$\text{P}$中的排列次序一致，一般情况下总是按照特征值降序排列特征向量。(很多博客没讲清楚降维时保留较大的特征值对应的特征矩阵的原因，通过改进二我们知道，这是为了保留较大的方差，$\text{D}$的对角线为坐标转换后各个影响因素的所对应的方差，同时对角线元素也是矩阵$\text{C}$的特征值，所以保留较大的特征值就是保留较大的方差)。
7).保留特征值较大的特征向量，不知道降低到多少维合适时一般按照经验判断保留到特征值所在比例为95%左右，有时候为了可视化比例可能远低于这个数值，令$\text{P}=\left[\begin{array}{cccc}{\text{e}}_1& {\text{e}}_2& \cdots & {\text{e}}_k\end{array}\right](k\le m)$，最终降维之后的坐标为${\text{y}}_j={\text{P}}^T{\text{x}}_j$。(舍弃的${\text{e}}_{m+1},\cdots ,{\text{e}}_n$是变换坐标之后的一些不重要的坐标轴，饶了一大圈还是回到了原始思想的粗暴方法中，也就是说如果原始数据的协方差矩阵符合矩阵D的形式，那么按照PCA的思想可以直接去除方差较小的因素)
PCA降维应对大量现实问题时比较合理，但是通过原始思想也可以看出，当求出的特征值(矩阵D对应的方差)变化不明显时，PCA表现效果不会特别如愿。

---

python代码

现在来处理文章开篇的引出的问题

附加python代码

import numpy as np 
def myPCA(a,k):
	print('调查表矩阵：')
	print(a)
	assert(isinstance(k,int))
	assert(isinstance(a,np.ndarray))
	ax,ay=a.shape
	assert(k>0 and k<=ay)
	arr_mean=np.array([np.mean(a,axis=1)]).T
	print('每行行平均值：')
	print(arr_mean)
	a=a-arr_mean#去中心化
	print('去中心化后矩阵：')
	print(a)
	#计算协方差矩阵
	arr_c=np.dot(a,a.T)
	print('协方差矩阵：')
	print(arr_c)
	#计算特征值与特征向量,numpy计算的特征向量都是单位化之后的
	eig_value,eig_vec=np.linalg.eig(arr_c)
	print('特征值：')
	print(eig_value)
	print('特征向量，按行排列：')
	print(eig_vec.T)
	#升序排列
	index=np.argsort(eig_value)
	#取最后k个特征向量
	eig_vec_k=eig_vec[:,index[:-k-1:-1]].T
	print('降到%d维后：'%k)
	return np.dot(np.array(eig_vec_k),a)

k=2
a = np.array([[2,6,7,9]
			,[10,60,30,60]
			,[8,8,9,7]])  
print(myPCA(a,k))

调用scikit-learn求解结果以及代码

第一行与自己写的代码中的结果有着正负上的区别，熟悉线性代数的都知道，一个特征值的单位特征向量有正负两个，两者都是对的。

import numpy as np
from sklearn.decomposition import PCA
a = np.array([[2,6,7,9]
			,[10,60,30,60]
			,[8,8,9,7]])
#n_components降低到的维度
pca=PCA(n_components=2)
#最终结果
print(pca.fit_transform(a.T).T)