主成分分析(Principle Component Analysis,PCA),是一种多变量统计分析方法,也是数据降维和可视化的常用方法。PCA的原理是将原矩阵投影到一个新的正交坐标系下,且通过依次选择新坐标轴的方向,使得矩阵在新坐标轴上的投影(主成分的样本值)的方差最大[1]。
本文以粗体表示矩阵和向量,标量则不加粗。
<关于向量的投影>
设向量
![]()
为单位向量,则向量
![]()
在向量
![]()
上的标量投影(scalar projection)为
![]()
或
![]()
,向量投影(vector projection)为
![]()
。
![]()
与其向量投影的差(vector rejection)为
![]()
。
https://en.wikipedia.org/wiki/Vector_projection
PCA模型有多种方式可以推导,现介绍其中三种:特征分解,奇异值分解,非线性迭代偏最小二乘法(NIPALS)。
1. 特征分解
设矩阵
![]()
为
![]()
行
![]()
列的矩阵,即包含
![]()
个观测点(样本),
![]()
个变量(特征)。通常要先对
![]()
进行预处理:去平均值,使其每列(每个特征)的均值为0,并缩放使其方差为1。预处理之后计算方差更简单。
设
![]()
已经过预处理,现在将
![]()
中的所有样本点都投影到单位向量
![]()
(要求的第一个坐标轴,或者说权重向量,weight vector)上,得到第一个主成分的样本值
![]()
。注意
![]()
的每一行即为前文中的
![]()
,所以
![]()
向
![]()
投影得到向量
![]()
。
![]()
的每个元素为
![]()
中相应样本点的score,是
![]()
个变量的线性组合,系数(权重)即为
![]()
各元素 。
现在要选择
![]()
的方向,使得
![]()
的方差最大。不难证明,若矩阵
![]()
中每列的均值为0,则
![]()
中各元素的均值也为0,故
![]()
的方差可简单地表示为
![]()
,所以求解
![]()
的过程可形成如下有约束的优化问题:
其中
![]()
,故优化问题变为
用拉格朗日乘子法求解,拉格朗日函数为:
将
![]()
对
![]()
求梯度,对
![]()
求偏导,并令其为0,得
由(1)式得
![]()
。这表明,极值点处
![]()
为实对称矩阵
![]()
的一个特征向量,而
![]()
为对应的特征值。故目标函数的极值为:
即
![]()
对应的特征值
![]()
。既然我们要求目标函数的最大值,所以
![]()
应取
![]()
的最大特征值(主特征值),而第一个坐标轴
![]()
即为
![]()
的最大特征值对应的特征向量(主特征向量)。
(注:
![]()
为
![]()
的方阵,正比于协方差矩阵
![]()
;且
![]()
为半正定实对称矩阵,其特征值非负。)
接下来计算第二个坐标轴
![]()
,并要求
![]()
和
![]()
正交,即
![]()
。将
![]()
向
![]()
投影,得到第二个主成分
![]()
。此处也可理解为将
![]()
中第一个主成分
![]()
去除后,再向
![]()
投影,即
![]()
,毕竟
![]()
和
![]()
正交,结果是一样的。
求解
![]()
的过程同样可以形成如下一个有约束的优化问题:
用拉格朗日乘子法求解,拉格朗日函数为:
将
![]()
对
![]()
求梯度,对
![]()
和
![]()
求偏导,并令其为0:
将(2)式中的向量
![]()
与
![]()
求点积:
计算得
![]()
,再代回(2)式,有
上式表明,极值点处
![]()
为
![]()
的一个特征向量,
![]()
为对应的特征值。化简得目标函数极值
![]()
。而既然我们要求目标函数的最大值,故应取矩阵
![]()
的第二大的特征值。
依此类推,求出矩阵
![]()
的特征值并由大到小排列,依次取其对应的特征向量作为坐标轴即可。这实际上是对
![]()
进行特征分解的过程。
假设一共提取出
![]()
个主成分,放在一个矩阵里(即
![]()
在新坐标系下的坐标),记为:
![]()
;将
![]()
个坐标轴(权重向量)放在一个矩阵里:
![]()
,则
![]()
,residual矩阵为
![]()
。
* 多重特征值的情况也适用。由于
![]()
为半正定实对称矩阵,其特征值非负,且不同特征值对应的特征向量相互正交;若存在
![]()
重特征值,则必存在
![]()
个线性无关的特征向量,可通过施密特正交化产生相互正交的特征向量,这些正交化后的特征向量即可作为新坐标系的坐标轴。
2. 奇异值分解(SVD)
根据前面的推导,PCA模型实际上是一个特征值问题。注意到
![]()
的特征值即为
![]()
的奇异值的平方,那么自然地想到,可以用矩阵的奇异值分解(singular value decomposition, SVD)来得到权重向量。
对矩阵
![]()
进行SVD分解:
![]()
,则
![]()
的主对角元素为
![]()
的奇异值,其平方即为
![]()
的特征值;而
![]()
的列则是
![]()
的特征向量,也就是我们要找的权重向量[2]。
https://en.wikipedia.org/wiki/Singular_value_decomposition
3. NIPALS与Power Method
既然要降维,我们往往不需要计算出
![]()
的全部特征向量,因为如果取全部特征向量作为新坐标轴,实际上只是将原矩阵进行了旋转而已。故只需取最大的几个特征值对应的特征向量即可。非线性迭代偏最小二乘法(Non-linear iterative partial least-squares, NIPALS)正是这样一种迭代的算法[1]。NIPALS可以看作是幂法(Power Method)的一种变体,在计算得到一个主成分后,进行矩阵收缩(matrix deflation),然后计算下一个主成分。
下面先简单介绍一下幂法。设
![]()
为
![]()
的方阵,
![]()
为其
![]()
个特征值中模最大的特征值,即:
则称
![]()
为
![]()
的主特征值(dominant eigenvalue),其对应的特征向量称为
![]()
的主特征向量(dominant eigenvector)。
假设
![]()
的
![]()
个特征向量
![]()
线性无关,
![]()
。
则任意选择一非零向量
![]()
,
![]()
可以表示为
![]()
的线性组合:
上式左乘
![]()
,得
左乘
![]()
次后,
当
![]()
足够大时,
![]()
。
为避免溢出(overflow/underflow),往往进行归一化处理(normalize)。记上式归一化处理后得到的结果为
![]()
,则
![]()
(如果
![]()
为单位向量)。
也就是说当迭代收敛后,得到的
![]()
为矩阵
![]()
的主特征向量。至于主特征值,可由Rayleigh商求出[3]:
![]()
。
不难看出,在PCA中要计算
![]()
的主特征值,将上面的
![]()
用
![]()
替换即可,因为
![]()
存在
![]()
个线性无关的特征向量。不过,若最大特征值为多重特征值,则可能会得到其对应的特征向量的线性组合。
以上是幂法的原理。下面介绍用NIPALS算法计算PCA主成分的流程:
0) 对
![]()
预处理;
1) 对
![]()
赋初值,可以是除零矢量以外的任意初值;文献中也常将其取作
![]()
的一列;
2) 计算
![]()
;(可看作将
![]()
的每一列对
![]()
回归)
3) 归一化
![]()
;
4) 计算
![]()
;(可看作将
![]()
的每一行对
![]()
回归)
5) 检查
![]()
是否收敛(两次迭代之间差的2范数);
如果已收敛,则得到
![]()
为权重向量,
![]()
为
![]()
在
![]()
上投影得到的主成分(样本值),接下来进行第6步;
否则回到第2步,继续迭代;
6) 矩阵收缩,记(residual)
![]()
;
令
![]()
,回到第1步,计算下一个主成分。
这个过程中,第2~4步的实际上是在循环执行
![]()
,然后将
![]()
归一化;这就是幂法中不断左乘矩阵的过程。因此根据幂法,最终
![]()
会收敛到
![]()
的主特征向量。
注意到这个过程中,也在循环执行
![]()
,这表明
![]()
正比于
![]()
的主特征向量。事实上,若
![]()
,则
![]()
。
如前文所说,幂法会收敛到主特征向量,这意味着当NIPALS算法第一次收敛后,只能找到第一对主成分/权重向量;那么接下来如何计算第二对呢?这其中的关键是第6步:矩阵的收缩。
假设第一次收敛后得到了
![]()
和
![]()
,residual
![]()
。在NIPALS的下一次收敛后,会计算出矩阵
![]()
的主特征向量;然而我们的目标是计算
![]()
的第二大特征值对应的特征向量。所以现在的问题是:
![]()
的主特征向量是否和
![]()
的 第二大特征值对应的特征向量相同呢?
答案是肯定的,证明如下:
将
![]()
代入化简得
记
![]()
的
![]()
个特征向量为
![]()
,
![]()
;则有
若
![]()
,则
![]()
,故
![]()
;
若
![]()
,则
![]()
,故
![]()
;
这表明,
![]()
的特征向量也是
![]()
的特征向量。除了最大特征值用0替代以外,
![]()
的其他特征值均与
![]()
相同。这种矩阵收缩的方法叫Hotelling‘s deflation。
https://wiki.math.uwaterloo.ca/statwiki/index.php?title=deflation_Methods_for_Sparse_PCA
http://www.robots.ox.ac.uk/~sjrob/Teaching/EngComp/ecl4.pdf
至此我们知道,当NIPALS算法第二次收敛后,得到的向量
![]()
对应于
![]()
的第二大特征值对应的特征向量,也就是第二个权重向量,而
![]()
则为第二个主成分。
依此类推,可以用NIPALS逐个计算主成分。至于需要取多少个主成分,可以用交叉验证(cross-validation)的方式和
![]()
值来确定。
参考文献
[1] Kevin Dunn. Process Improvement Using Data. (Chapter 6, 6.5)
Process Improvement Using Data
[2] Jolliffe I.T. Principal Component Analysis, Second Edition, Springer Series in Statistics. Springer, NY, 2002. (Chapter 3, 3.5)
[3] Uri M. Ascher, Chen Greif. A First Course in Numerical Methods, SIAM, 2011. (Chapter 8, 8.1).
[4] https://en.wikipedia.org/wiki/Principal_component_analysis