1500字范文 > 【机器学习】深入剖析主成分分析（PCA）与协方差矩阵

【机器学习】深入剖析主成分分析（PCA）与协方差矩阵

时间：2022-07-20 04:25:37

其他机器学习系列文章见于专题：机器学习进阶之路——学习笔记整理，欢迎大家关注。

1. 线性降维方法

在高维情形下出现的数据样本稀疏、距离计算困难问题，是所有机器学习方法共同面临的严重障碍，被称为“维数灾难”（curse of dimensionality）。

缓解维数灾难的一个重要途径是降维（dimension reduction），即通过某种数学变换将原始高维属性空间转变为一个低维“子空间”，在这个子空间中样本密度大幅提高，距离计算也变得更为容易。

对原始高维空间进行线性变换，来进行降维的方法称为线性降维方法。它们都符合以下基本形式：

Z=WTX\mathbf { Z } = \mathbf { W } ^ { \mathrm { T } } \mathbf { X } Z=WTX

其中，样本矩阵为X=(x1,x2,…,xm)∈Rd×m\mathbf {X}= \left( \boldsymbol { x } _ { 1 } , \boldsymbol { x } _ { 2 } , \ldots , \boldsymbol { x } _ { m } \right) \in \mathbb { R } ^ { d \times m }X=(x1,x2,…,xm)∈Rd×m，投影后样本矩阵为Z=WTXZ = W ^ { T } XZ=WTX，W∈Rd×d′\mathbf {W} \in \mathbb { R } ^ { d \times d ^ { \prime } }W∈Rd×d′是投影矩阵，将维度从ddd降到d′d'd′，Z∈Rd′×m\mathbf { Z } \in \mathbb { R } ^ { d ^ { \prime } \times m }Z∈Rd′×m是样本在新坐标系中的表达。

不同线性降维方法的区别在于对低维子空间的性质有不同的要求，相当于对投影矩阵WWW施加不同的约束。我们下面将详细介绍的PCA就是在WTW=IW^TW=IWTW=I的约束下的线性降维方法。

2. PCA概念

主成分分析（Principal Component Analysis，PCA）是最常用的一种降维方法，通过一个投影矩阵将可能存在相关性和冗余的特征转换为一组更低维度的线性不相关的特征，转换后的特征就叫做主成分。

3. PCA原理

在降维的过程中，我们希望损失的信息尽可能少，也就是希望保留的信息尽可能多。PCA用方差来度量信息量，在某个维度上，样本分布越分散，方差越大，信息越多。因此，PCA对投影矩阵的第一个要求是使投影后的样本在各维度上方差尽可能大。

然而，如果单纯只选择方差最大的方向，会导致选择的基向量方向差不多，彼此相关性大，表示的信息几乎是重复的。所以为了使降维后的维度能尽可能地表达信息，第二个要求是不希望投影后的特征之间存在（线性）相关性。

综上，PCA的优化目标是：

降维后在新维度上的同一维度方差最大；不同维度之间相关性为0。

根据线性代数可知，投影值的协方差矩阵的对角线代表了投影值在各个维度上的方差，其他元素代表各个维度之间的相关性（协方差）。基于优化目标，我们希望协方差矩阵是个对角线上值很大的对角矩阵（对角矩阵意味着非对角元素为0）。

相同维度的多维随机变量X=[X1,X2,…,Xn]T\mathbf { X } = \left[ X _ { 1 } , X _ { 2 } , \dots , X _ { n } \right] ^ { T }X=[X1,X2,…,Xn]T和Y=[Y1,Y2,…,Yn]T\mathbf { Y } = \left[ Y _ { 1 } , Y _ { 2 } , \dots , Y _ { n } \right] ^ { T }Y=[Y1,Y2,…,Yn]T的协方差矩阵的第(i,j)(i,j)(i,j)项（第(i,j)(i,j)(i,j)项是一个协方差）的定义为：

