主成分分析|PCA算法大全

主成分分析|PCA算法大全

主元分析提供了一种用低维数据来表示高维复杂数据最主要特征的途径。PCA的思想是将$n$维特征映射到$k$维空间上$k<n$，这$k$维特征是全新的正交特征，是新构造出来的$k$维特征，而不是简单地从$n$维特征中去除其余$n-k$维特征。在原始PCA的基础上，根据数据的特点，研究人员提出各种版本的PCA，但是核心思想一样。例如，针对非线性数据的核PCA，针对序列数据的CCIPCA，针对二维图片数据的2DPCA、2D2DPCA以及BDPCA等。下面首先详细介绍PCA的思想，然后简单介绍其他版本的PCA算法。

1. PCA原理

（本部分参考从零开始实现主成分分析(PCA)算法 ）
PCA可以被定义为数据在低维线性空间上的正交投影，这个线性空间被称为主子空间（principal subspace），使得投影数据的方差被最大化（Hotelling, 1933），即最大方差理论。等价地，它也可以被定义为使得平均投影代价最小的线性投影，即最小误差理论。平均投影代价是指数据点和它们的投影之间的平均平方距离（Pearson, 1901）。

1.1 最大方差理论

在信号处理中，认为信号具有较大的方差，噪声有较小的方差，信噪比就是信号与噪声的方差比，越大越好。因些我们认为，最好的$k$维特征是将$n$维样本点变换为$k$后，每一维上的样本方差都尽可能的大。

首先，考虑在一维空间$(k=1)$上的投影。我们可以使用$n$维向量$u$定义这个空间的方向。为了方便并不失一般性，我们假定选择一个单位向量，从而$u^{T}u=1$。

(假设数据是零均值化后的)
如上图所示，红色点表示原样本点$x^{(i)}$，$u$是蓝色直线的斜率也是直线的方向向量，而且是单位向量，直线上蓝色点表示原样本点$x^{(i)}$在$u$上的投影。容易知道投影点离原点的距离是$x^{(i)T}u$,由于这些原始样本点的每一维特征均值都为0，因此投影到$u$上的样本点的均值仍然是0。

假设原始数据集为$X_{m\times n}$，我们的目标是找到最佳的投影空间$W_{n\times k}=(w_1,w_2,...,w_k)$，其中$w_i$是单位向量且$w_i$是单位向量且$w_i$与$w_j(i\ne j)$正交，何为最佳的$W$？就是原始样本点投影到$W$上之后，使得投影后的样本点方差最大。

由于投影后均值 为0，因些投影后的总方差为：$$\frac{1}{m}\sum _{i=1}^m(x^{(i)T}w{)}^2=\frac{1}{m}\sum _{i=1}^m{w}^T{x}^{(i)}{x}^{(i)T}w=\sum _{i=1}^m{w}^T(\frac{1}{m}{x}^{(i)}{x}^{(i)T})w$$

其中$\frac{1}{m}{x}^{(i)}{x}^{(i)T}$就是原始数据集$X$的协方差矩阵（因为$x^{(i)}$的均值为0，因为无偏估计的原因，一般协方差矩阵除以$m-1$，这是用m）。

记$\lambda =\frac{1}{m}{\sum }_{i=1}^m(x^{(i)T}w{)}^2,\Sigma =\frac{1}{m}{x}^{(i)}{x}^{(i)T}$,则有$$\lambda =w^T\Sigma w$$

上式两边同时左乘$w$，注意到$w{w}^T=1$（单位向量），则有$$\lambda w=\Sigma w$$
所以$w$是矩阵$\Sigma$的特征值所对应的特征向量。

欲使投影后的总方差最大，即$\lambda$最大，因此最佳的投影向量$w$是特征值$\lambda$最大时对应的特征向量，因些，当我们将$w$设置为与具有最大的特征值向量相等时，方差会达到最大值。这个特征向量被称为第一主成分。

我们可以用一种增量的方式定义额外的主成分，方法为：在所有与那些已经考虑过的方向正交的所有可能的方向中，将新的方向选择为最大化投影方差的方向。如果我们考虑$k$维投影空间的一般情形，那么最大化投影数据方差的最优线性投影由数据协方差矩阵$\Sigma$的$k$个特征向量$w_1,...,w_k$定义，对应于$k$个最大的特征值 ${\lambda }_1,...,{\lambda }_k$。可以通过归纳法很容易地证明出来。

