利用主成分分析对声速剖面降维_数据降维: 主成分分析(PCA)与NIPALS算法数学原理...

主成分分析(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