PCA计算流程详解与实现（Python详细编码，全部测试正确，与sklearn完全一致，只有7行代码）

如果无法理解一种算法的原理，那么就整理出它的流程，再在直接的计算过程中慢慢体会它的思想。

– 我说的。

许多问题在一开始都很难理解其原理。像PCA这种经典的算法，其实从头至尾的每个步骤都非常严密，细节的知识特别多。但最近看到不少的书、网上的资料都是在通篇扯原理，讲道理，到最后就只列一个调包的程序，读完之后还是总觉得不清楚它到底是怎么玩的，这就没多大意思了。

理解原理是进行深入研究的必要前提，但要做到这一点的确很难。但这个其实不是我们的问题，而是因为经典的方法通常是被许多优秀的科学家和工程师不断完善才有了今天这副完美的样子。作为凡人的我们一时半会理解不透是非常正常的事情。

在这种时候最好的办法就是不去纠结过深的推导和证明（因为已经被证明是对的了），从复杂的理论中暂时抽离出来单纯地梳理其操作流程，在这些具体的计算过程中再去体会那些理论本身的意义。这是一个屡试不爽的方法。在这个过程中，我们需要自觉地避开那些与我们主要目的无关的信息，梳理出最需要的技术脉络。这自然需要的是阅读理解能力、知识转换能力（动手能力），长期坚持下去又可以不断提升这些方面的本事，形成一个非常良性又高效的循环。

本文就以PCA为例，展示一遍上述思路的具体操作步骤。

1、问题描述

注意，为了和程序数据结构匹配这里使用行向量。

原特征
$$\mathbf{x}=\left(x_1,x_2,\dots ,x_d\right)$$
新特征

$$\mathbf{w}=\left(w_1,w_2,\dots ,w_d\right)$$
记样本为：
$$\mathbf{X}=\left(\begin{array}{c}{\mathbf{x}}^1\\ {\mathbf{x}}^2\\ ⋮\\ {\mathbf{x}}^n\end{array}\right)$$
问题目标：

求出特征变换矩阵$A$使得：
$${\mathbf{X}}_{n\times d}{A}_{d\times d}={\mathbf{W}}_{n\times d}$$

2、计算过程

Step 1: 计算原特征协方差矩阵
$$C=(\mathbf{X}-\mathbf{u}{)}^T(\mathbf{X}-\mathbf{u})$$
Step 2: 计算协方差矩阵特征分解:
$$C\xi =\lambda \xi$$
Step 3: 特征值排序
$${\lambda }_1\ge {\lambda }_2\ge {\lambda }_d$$
Step 4: 对应特征向量排序
$${\xi }_1,{\xi }_2,\dots ,{\xi }_d$$
Step 5: 获得变换矩阵：
$$A=\left(\begin{array}{c}{\xi }_1\\ {\xi }_2\\ ⋮\\ {\xi }_d\end{array}\right)$$
Step 6: 计算特征贡献率(第$k$个新特征的贡献率):
$$\mathbf{r}=\frac{{\lambda }_k}{{\sum }_{i=1}^d{\lambda }_i}$$
回顾：若要计算前 $k$ 个主成分，则只需要取 $\mathbf{W}=\mathbf{X}A$ 的前 $k$ 列即可。

3、Python代码

先上代码，直接对照公式一步步来：

import numpy as np
import scipy.linalg as lina
x = np.random.rand(10,5) #随机生成一组样本
x -= x.mean(axis=0) 	# 见详注1
C = x.T.dot(x)    		# 计算自协方差矩阵
lam,v= lina.eig(C)    	# 特征分解，v是
new_index = np.argsort(lam)[::-1]    # 特征值排序，见详注2
A = -v[:,new_index]					 # 得到A
w = x.dot(A)						 # 计算变换后的特征
r = lam[new_index]/lam.sum()		 # 计算所有特征对应的贡献率

测试结果（只看关键结果）：

w[:,:2]		# 新特征的前2个

>>> array([[-0.3939524518, -0.4184678305],
       [-0.5907434013,  0.2033346207],
       [-0.4585388051, -0.111367225 ],
       [ 0.4552495673, -0.0405062598],
       [-0.2335902798, -0.4260334862],
       [ 0.4523182771,  0.039755097 ],
       [ 0.0902288594,  0.1869543779],
       [ 0.089419155 ,  0.7656098218],
       [ 0.7645053936, -0.3353675658],
       [-0.1748963144,  0.1360884499]])

r			# 各个特征值对应的贡献率

>>> array([0.4026073116, 0.2589988934, 0.2088275432, 0.0902665298,
       0.0392997221])

4、与sklearn结果对比

pca = PCA(n_components=2)

pca.fit(x)	# x还是最开始那个x

pca.explained_variance_ratio_

>>> array([0.4026073116, 0.2589988934])	# 前2个特征对应的贡献率（完全一致）

pca.transform(x)								# 降维变换（完全一致）

>>> array([[-0.3939524518, -0.4184678305],
       [-0.5907434013,  0.2033346207],
       [-0.4585388051, -0.111367225 ],
       [ 0.4552495673, -0.0405062598],
       [-0.2335902798, -0.4260334862],
       [ 0.4523182771,  0.039755097 ],
       [ 0.0902288594,  0.1869543779],
       [ 0.089419155 ,  0.7656098218],
       [ 0.7645053936, -0.3353675658],
       [-0.1748963144,  0.1360884499]]) 

再仔细计算一下其相差的值：

 w1 = w[:,0]
 w2 = pca.transform(c)[:,0]
((w1-w2)**2).sum()
>>> 2.2980196096428498e-30

相差只有1e-30的量级，注意这里是平方和，那么完全可以视作只是一些系统误差（包括最后一位的舍入误差）造成的。

再看r的值：

pca.explained_variance_ratio_ - r[:2]
>>> array([ 1.1102230246e-16, -1.6653345369e-16])

仍然只有1e-16的量级。

因此上述方法和sklearn中的方法完全一致。

5、详注

• 详注1：x -= x.mean(axis=0)

这里x.mean(axis=0) 表示求出x中每列的平均值，返回一个一维数组。这里之所以可以让不同形状的数组做减法是用到了python自带的broadcasting机制（广播机制），它会自动将一维数组扩充至二维，使其变成每一行都完全相同的矩阵，行数与x相同。因此不需要单独再做其它处理。

• 详注2:

这里首先要明白一点的是，利用eig函数得到的特征值和特征向量是未排序的，同时特征值与特征向量的位置对应，也就是说特征值的下标与特征向量第二维的下标（也可以理解成列标）是完全一样的。因此只需要对特征值进行排序，同时获得其降序排列后新的下标位置即可。

其次，np.argsort只返回按升序排列得到的新下标集。那么如果对新下标集进行反向处理，即可得到降序排列的新下标。（有点绕，但这个逻辑是层次清楚的）。而符号[::-1]正好是一维数组倒序，因此所有操作放在一起就有了这段代码new_index = np.argsort(lam)[::-1]。