cov⁡(Xi,Yj)=E[(Xi−E[Xi])(Yj−E[Yj])]=E[XiYj−XiE[Yj]−E[Xi]Yj+E[Xi]E[Yj]]=E[XiYj]−E[Xi]E[Yj]−E[Xi]E[Yj]+E[Xi]E[Yj]=E[XiYj]−E[Xi]E[Yj]\begin{aligned} \operatorname { cov } ( X_i , Y_j ) & = \mathrm { E } [ ( X_i - \mathrm { E } [ X_i ] ) ( Y_j - \mathrm { E } [ Y_j ] ) ] \\ & = \mathrm { E } [ X_i Y_j - X _i\mathrm { E } [ Y_j ] - \mathrm { E } [ X_i ] Y_j + \mathrm { E } [ X_i ] \mathrm { E } [ Y_j ] ] \\ & = \mathrm { E } [ X_i Y_j ] - \mathrm { E } [ X_i ] \mathrm { E } [ Y_j ] - \mathrm { E } [ X_i ] \mathrm { E } [ Y_j ] + \mathrm { E } [ X_i ] \mathrm { E } [ Y_j ] \\ & = \mathrm { E } [ X_i Y_j ] - \mathrm { E } [ X_i ] \mathrm { E } [ Y_j ] \end{aligned} cov(Xi,Yj)=E[(Xi−E[Xi])(Yj−E[Yj])]=E[XiYj−XiE[Yj]−E[Xi]Yj+E[Xi]E[Yj]]=E[XiYj]−E[Xi]E[Yj]−E[Xi]E[Yj]+E[Xi]E[Yj]=E[XiYj]−E[Xi]E[Yj]

如果能对X\mathrm { X}X和Y\mathrm {Y}Y分别进行中心化，那么E[X]=0\mathrm { E } [ X ]=0E[X]=0，E[Y]=0\mathrm { E } [ Y ]=0E[Y]=0，将大大简化协方差的计算。这解释了为什么要对数据样本进行中心化。

假设数据样本已经进行了中心化，得到样本矩阵为X=(x1,x2,…,xm)∈Rd×m\mathbf {X}= \left( \boldsymbol { x } _ { 1 } , \boldsymbol { x } _ { 2 } , \ldots , \boldsymbol { x } _ { m } \right) \in \mathbb { R } ^ { d \times m }X=(x1,x2,…,xm)∈Rd×m，投影后样本矩阵为Z=WTXZ = W ^ { T } XZ=WTX，W∈Rd×d′\mathbf {W} \in \mathbb { R } ^ { d \times d ^ { \prime } }W∈Rd×d′是投影矩阵，将维度从ddd降到d′d'd′，Z∈Rd′×m\mathbf { Z } \in \mathbb { R } ^ { d ^ { \prime } \times m }Z∈Rd′×m是样本在新坐标系中的表达。投影后的协方差矩阵为：

C=1m−1∑i=1mziziT=1m−1ZZT=1m1(WTX)(WTX)T=1m−1WTXXTW\begin{array} { l } { C = \frac { 1 } { m - 1 } \sum _ { i = 1 } ^ { m } z _ { i } z _ { i } ^ { T } } \\ { = \frac { 1 } { m - 1 } Z Z ^ { T } } \\ { = \frac { 1 } { m _ { 1 } } \left( W ^ { T } X \right) \left( W ^ { T } X \right) ^ { T } } \\ { = \frac { 1 } { m - 1 } W ^ { T } X X ^ { T } W } \end{array} C=m−11∑i=1mziziT=m−11ZZT=m11(WTX)(WTX)T=m−11WTXXTW

在这里必须要声明的是，如果我们用mmm代表样本数，ddd代表维度数。
1）如果XXX的维度如果是m×dm \times dm×d的话，WWW的维度依然是d×d′d \times d'd×d′，投影后样本矩阵为Z=XWZ = XWZ=XW，那么PCA要计算的协方差矩阵是XTXX^TXXTX。
2）如果XXX的维度如果是d×md \times md×m的话，WWW的维度依然是d×d′d \times d'd×d′，投影后样本矩阵为Z=WTXZ = W ^ { T } XZ=WTX，那么PCA要计算的协方差矩阵是XXTXX^TXXT。
本文按照《机器学习（西瓜书）》里面的写法，采用第二种表达方式，如果要进行代码实现，这种表达方式下的XXX是我们通常理解的样本矩阵的转置。

