【2022版】基于矩阵分解的PCA 白化&ZCA白化

【2022版】基于矩阵分解的PCA 白化&ZCA白化

@author: Heisenberg

主成分分析（principal component analysis, PCA）是一种常用的数据降维算法。主要思想是将$n$维的向量映射到$k$维上，这$k$维全新的正交特征也被成为主成分。通过计算数据矩阵$X$的协方差矩阵$\Sigma$，得到协方差矩阵的特征值特征向量，选择方差最大（也即特征值最大）的k个特征（坐标轴方向）组成的矩阵，以将原来的$X$变换到新的空间上。

做whitening白化（亦作sphering,球化）的处理的条件：

1. 使数据的不同维度去相关；
2. 使数据每个维度的方差为1；

1、协方差矩阵及散度矩阵

在数据处理之前，先对数据矩阵$X$做中心化，则样本方差为

$S^2=\frac{1}{n-1}{\Sigma }_{i=1}^N(x_i-\bar{x}{)}^2$，用矩阵描述即$S^2=\frac{1}{n-1}\mathbf{X}{\mathbf{X}}^{\mathbf{T}}$

而散度矩阵为

$S^2={\Sigma }_{i=1}^n({\mathbf{x}}_{\mathbf{i}}-\bar{\mathbf{x}})({\mathbf{x}}_{\mathbf{i}}-\bar{\mathbf{x}}{)}^T$，即$S^2=\mathbf{X}{\mathbf{X}}^{\mathbf{T}}$

可以看到散度矩阵即协方差矩阵乘以（样本量-1）。二者的特征值和特征向量一样的。

基本思路即对$S^2$进行分解。

补充

如果是对列进行降维，则是对$S^2={\mathbf{X}}^T\mathbf{X}$进行矩阵分解。

2、特征值矩阵分解

2.1 特征值与特征向量

若一个向量$v$是矩阵$\mathbf{A}$的特征向量，则可以表示为$Av=\lambda v$，其中$\lambda$称之为特征值，$v$是正交向量。

2.2 特征值分解矩阵

对于矩阵$A$，有一组特征向量$v$，将这组向量正交化单位化后，能得到一组正交单位向量。

$\mathbf{A}=Q\Sigma Q^{-1}$.

其中$Q$是A的特征向量组成的矩阵，$\Sigma$是一个对角阵，对角线上的元素就是特征值。实对称矩阵都可以分解为实特征向量和实特征值。且只有方阵才可以进行特征值分解。

3、SVD矩阵分解

对于任意矩阵$A$总是存在如下的奇异值分解：

$A=U\Sigma V^T$

也即$A_{m\times n}=U_{m\times m}{\Sigma }_{m\times n}{V}_{n\times n}^T$, 展开形式如下所示：
$$\left(\begin{array}{cccc}{x}_{11}& x_{12}& & x_{1n}\\ & & \ddots & \\ x_{m1}& & & x_{mn}\end{array}\right)=\left(\begin{array}{ccc}{u}_{11}& & u_{m1}\\ & \ddots & \\ u_{1m}& & u_{mm}\end{array}\right)\left(\begin{array}{ccccc}{\sigma }_1& & & & 0\\ & \ddots & & & \\ & & {\sigma }_r& & \\ & & & \ddots & \\ 0& & & & 0\end{array}\right)\left(\begin{array}{ccc}{v}_{11}& & v_{1n}\\ & \ddots & \\ v_{n1}& & v_{nn}\end{array}\right)$$
其中，

$U$是左奇异向量，是$A{A}^T$单位化后的特征向量构成的正交矩阵，有$U^{T}U=I$,$U^T=U^{-1}$；

$V$是右奇异向量，是$A^{T}A$单位化后的特征向量构成的正交矩阵，有$V^{T}V=I$,$V^T=V^{-1}$，所以奇异值分解也可以写为$AV=U\Sigma$；

$\Sigma$是$A^{T}A$(或$A{A}^T$)的特征值求平方根后的对角矩阵。

Remarks：

所以SVD分解能够同时按行按列进行降维：
$$X_{k\times r}^{\prime }=U_{k\times m}^T{X}_{m\times n}{V}_{n\times r}$$

图解：

图源自维基百科，可见左右奇异矩阵$U,V$都是旋转坐标轴，对角矩阵$\Sigma$是为了压缩数据矩阵。

4、基于特征值分解的PCA降维

4.1 基本步骤