因此，我们只需要对协方差矩阵进行特征值分解，得到的前$k$大特征值对应的特征向量就是最佳的$k$维新特征，而且这$k$维新特征是正交的。得到前$k$个$u$以后，原始数据集$X$通过变换可以得到新的样本。

1.2 最小平方误差理论

如上图所示，假设有这样的二维样本点（红色点），按照前文我们讲解的最大方差理论，我们的目标是是求一条直线，使得样本点投影到直线或者平面上的点的方差最大。本质是求直线或者平面，那么度量直线求的好不好，不仅仅只有方差最大化的方法。再回想我们最开始学习的线性回归等，目的也是求一个线性函数使得直线能够最佳拟合样本点，那么我们能不能认为最佳的直线就是回归后的直线呢？回归时我们的最小二乘法度量的是样本点到直线的坐标轴距离。比如这个问题中，特征是$x$，类标签是$y$。回归时最小二乘法度量的是距离$d$。如果使用回归方法来度量最佳直线，那么就是直接在原始样本上做回归了，跟特征选择就没什么关系了。

因此，我们打算选用另外一种评价直线好坏的方法，使用点到直线的距离$d\prime$来度量。

现在有$m$个样本点$x^{(1)},...,x^{(m)}$，每个样本点为$n$维。将样本点$x^{(i)}$在直线上的投影记为$x^{(1{)}^{\prime \prime }}$，那么我们就是要最小化$$\sum _{i=1}^m(x^{(i{)}^{\prime }}-x^{(i)}{)}^2$$

这个公式称作最小平方误差(Least Squared Error)。

首先，我们确定直线经过的点，假设要在空间中找一点$x_0$来代表这$m$个样本点，“代表”这个词不是量化的，因此要量化的话，我们就是要找一个$n$维的点$x_0$，使得$$J_0(x_0)=\sum _{i=1}^m(x^{(i{)}^{\prime }}-x^{(i)}{)}^2$$最小。其中$J_0(x_0)$是平方错误评价函数。假设$\overline{x}$为$m$个样本点的均值，即$$\overline{x}=\frac{1}{m}\sum _{i=1}^m{x}^{(i)}$$
则$$\begin{gathered} J_0(x_0)=\sum _{i=1}^m(x_0-x^{(i)}{)}^2 \\ =\sum _{i=1}^m((x_0-\overline{x})-(x^{(i)}-\overline{x}){)}^2 \\ =\sum _{i=1}^m(x_0-\overline{x}{)}^2-2\sum _{i=1}^m(x_0-\overline{x}{)}^T(x^{(i)}-\overline{x})+\sum _{i=1}^m(x^{(i)}-\overline{x}{)}^2 \\ =\sum _{i=1}^m(x_0-\overline{x}{)}^2-2(x_0-\overline{x}{)}^T\sum _{i=1}^m(x^{(i)}-\overline{x})+\sum _{i=1}^m(x^{(i)}-\overline{x}{)}^2 \\ =-\sum _{i=1}^m(x_0-\overline{x}{)}^2+\sum _{i=1}^m(x^{(i)}-\overline{x}{)}^2 \end{gathered}$$

显然，上式的第二项与$x_0$无关，因此，$J_0(x_0)$在$\overline{x}$处有最小值。

接下来，我们确定直线的方向向量。我们已经知道直线经过点$\overline{x}$，假设直线的方向是单位向量$\overrightarrow{e}$。那么直线上任意一点$x^{(i{)}^{\prime }}$有：$x^{(i{)}^{\prime }}=\overline{x}+a_i\overrightarrow{e}$,其中，$a_i$是$x_{(i{)}^{\prime }}$到点$\overline{x}$的距离。

我们重新定义最小平方误差：

这时候点都聚集在新的坐标轴周围，因为我们使用的最小平方误差的意义就在此。另外，PRML书上从线性子空间的角度进行了详细的阐述，有兴趣的读者可以看看。

优点：

它是无监督学习，完全无参数限制的。在PCA的计算过程中完全不需要人为的设定参数或是根据任何经验模型对计算进行干预，最后的结果只与数据相关，与用户是独立的。

用PCA技术可以对数据进行降维，同时对新求出的“主元”向量的重要性进行排序，根据需要取前面最重要的部分，将后面的维数省去，可以达到降维从而简化模型或是对数据进行压缩的效果。同时最大程度的保持了原有数据的信息。