要使投影后的协方差矩阵为对角矩阵，就要找到能使投影前的协方差矩阵XXTXX ^ { T }XXT对角化的矩阵WWW（由于协方差矩阵XXTXX ^ { T }XXT是一个实对称矩阵，那么必然存在一个可逆矩阵可使其对角化，且相似对角阵上的元素即为矩阵本身特征值）。然后对投影前协方差矩阵XXTXX ^ { T }XXT进行特征值分解，将求得的特征值排序：λ1≥λ2≥…≥λd\lambda _ { 1 } \geq \lambda _ { 2 } \geq \ldots \geq \lambda _ { d }λ1≥λ2≥…≥λd，取前d′d'd′个特征值对应的特征向量构成W=(w1,w2,…,wd′)W = \left( w _ { 1 } , w _ { 2 } , \dots , w _ { d ^ { \prime } } \right)W=(w1,w2,…,wd′)，这就是主成分分析的解。

严格来说，协方差矩阵是1m−1XXT\frac{1}{m-1}X{X^T}m−11XXT，但前面常数项不影响，为方便描述我们指的的协方差矩阵是XXTX{X^T}XXT。

4. 严格推导过程

上一节中提到，PCA的优化目标为：

目标1：降维后同一维度方差最大；

目标2：不同维度之间相关性为0。

我们可以将PCA的优化目标转化为在不同维度之间相关性为0的约束条件下，求解使同一维度方差最大化的投影矩阵的问题。

从两个角度，可以等价地解释方差最大化：

1. 最近重构性：样本点到投影的超平面（直线的高维推广）的距离都足够近；

2. 最大可分性：样本点在这个超平面上的投影都尽可能分开。

根据最近重构性和最大可分性，能够得到主成分分析的两种等价推导。

（一）基于最近重构性推导PCA