对于数据矩阵$X=[x_1,x_2,\cdots ,x_n]$。其中$x_i$都是列向量

1. 去中心化。i.e. $x_i-\bar{x}$
2. 就算协方差矩阵$S^2=\frac{1}{n}\mathbf{X}{\mathbf{X}}^{\mathbf{T}}$（按行），$S^2=\frac{1}{n}{\mathbf{X}}^{\mathbf{T}}\mathbf{X}$（按列）
3. 通过特征值分解法求$S^2$的特征值与特征向量。
4. 对特征值按照从大到小的顺序，选取最大的$K$个。
   • 如果要按行降维，将对应的K个特征向量分别作为行向量组成特征向量矩阵$P$.
   • 如果要按列降维，将对应的K个特征向量分别作为列向量组成特征向量矩阵$P
5. 将数据转换到K个特征向量构建的新空间中，i.e. $Y_{k\times n}=P_{k\times m}{X}_{m\times n}$（按行）或$Y_{m\times k}=X_{m\times n}{P}_{n\times k}$

4.2 举例

举个栗子更形象些：

For 矩阵$X=\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)$

1. 对列做降维，类似于将特征整合。

   • 求其协方差矩阵
     $$S^2=\frac{1}{5}\left(\begin{array}{cc}-1& -2\\ -1& 0\\ 0& 0\\ 2& 1\\ 0& 1\end{array}\right)\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)=\left(\begin{array}{ccccc}1& -\frac{1}{5}& 0& -\frac{4}{5}& -\frac{2}{5}\\ -\frac{1}{5}& \frac{1}{5}& 0& \frac{2}{5}& 0\\ 0& 0& 0& 0& 0\\ -\frac{4}{5}& \frac{2}{5}& 0& 1& \frac{1}{5}\\ -\frac{2}{5}& 0& 0& \frac{1}{5}& \frac{1}{5}\end{array}\right)$$

   • 求解得到特征值为${\lambda }_1=2,{\lambda }_2=\frac{2}{5}$，${\lambda }_3={\lambda }_4={\lambda }_5=0$，对应的标准化后的特征向量为
     $$q_1=\left(\begin{array}{c}-\frac{3}{\sqrt{20}}\\ \frac{1}{\sqrt{20}}\\ 0\\ \frac{3}{\sqrt{20}}\\ \frac{1}{\sqrt{20}}\end{array}\right),q_2=\left(\begin{array}{c}-\frac{1}{2}\\ -\frac{1}{2}\\ 0\\ -\frac{1}{2}\\ \frac{1}{2}\end{array}\right),q_3=q_4=q_5=0\in {\mathbf{R}}^{5\times 1}$$

   • 选取K=2，则线性变换矩阵如下：
     $$\mathbf{P}=\left(\begin{array}{cc}-\frac{3}{\sqrt{20}}& -\frac{1}{2}\\ \frac{1}{\sqrt{20}}& -\frac{1}{2}\\ 0& 0\\ \frac{3}{\sqrt{20}}& -\frac{1}{2}\\ \frac{1}{\sqrt{20}}& \frac{1}{2}\end{array}\right)\in {\mathbf{R}}^{5\times 2}$$

   • 则可得到$X^{\text{PCA}}=X_{2\times 5}{P}_{5\times 2}=\left(\begin{array}{cc}\sqrt{5}& -1\\ \sqrt{5}& 1\end{array}\right)$

2. 对行做降维，类似于将样本分类、总结。

   • 求其协方差矩阵
     $$S^2=\frac{1}{5}\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)\left(\begin{array}{cc}-1& -2\\ -1& 0\\ 0& 0\\ 2& 1\\ 0& 1\end{array}\right)=\left(\begin{array}{cc}\frac{6}{5}& \frac{4}{5}\\ \frac{4}{5}& \frac{6}{5}\end{array}\right)$$

   • 求解得到特征值为${\lambda }_1=2,{\lambda }_2=\frac{2}{5}$，对应的特征向量为$q_1=[1,1{]}^T,q_2=[-1,1{]}^T$

   • 将特征向量标准化得到$q_1=[\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}{]}^T,q_2=[-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}{]}^T$

   • 则线性变换矩阵$P=\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\\ -\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right)$

   • 选取K=1，则可得到