各主成分之间正交，可消除原始数据成分间的相互影响。

计算方法简单，易于在计算机上实现。

缺点：

如果用户对观测对象有一定的先验知识，掌握了数据的一些特征，却无法通过参数化等方法对处理过程进行干预，可能会得不到预期的效果，效率也不高。

贡献率小的主成分往往可能含有对样本差异的重要信息。

特征值矩阵的正交向量空间是否唯一有待讨论。

在非高斯分布的情况下，PCA方法得出的主元可能并不是最优的，此时在寻找主元时不能将方差作为衡量重要性的标准。

1.3 高维数据下的特征值分解

在现实生活中，由于样本的高维特征（图像识别每一个像素点都是一个维度），将会导致协方差矩阵$\Sigma$非常大，计算机难以存储与计算，那如何针对高维数据做特征值分解呢？我们知道$\Sigma =X{X}^T$，考虑替代矩阵$P=X^{T}X$，如果有100个样本，10000个维度，那么$\Sigma$就是一个10000维的方阵，而$P$只是一个100维的方阵。我们做如下推导：$$\begin{gathered} Pv=\lambda v \\ {X}^{T}Xv=\lambda v \\ X{X}^{T}Xv=X\lambda v \\ \Sigma (Xv)=\lambda (Xv) \end{gathered}$$
这样通过求$P$的特征值和特征向量，其特征值即$\Sigma$的特征值，其特征向量右乘随机向量$X$即为$\Sigma$的特征向量，从而完成高维数据的特征值分解（只能得到部分前几的主成分）。

求证：$A^{T}A$与$A{A}^T$的特征值相等？
补充：一般情况下，$AB$与$BA$的特征值只差$0$特征值 的个数。

若存在可逆矩阵$P$使得，$P^{-1}AP=B$，则矩阵$A$与矩阵$B$相似。若两矩阵相似，则它们的秩相等，行列式相等，迹（所有特征值的和，或主对角线元素的总和）相等，两者拥有同样的特征值（可能会多出一些为零的特征值）。

2. CCIPCA增量主元分析算法[1]

candid covariance-free incremental principle component analysis直接协方差无关的增量式主成分分析。

假设输入向量序列为$u^{\prime }(t),t=1,2,...$;第$n$幅图像输入时均值为$m(n)=\frac{1}{n}{\sum }_{t=1}^n{u}^{\prime }(t)$,它的协方差矩阵为$$A(n)=\frac{1}{n}\sum _{t=1}^n[u^{\prime }(t)-m(n)][u^{\prime }(t)-m(n){]}^T=\frac{1}{n}\sum _{i=1}^{n}u(t)u(t{)}^T\text{\hspace{1em}(1)}$$这里$u(t)=u^{\prime }(t)-m(n)$

$u(n)$的第$i$个特征值和特征值和特征向量的计算公式为${\lambda }_i{x}_i(n)=A(n)x_i(n)$,其中$x_i(n)$为第$n$输入时待求的第$i$个特征向量，${\lambda }_i$为对应的特征值 。CCIPCA算法为了加快迭代的速度，整个迭代是对特征值和特征向量的乘积${\lambda }_i{x}_i$进行的，设第$n$个输入时有$$v_i(n)={\lambda }_i{x}_i(n)=A(n)x_i(n)\text{\hspace{1em}(2)}$$
把（1）式代入（2）式，可得$$v_i(n)=\frac{1}{n}\sum _{t=1}^{n}u(t)u(t{)}^T{x}_i(n)\text{\hspace{1em}(3)}$$
若通过迭代获得特征值和特征向量的乘积$v_i$，因特征向量是归一的，只要对（2）式求模（内积，开根），可求得${\lambda }_i=||v_i||,x_i=\frac{v_i}{||v_i||}$。迭代采用（3）式，把$\frac{v_i(n-1)}{||v_i(n-1)||}$近似为$x_i(n)$代入（3）式，经变换可得CCIPCA的基本迭代公式，$$v_i(n)=\frac{n-1}{n}{v}_i(n-1)+\frac{1}{n}u(n)u(n{)}^T\frac{v_i(n-1)}{||v_i(n-1)||}\text{\hspace{1em}(4)}$$