假定数据样本已经进行了中心化，即∑i=1mxi=0\sum _ { i = 1 } ^ { m } x _ { i } = 0∑i=1mxi=0，样本点xi=(xi1,xi2,…,xid)T{x_i} = {\left( {{x_{i1}},{x_{i2}}, \ldots ,{x_{id}}} \right)^T}xi=(xi1,xi2,…,xid)T是d×1d \times 1d×1维的列向量，X=(x1,…,xi,…,xm)X = \left( {{x_1}, \ldots ,{x_i}, \ldots ,{x_m}} \right)X=(x1,…,xi,…,xm)由mmm个ddd维列向量xix_ixi构成，XXX的维度为d×md \times md×m。如果丢弃部分坐标维度，将维度从ddd降低到d′d'd′，则标准正交基构成的投影矩阵W=(w1,…,wi,…,wd′)W = \left( {{w_1}, \ldots ,{w_i}, \ldots ,{w_{d'}}} \right)W=(w1,…,wi,…,wd′)的维度为d×d′d \times d'd×d′，由d′d'd′个ddd维列向量wiw_iwi构成，wi=(wi1,wi2,…,wid)T{w_i} = {\left( {{w_{i1}},{w_{i2}}, \ldots ,{w_{id}}} \right)^T}wi=(wi1,wi2,…,wid)T是标准正交基向量，∥wi∥2=1,wiTwj=0(i≠j)\left\| w _ { i } \right\| _ { 2 } = 1, w _ { i } ^ { T } w _ { j } = 0 ( i \neq j )∥wi∥2=1,wiTwj=0(i̸=j)。投影变换后的新坐标系为(w1,…,wi,…,wd′)\left( {{w_1}, \ldots ,{w_i}, \ldots ,{w_{d'}}} \right)(w1,…,wi,…,wd′)，样本点xix_ixi在低维坐标系中的投影为zi=(zi1,zi2,…,zid′)z _ { i } = \left( z _ { i 1 } , z _ { i 2 } , \dots , z _ { i d ^ { \prime } } \right)zi=(zi1,zi2,…,zid′)，zij=wjTxiz _ { i j } = w _ { j } ^ { T } x _ { i }zij=wjTxi是xix_ixi在低维坐标系下第jjj维的坐标。若基于zjz_jzj来重构xix_ixi，则得到x^i=∑j=1d′zijwj=Wzi\hat { x } _ { i } = \sum _ { j= 1 } ^ { d' } z _ { i j } w _ { j } = W z _ { i }x^i=∑j=1d′zijwj=Wzi。

考虑整个训练集，我们想要使原样本点xix_ixi与基于投影重构的样本点x^i\hat { x } _ { i }x^i之间的距离为：

(1)∑i=1m∥x^i−xi∥22=∑i=1m∥∑j=1d′zijwj−xi∥22=∑i=1m(Wzi−xi)2=∑i=1m[(Wzi)T(Wzi)−2(Wzi)Txi+xiTxi]=∑i=1m(ziTWTWzi−2ziTWTxi+xiTxi)=WTW=I,zi=WTxi∑i=1m(ziTzi−2ziTzi+xiTxi)=∑i=1m(−ziTzi+xiTxi)=zi=WTxi−∑i=1m(WTxi)T(WTxi)+∑i=1mxiTxi=−∑i=1mxiTWWTxi+∑i=1mxiTxi\begin{array}{l} \sum\limits_{i = 1}^m {\left\| {{{\widehat x}_i} - {x_i}} \right\|_2^2} = \sum\limits_{i = 1}^m {\left\| {\sum\limits_{j = 1}^{d'} {{z_{ij}}{w_j} - {x_i}} } \right\|} _2^2\\ = \sum\limits_{i = 1}^m {{{\left( {W{z_i} - {x_i}} \right)}^2}} \\ = \sum\limits_{i = 1}^m {\left[ {{{\left( {W{z_i}} \right)}^T}\left( {W{z_i}} \right) - 2{{\left( {W{z_i}} \right)}^T}{x_i} + x_i^T{x_i}} \right]} \\ = \sum\limits_{i = 1}^m {\left( {z_i^T{W^T}W{z_i} - 2z_i^T{W^T}{x_i} + x_i^T{x_i}} \right)} \\ \mathop = \limits^{{W^T}W = I,{z_i} = {W^T}{x_i}} \sum\limits_{i = 1}^m {\left( {z_i^T{z_i} - 2z_i^T{z_i} + x_i^T{x_i}} \right)} \\ = \sum\limits_{i = 1}^m {\left( { - z_i^T{z_i} + x_i^T{x_i}} \right)} \\ \mathop = \limits^{{z_i} = {W^T}{x_i}} - \sum\limits_{i = 1}^m {{{\left( {{W^T}{x_i}} \right)}^T}\left( {{W^T}{x_i}} \right)} + \sum\limits_{i = 1}^m {x_i^T{x_i}} \\ {\rm{ = }} - \sum\limits_{i = 1}^m {x_i^TW{W^T}x_i} + \sum\limits_{i = 1}^m {x_i^T{x_i}} \end{array}\tag {1} i=1∑m∥xi−xi∥22=i=1∑m∥∥∥∥∥j=1∑d′zijwj−xi∥∥∥∥∥22=i=1∑m(Wzi−xi)2=i=1∑m[(Wzi)T(Wzi)−2(Wzi)Txi+xiTxi]=i=1∑m(ziTWTWzi−2ziTWTxi+xiTxi)=WTW=I,zi=WTxii=1∑m(ziTzi−2ziTzi+xiTxi)=i=1∑m(−ziTzi+xiTxi)=zi=WTxi−i=1∑m(WTxi)T(WTxi)+i=1∑mxiTxi=−i=1∑mxiTWWTxi+i=1∑mxiTxi(1)

有人可能会有疑问，该式中的WWTWW^TWWT不是等于III吗，那么式（1）不就等于0了吗？注意，WWW是d×d′d \times d'd×d′维度的矩阵，不是正交矩阵，尽管WTW=I{W^T}W = IWTW=I，但WWT≠IW{W^T} \ne IWWT̸=I，不能和式（1）的第二项抵消。

由于xiTWWTxi{x_i^TW{W^T}x_i}xiTWWTxi是一个标量，不是一个向量，可以用迹tr(xiTWWTxi)tr\left( {x_i^TW{W^T}{x_i}} \right)tr(xiTWWTxi)来代替，(1)式就变成了：

(2)−tr(∑i=1mxiTWWTxi)+∑i=1mxiTxi=tr(AB)=tr(BA)−tr(WT(∑i=1mxixiT)W)+∑i=1mxiTxi=−tr(WTXXTW)+∑i=1mxiTxi\begin{array}{l}- tr\left( {\sum\limits_{i = 1}^m {x_i^TW{W^T}{x_i}} } \right) + \sum\limits_{i = 1}^m {x_i^T{x_i}} \\ \mathop = \limits^{tr\left( {AB} \right) = tr\left( {BA} \right)} - tr\left( {{W^T}\left( {\sum\limits_{i = 1}^m {{x_i}x_i^T} } \right)W} \right) + \sum\limits_{i = 1}^m {x_i^T{x_i}} \\ = - tr\left( {{W^T}X{X^T}W} \right) + \sum\limits_{i = 1}^m {x_i^T{x_i}} \end{array} \tag {2} −tr(i=1∑mxiTWWTxi)+i=1∑mxiTxi=tr(AB)=tr(BA)−tr(WT(i=1∑mxixiT)W)+i=1∑mxiTxi=−tr(WTXXTW)+i=1∑mxiTxi(2)

由于xix_ixi的维度是d×1d \times 1d×1，∑i=1mxiTxi\sum\limits_{i = 1}^m {x_i^T{x_i}}i=1∑mxiTxi是一个常量，且WWW的每一个向量wjw_jwj是标准正交基，因此，最小化上式中的距离等价于：

(3)arg⁡min⁡⎵W−tr⁡(WTXXTW)s.t.WTW=I\underbrace { \arg \min } _ { W } - \operatorname { tr } \left( W ^ { T } X X ^ { T } W \right) \text { s.t. } W ^ { T } W = I\tag {3} Wargmin−tr(WTXXTW)s.t.WTW=I(3)

（二）基于最大可分性推导PCA

样本点xix_ixi在空间中超平面上的投影为WTxiW^Tx_iWTxi，如果要让所有样本点的投影尽可能分开，应使投影后的方差最大化，如图所示。

使所有样本的投影尽可能分开，则需最大化投影点的方差

由于投影值的协方差矩阵的对角线代表了投影值在各个维度上的方差，则所有维度上的方差和可写成协方差矩阵的迹：（这个过程与第2节中描述的一致）

(4)∑i=1m(WTxi)(WTxi)T=∑i=1mWTxixiTW=tr(WTXXTW)\begin{array}{l} \sum\limits_{i = 1}^m {\left( {{W^T}{x_i}} \right){{\left( {{W^T}{x_i}} \right)}^T}} = \sum\limits_{i = 1}^m {{W^T}{x_i}x_i^TW} \\ = tr\left( {{W^T}X{X^T}W} \right) \end{array}\tag {4} i=1∑m(WTxi)(WTxi)T=i=1∑mWTxixiTW=tr(WTXXTW)(4)

因此，最大化方差等价于：

(5)arg⁡max⁡⎵Wtr⁡(WTXXTW)s.t.WTW=I\underbrace { \arg \max } _ { W } \operatorname { tr } \left( W ^ { T } X X ^ { T } W \right) \text { s.t. } W ^ { T } W = I\tag {5} Wargmaxtr(WTXXTW)s.t.WTW=I(5)

(5)式与(3)式等价。

5. PCA求解

对(3)式和(5)式中的优化目标，利用拉格朗日乘子法可得：

(6)J(W)=−tr(WTXXTW)+λ(WTW−I)J(W) = -tr \left( {{W^T}X{X^T}W} \right) + \lambda \left( {{W^T}W - I} \right)\tag {6} J(W)=−tr(WTXXTW)+λ(WTW−I)(6)

对XXX求导，由于∂tr(ATB)∂A=B\frac{{\partial tr\left( {A^TB} \right)}}{{\partial A}} = B∂A∂tr(ATB)=B，可得：

(7)∂J(W)∂W=−XXTW+λW\frac{{\partial J(W)}}{{\partial W}} = -X{X^T}W + \lambda W\tag {7} ∂W∂J(W)=−XXTW+λW(7)

令导数为0，得：

(8)XXTW=λWX{X^T}W = \lambda W\tag {8}XXTW=λW(8)

根据线性代数中的特征值分解Ax=λxAx = \lambda xAx=λx可知上式是一个类似的问题。于是，只需要对协方差矩阵XXTX{X^T}XXT进行特征值分解，将求得的特征值排序：λ1≥λ2≥…≥λd\lambda _ { 1 } \geq \lambda _ { 2 } \geq \ldots \geq \lambda _ { d }λ1≥λ2≥…≥λd，取前d′d'd′个特征值对应的特征向量构成W=(w1,w2,…,wd′)W = \left( w _ { 1 } , w _ { 2 } , \dots , w _ { d ^ { \prime } } \right)W=(w1,w2,…,wd′)，这就是主成分分析的解。

降维后低维空间的维数d′d'd′通常是事先指定的，还可以设置一个重构阈值，例如t=95t=95%t=95，然后选取使下式成立的最小d′d'd′值：

(9)∑i=1d′λi∑i=1dλi≥t\frac { \sum _ { i = 1 } ^ { d ^ { \prime } } \lambda _ { i } } { \sum _ { i = 1 } ^ { d } \lambda _ { i } } \geq t\tag {9} ∑i=1dλi∑i=1d′λi≥t(9)

6. PCA算法描述

输入：训练数据集D={x1,x2,⋯ ,xm}D= \left\{ x _ { 1 } , x _ { 2 } , \cdots , x _ { m } \right\}D={x1,x2,⋯,xm}，低维空间维数d′d'd′。

过程：

（1）对所有样本进行中心化：

xi=xi−1m∑j=1mxjx_i = x_i - \frac { 1 } { m } \sum _ { j = 1 } ^ { m } x _j xi=xi−m1j=1∑mxj

（2）计算样本的协方差矩阵XXTX{X^T}XXT；

（3）对协方差矩阵XXTX{X^T}XXT做特征值分解，求出其特征值及其对应的特征向量；

（4）取最大的d′d'd′个特征值对应的特征向量(w1,w2,…,wd′)\left( w _ { 1 } , w _ { 2 } , \dots , w _ { d ^ { \prime } } \right)(w1,w2,…,wd′)；

（5）对样本集DDD中的每一个样本xix_ixi，转化为新样本zi=WTxiz_i=W^Tx_izi=WTxi，得到输出样本集D′={z1,z2,⋯ ,zm}D'= \left\{ z _ { 1 } , z _ { 2 } , \cdots , z_ { m } \right\}D′={z1,z2,⋯,zm}。

输出：投影矩阵W=(w1,w2,…,wd′)W = \left( w _ { 1 } , w _ { 2 } , \dots , w _ { d ^ { \prime } } \right)W=(w1,w2,…,wd′)和降维后的样本集D′={z1,z2,⋯ ,zm}D'= \left\{ z _ { 1 } , z _ { 2 } , \cdots , z_ { m } \right\}D′={z1,z2,⋯,zm}。