$$\begin{gathered} Y_{1\times 5}=P_{1\times 2}{X}_{2\times 5} \\ =\left(\begin{array}{ll}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right)\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right) \\ =\left(\begin{array}{lllll}-\frac{3}{\sqrt{2}}& -\frac{1}{\sqrt{2}}& 0& \frac{3}{\sqrt{2}}& -\frac{1}{\sqrt{2}}\end{array}\right) \end{gathered}$$

5、基于SVD分解的PCA降维

5.1 基本步骤

对于数据矩阵$X=[x_1,x_2,\cdots ,x_n]$。其中$x_i$都是列向量

1. 去中心化。i.e. $x_i-\bar{x}$

2. 就算协方差矩阵$S_{\text{row}}^2=\frac{1}{n-1}\mathbf{X}{\mathbf{X}}^{\mathbf{T}}$和$S_{\text{col}}^2=\frac{1}{n-1}{\mathbf{X}}^{\mathbf{T}}\mathbf{X}$（SVD分解能同时按行、按列降维。）

3. 通过SVD分解法求两个$S^2$的特征值与特征向量。

4. 对特征值按照从大到小的顺序，选取最大的$K$个。

   将对应的K个特征向量中，U分别作为行向量，V分别作为列向量组成特征向量矩阵$P$

5. 将数据转换到K个特征向量构建的新空间中。

5.2 优点

• 某些SVD算法不用先求$X{X}^T$就求出右奇异矩阵，降低运算开销。
• 左右两个奇异矩阵可以在行和列两个方向对数据进行压缩。i.e. $Y_{m\times k}=X_{m\times n}{V}_{n\ast k}^T$ ，右奇异矩阵对列进行压缩; $Y_{k\times n}=U_{k\times m}^T{X}_{m\times n}$，左奇异矩阵对行数进行特征压缩。(在此说的维度转置与否都是要与原矩阵对应的)

5.3 举例

举个栗子：

对于矩阵$A=\left(\begin{array}{ccc}3& 2& 2\\ 2& 3& -2\end{array}\right)$，可得到：

$$A{A}^T=\left(\begin{array}{cc}17& 8\\ 8& 17\end{array}\right),A^{T}A=\left(\begin{array}{ccc}13& 12& 2\\ 12& 13& -2\\ 2& -2& 8\end{array}\right)$$

1. 对$A{A}^T$求特征值为${\lambda }_1=25,{\lambda }_2=9$; 特征向量为$u_1=\left(\begin{array}{c}\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\end{array}\right)u_2=\left(\begin{array}{c}\frac{1}{\sqrt{2}},-\frac{1}{\sqrt{2}}\end{array}\right)$
2. 对$A^{T}A$求特征值为${\lambda }_1=25,{\lambda }_2=9,{\lambda }_3=0$; 特征向量为$v_1=\left(\begin{array}{c}1/\sqrt{2}\\ 1/\sqrt{2}\\ 0\end{array}\right)v_2=\left(\begin{array}{c}1/\sqrt{18}\\ -1/\sqrt{18}\\ 4/\sqrt{18}\end{array}\right)v_3=\left(\begin{array}{c}2/3\\ -2/3\\ -1/3\end{array}\right)$
3. 则A的SVD分解可以写作如下：

$$A=US{V}^T=\left(\begin{array}{cc}1/\sqrt{2}& 1/\sqrt{2}\\ 1/\sqrt{2}& -1/\sqrt{2}\end{array}\right)\left(\begin{array}{ccc}5& 0& 0\\ 0& 3& 0\end{array}\right)\left(\begin{array}{rrr}1/\sqrt{2}& 1/\sqrt{2}& 0\\ 1/\sqrt{18}& -1/\sqrt{18}& 4/\sqrt{18}\\ 2/3& -2/3& -1/3\end{array}\right)$$

• 选取K=2，则对A进行列的压缩可以得到：

$Y_{(2\times 2)}=A_{(2\times 3)}{V}_{(3\ast 2)}$

也即
$$Y=\left(\begin{array}{ccc}3& 2& 2\\ 2& 3& -2\end{array}\right)\left(\begin{array}{cc}1/\sqrt{2}& 1/\sqrt{18}\\ 1/\sqrt{2}& -1/\sqrt{18}\\ 0& 4/\sqrt{18}\end{array}\right)=\left(\begin{array}{cc}5/\sqrt{2}& 3/\sqrt{2}\\ 5/\sqrt{2}& -3/\sqrt{2}\end{array}\right)$$

• 选取K=1，对A进行行的压缩可以得到：