其中$\frac{n-1}{n}$为上一步的迭代值 $v_i(n-1)$的权值，第2项的$\frac{1}{n}$相当于迭代的调整步长。$u(n)$作为第$n$幅新输入数据对迭代向量$v_i(n)$的调整，在迭代中$v_i(n)$逐步收敛到所求的第$i$个特征向量。对不同序号的特征向量，都可以用（4）式迭代，只是输入的向量$u(n)$不同。求最大的特征值对应的特征向量时，$u(n)$为机器人眼睛直接采到的第$n$个数据。在求第2，第3乃至更高维特征向量时，须作以下处理，已经通过迭代得到第1个特征向量，先设$u_1(n)=u(n)$，并把$u_1(n)$投影到上一个已经求到的特征向量上（现为第1个特征向量），求出残差图像$u_2(n)$，如下式表示，$$u_2(n)=u_1(n)-u_1^T(n)\frac{v_1(n)}{||v_1(n)||}\frac{v_1(n)}{||v_1(n)||}$$
$u_2(n)$作为求第2个特征向量的输入，类似的可求出第3，第4，…个特征向量。因残差图像和上1个特征向量所恢复的图像正交，从而可求出所有相互正交的特征向量。另外，每输入1幅新的数据时，均值也要更新，对（1）式，输入第$n$幅图像时的均值 采用如下迭代式，$$m(n)=\frac{n-1}{n}m(n-1)+\frac{1}{n}{u}^{\prime }(n)$$

3. KPCA

一篇关于KPCA的超级全的博客

• 原始PCA
  $X{X}^T{w}_i={\lambda }_i{w}_i$

• 利用函数$\Phi$将数据$X$映射到高维空间$\Phi (X)$
  $\Phi (X)\Phi (X{)}^T{w}_i={\lambda }_i{w}_i$

• 令$w_i={\sum }_{k=1}^N{\alpha }_i\Phi (x_i)=\Phi (X)\alpha$得
  $\Phi (X)\Phi (X{)}^T\Phi (X)\alpha ={\lambda }_i\Phi (X)\alpha$

• 两边左乘$\Phi (X{)}^T$得
  $\Phi (X{)}^T\Phi (X)\Phi (X{)}^T\Phi (X)\alpha ={\lambda }_i\Phi (X{)}^T\Phi (X)\alpha$

• 由核函数的性质可得$K=\Phi (X{)}^T\Phi (X)$
  $$\begin{gathered} KK\alpha ={\lambda }_{i}K\alpha \\ \Rightarrow K\alpha ={\lambda }_i\alpha \end{gathered}$$

此处注意，特征向量$w_i$应该进一步由$\Phi (X)\alpha$计算得出。但是，$\Phi (X)$是无法计算的，从而也就没办法求得特征向量。

此处，我们分析是否可以直接利用$\alpha$来对数据进行降维呢，答案是肯定的。

假如新来一个数据$x$，需要利用以上结果来得到降维的数据$\hat{x}$：

$$\begin{gathered} \hat{x}=w_i^{T}x \\ =(\Phi (x)\alpha {)}^T\Phi (x) \\ ={\alpha }^T\Phi (X{)}^T\Phi (x) \\ =[{\alpha }_1,...,{\alpha }_N][k(x_1,x),...,k(x_N,x){]}^T \end{gathered}$$

PCA与KPCA的区别

• PCA 仍然不失为一种好分析方法，数据呈非线性流形分布，或者说是各指标呈非线性关系时，对于线性分析分析方法来说可能效果不是特别好，但同时应注意的是它也是一种统计分析方法，实际经济指标中都存在线性相关性（信息冗余），这是符合统计规律的，完全不相关的经济数据是极其少见，也就是说要求的数据只要大致呈线性分布，而且有PCA计算简单，无需先验知识，无需参数设置等优点。
• PCA与线性核KPCA不完全一样，对于有P个指标的N个数据样本。PCA计算协方差阵为PXP维矩阵，它可以提取的主成分数为P，而KPCA是从核矩阵出发计算的，最大可以提取的主成分为N，为满足在特征空间中的样本均值为零，还要对核矩阵K进行特殊处理，这也是导致与线性核与原样本内积的不一致的原因。
• KPCA核函数与核参数难于选择。PCA的协方差矩阵的特征向量对应于各指标的主成分比重，从而能用原指标云解释主成分，而KPCA的是基于核矩阵的特征向量，与原指标没有对应关系，从而核主成分解释困难，其次KPCA将指标投影到高维特征空间后，而其实际的数据又是在原空间处理的，其数值在原空间中是否均有排序意义也值得进一步研究。

4. 2DPCA[2]

设$A_i,i=1,2,...,N$是$N$个样本图像，$A_i\in R^{m\times n}$。首先计算协方差矩阵，也就是总体散布矩阵$G_t$ $$G_t=\frac{1}{N}\sum _{i=1}^N(A_i-\hat{A}{)}^T(A_i-\hat{A})$$其中$\hat{A}=\frac{1}{N}{\sum }_{i=1}^N{A}_i$为全体样本的均值。

计算$G_t$的特征值和特征向量，取特征值累积项献率$\alpha =0.90.99$所对应的特征向量组成投影矩阵$U=[u_1,u_2,...,u_k]\in R^{n\times k}$ 。则，$F_i=A_{i}U\in R^{m\times k}$就是$A_j$的特征。可知原来二维的图像大小$m\times n$，现在降维为$m\times k$，$k$是根据$\alpha$来确定的。也就是说，实行特征提取后，只是压缩了图像矩阵列向量的位数，行向量位数不变。

5. 2D2DPCA

同于2DPCA只在列方向降了维数，降维的效果不理想。为了更好的降维，D.Q.Zhang和Z.H.Zhou提出了双向的二维主成分分析方法(2D2DPCA)，也就是在行和列两个方向都进行2DPCA处理。

对所有训练样本利用上述2DPCA处理之后得到新的训练样本$F_i,i=1,2,...,N$，其中$Y_i\in R^{m\times k}$。在新样本上构造协方差矩阵$G_t^{\ast }$ $$G_t^{\ast }=\frac{1}{N}\sum _{i=1}^N(F_i-\hat{F})(F_i-\hat{F}{)}^T$$

同理，求$G_t^{\ast }$的特征向值与特征向量，取特征值累计项献率为$\alpha$的特征向量，得到行方向的投影矩阵$V=[v_1,v_2,...,v_d]\in R^{m\times d}$。$F_i$的投影为$V^T{F}_i\in R^{d\times k}$。

至此，两个投影方向的最优投影矩阵$U$和$V$都求得了，对于图像$A_i\in R^{m\times n}$最终降维矩阵为$Y_i=V^T{A}_{i}U\in R^{d\times k}$，其重构图像为$A_i^{approx}=V{Y}_i{U}^T$

6. BDPCA

BDPCA分别对行和列方向进行数据降维。
将图像样本矩阵$X_i^{\prime }$分解成$p$个$1\times q$的行向量，则行方向总体散度矩阵为:

$$S_r=\frac{1}{np}\sum _{i=1}^n{\left[X_i^{\prime }-\overline{X}(n)\right]}^{\mathrm{T}}\left[X_i^{\prime }-\overline{X}(n)\right]$$

其中，$\overline{X}(n)$为图像样本矩阵的均值。取前$k_r$个最大特征值所对应的特征向量组成行方向投影矩阵为
$${\mathbf{W}}_r=\left[w_{r1},w_{r2},\cdots ,w_{r{k}_r}\right]$$

同理，列方向总体散度矩阵及前$k_c$个最大特征值所对应的特征向量组成列方向投影矩阵分别为
$$S_c=\frac{1}{nq}\sum _{i=1}^n\left[X_i^{\prime }-\overline{X}(n)\right]{\left[X_i^{\prime }-\overline{X}(n)\right]}^{\mathrm{T}}$$
$${\mathbf{W}}_c=\left[w_{c1},w_{c2},\cdots ,w_{c{k}_c}\right]$$
图像样本矩阵$X_i^{\prime }$所对应的特征矩阵为
$$Y=W_c^{\mathrm{T}}{X}_i^{\prime }{W}_c$$

BDPCA的特征矩阵维数仅为$k_c\times k_r$，因此其运算量要远小于CCIPCA。但BDPCA是2维的批处理计算。

参考文献

[1]Juyang Weng, Yilu Zhang, and Wey-Shiuan Hwang, “Candid covariance-free incremental principal component analysis,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 25, no. 8, pp. 1034–1040, Aug. 2003.
[2] Jian Y, David Z, Frangi A F, et al. Two-dimensional PCA: a new approach to appearance-based face representation and recognition.[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2004, 26(1):131-137.