$Y_{(1\times 3)}=U_{(1\times 2)}^T{A}_{(2\times 3)}$
$$Y=\left(\begin{array}{cc}1/\sqrt{2}& 1/\sqrt{2}\end{array}\right)\left(\begin{array}{ccc}3& 2& 2\\ 2& 3& -2\end{array}\right)=\left(\begin{array}{ccc}5/\sqrt{2}& 5/\sqrt{2}& 0\end{array}\right)$$
其中$u_1=\left(\begin{array}{c}1/\sqrt{2}\\ 1/\sqrt{2}\end{array}\right)$称之为第一主成分，以此类推。

5.4 方差贡献率

对于特征值对角矩阵$\Sigma ={\left(\begin{array}{ccc}{\sigma }_1& & \\ & \ddots & \\ & & {\sigma }_r\end{array}\right)}_{(r\times r)}$

如果我们想保留99%的方差贡献率则使得${\alpha }_i=\frac{{\Sigma }_{i=1}^k{\sigma }_i}{{\Sigma }_{j=1}^r{\sigma }_j}\ge 0.99$即可.

Generally，75%以上就合格。

Remarks

原则上来说$A$矩阵的特征值分解的$Q$矩阵和A矩阵的$V^T$矩阵是相等的。但是由于某些包算法上的区别使得列的排列顺序、正负号、数值大小有细微的区别，但仍是相等的。（排列顺序、正负号不改变矩阵相等是因为Q,V都是orthography matrix，正交矩阵保证该性质。数值差异大概率是因为下面PArt6中正则化导致的。）

如该例中

In [2]: a = np.array([[3,2,2],[2,3,-2]])

In [3]: a
Out[3]:
array([[ 3,  2,  2],
       [ 2,  3, -2]])

In [5]: eig_vals,eig_vecs=np.linalg.eig(a.T@a)

In [6]: eig_vals
Out[6]: array([2.50000000e+01, 3.44694725e-15, 9.00000000e+00])

In [11]: eig_vecs
Out[11]:
array([[-7.07106781e-01, -6.66666667e-01,  2.35702260e-01],
       [-7.07106781e-01,  6.66666667e-01, -2.35702260e-01],
       [-1.16847159e-17,  3.33333333e-01,  9.42809042e-01]])

In [12]: u,s,vh = np.linalg.svd(a)

In [23]: vh.T
Out[23]:
array([[-7.07106781e-01, -2.35702260e-01, -6.66666667e-01],
       [-7.07106781e-01,  2.35702260e-01,  6.66666667e-01],
       [-6.47932334e-17, -9.42809042e-01,  3.33333333e-01]])

引自PCA SVD原理详解及应用

最后说一下两者的联系，其实svd可以用来做pca的，因为在pca阶段我们主要求得是$X{X}^T$特征值和特征向量，对吧（暂不考虑减去均值），而在求svd的左奇异矩阵（U）过程中，也是求得是$X{X}^T$特征值和特征向量对吧，也就是说这里的SVD奇异值的平方就是pca中的特征值对吧，其实在SVD的分解过程中，不用像我们上面讨论的那样先求 ， 特征值和特征向量方法进行分解，而是有很多快速算法，但我们可以根据SVD的结果进行PCA对吧，其实scikit-learn内部的PCA也是用SVD做的

6、正则化

若在举例5.3的步骤2中不对特征向量做标准化

即$X_{PCA-White}={\Sigma }^{-\frac{1}{2}}{U}^{T}X=\left[\begin{array}{ccc}\frac{1}{\sqrt{{\lambda }_1}}& 0& 0\\ 0& \ddots & 0\\ 0& 0& \frac{1}{\sqrt{{\lambda }_n}}\end{array}\right]X_{rotate}$

Proof:
$$\begin{array}{rl}{S}_{PCA\text{ hite }}^2& =\frac{1}{m}{X}_{\text{PCAwhite }}{X}_{\text{PCAwhite }}^T\\ & ={\Sigma }^{-\frac{1}{2}}{U}^T\left(\frac{1}{m}X{X}^T\right)U{\left({\Sigma }^{-\frac{1}{2}}\right)}^T\\ & ={\Sigma }^{-\frac{1}{2}}{U}^T{S}^{2}U{\Sigma }^{-\frac{1}{2}}\because \text{对角阵}\\ & ={\Sigma }^{-\frac{1}{2}}\left(U^{T}U\right)\Sigma \left(U^{T}U\right){\Sigma }^{-\frac{1}{2}}\because \text{SVD分解}\\ & ={\Sigma }^{-\frac{1}{2}}I\Sigma I{\Sigma }^{-\frac{1}{2}}\because \text{酉矩阵}U{U}^T=U^{T}U=I\\ & =I\end{array}$$
通常在实践中，一般选$X_{PCA-White,i}=\frac{X_{rot,i}}{\sqrt{{\lambda }_i+ϵ}}$，$ϵ=10^{-5}$.

原因：

1. 有时一些特征值在数值上接近0，在缩放步骤时将导致除以一个接近0的值；这可能使数据上溢或造成数值不稳定；
2. 对图像来说，正则化项对输入图像也有一些平滑去噪（或低通滤波）的作用，可改善学习到的特征。

7、ZAC白化

Zero-phase Component Analysis Whitening 零相位成分分析白化

PCA and ZCA whitening differ only by a rotation.

——《Neural Networks: Tricks of the Trade》

ZCA白化只是在PCA白化的基础上做了一个旋转操作，使得白化之后的数据映射回源空间更加的接近原始数据。

$Y_{ZAC,(m\times n)=U_{(m\times k)}{Y}_{PCA,(k\times n)}}=U_{(m\times k)}{U}_{(k\times m)}^T{X}_{m\times n}$

下图可以更有助加深理解。

图片源自StackChange.

图1为原始数据，图2为PCA选取的最大方差方向的新坐标轴，下方对应的为映射到新坐标轴后，原始阴影部分数据在新坐标空间的最东边，而图3进行ZAC坐标旋转之后，将阴影部分数据转到了新坐标空间的东北部分。

Proof：
$$\begin{array}{rl}\sum _{\text{ZCAwhite }}& =\frac{1}{m}{X}_{\text{ZCAwhite }}{X}_{\text{ZCAwhite }}^T\\ & =U\frac{1}{m}{X}_{\text{PCAwhile }}{X}_{\text{PCAwhite }}{U}^T\\ & =UI{U}^T\\ & =I\end{array}$$
可见其仍能保证白化的合法化，协方差矩阵是对角阵的单位矩阵。

举例1（特征值分解）：

$Y_{ZCA,(2\times 5)}=U_{(2\times 1)}{Y}_{PCA,(1\times 5)}$
$$Y_{ZAC}=\left(\begin{array}{l}\frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\end{array}\right)\left(\begin{array}{lllll}-\frac{3}{\sqrt{2}}& -\frac{1}{\sqrt{2}}& 0& \frac{3}{\sqrt{2}}& -\frac{1}{\sqrt{2}}\end{array}\right)=\left(\begin{array}{ccccc}-3/2& -1/2& 0& 3/2& -1/2\\ -3/2& -1/2& 0& 3/2& -1/2\end{array}\right)$$
举例2（SVD分解）：

$Y_{ZCA,(2\times 3)}=U_{(2\times 1)}{Y}_{PCA,(1\times 3)}$
$$Y_{ZCA}=\left(\begin{array}{c}1/\sqrt{2}\\ 1/\sqrt{2}\end{array}\right)\left(\begin{array}{ccc}5/\sqrt{2}& 5/\sqrt{2}& 0\end{array}\right)=\left(\begin{array}{ccc}5/2& 5/2& 0\\ 5/2& 5/2& 0\end{array}\right)$$

8、应用

LSI应用

在机器学习中的数学(5)-强大的矩阵奇异值分解(SVD)及其应用中引用《数学之美》中吴军老师的话：

“三个矩阵有非常清楚的物理含义。第一个矩阵X中的每一行表示意思相关的一类词，其中的每个非零元素表示这类词中每个词的重要性（或者说相关性），数值越大越相关。最后一个矩阵Y中的每一列表示同一主题一类文章，其中每个元素表示这类文章中每篇文章的相关性。中间的矩阵则表示类词和文章雷之间的相关性。因此，我们只要对关联矩阵A进行一次奇异值分解，w 我们就可以同时完成了近义词分类和文章的分类。（同时得到每类文章和每类词的相关性）。”

参考资料

Machine Learning — Singular Value Decomposition (SVD) & Principal Component Analysis (PCA)

主成分分析（PCA）原理详解

斯坦佛UFLDL PCA Whitening

What is the difference between ZCA whitening and PCA whitening?

LDA、PCA、ZCA、ICA回顾

PCA SVD原理详解及应用

机器学习中的数学(5)-强大的矩阵奇异值分解(SVD)及其应用