PCA主成分分析

主成分分析(PCA)是一种常用的无监督学习方法，这一方法利用正交变换把由线性相关变量表示的观测数据转换为少数几个由线性无关变量表示的数据，线性无关的变量称为主成分。

主成分的个数通常小于原始变量的个数，所以主成分分析属于降维方法。

主成分分析主要用于发现数据中的基本结构，即数据中变量之间的关系。

总体主成分分析

主成分分析中，

1. 首先对给定数据进行规范化，使得数据每一变量的平均值为0，方差为1

2. 对数据进行正交变换。原来由线性相关变量表示的数据，通过正交变换变成由若干个线性无关的新变量表示的数据。新变量是可能的正交变换中变量的方差的和(信息保存)最大的，方差表示在新变量上的信息的大小。

   将新变量依次称为第一主成分、第二主成分等。

通过主成分分析，可以利用主成分近似地表示原始数据(可理解为发现数据的基本结构)，也可以把数据由少数主成分表示(这可理解为对数据降维)。

数据集合中的样本由实数空间(正交坐标系)中的点表示，空间的一个坐标轴表示一个变量，规范化处理后得到的数据分布在原点附近。对原坐标系中的数据进行主成分分析等价于进行坐标系旋转变换，将数据投影到新坐标系的坐标轴上；新坐标系的第一坐标轴、第二坐标轴等分别表示第一主成分、第二主成分等。数据在每一轴上的坐标值的平方表示相应变量的方差；并且，这个坐标系是在所有可能的新的坐标系中，坐标轴上的方差的和最大的。

如下图所示，原坐标系(X1,X2)，进行正交变换后新坐标系(Y1,Y2)：

接下来，看看方差最大的解释：

假设有两个变量X1,X2，三个样本点A,B,C，样本分布在由X1,X2轴组成的坐标系中，如下图所示：

对坐标系进行旋转变换，得到新的坐标轴Y1，表示新的变量Y1。样本点A,B,C在Y1轴上投影，得到Y1轴的坐标值A′,B′,C′。坐标值的平方和OA′2+OB′2+OC′2表示样本在变量y1上的方差。

主成分分析旨在选取正交变换中方差最大的变量，作为第一主成分，也就是旋转变换中的坐标值的平方和最大的轴。

如何确保平方和最大？首先旋转变换中样本点到原点的距离的平方和OA2+OB2+OC2保持不变，根据勾股定理，新坐标值的平方和OA′2+OB′2+OC′2最大等价于样本点到Y1轴的距离的平方和AA′2+BB′2+CC′2最小。也就是说，主成分分析在旋转变换中选取离样本点的距离平方和最小的轴，作为第一主成分。

第二主成分等的选取，在保证与已选新坐标轴(已选主成分)正交的条件下，类似的进行。

在数据总体上进行主成分分析称为总体主成分分析。

在有限样本上进行主成分分析称为样本主成分分析。

定义和导出

假设x=(x1,x2,…,xm)⊤是m维随机变量，其均值向量是μ

μ=E(x)=(μ1,μ2,…,μm)⊤

协方差矩阵是Σ

Σ=cov(x,x)=E[(x−μ)(x−μ)⊤]

考虑由m维随机变量x 到m维随机变量y=(y1,y2,…,ym)⊤的线性变换：

yi=α⊤ix=α1ix1+α2ix2+⋯+αmixm

其中α⊤i=(α1i,α2i,…,αmi) ,i=1,2,…,m。

由随机变量的性质可知：

E(yi)=α⊤iμ ,i=1,2,…,mvar(yi)=α⊤iΣαi ,i=1,2,…,mcov(yi,yj)=α⊤iΣαj ,i=1,2,…,m ;j=1,2,…,m

定义16.1(总体主成分) 给定一个如式(16-1)所示的线性变换，如果它们满足下列条件

1. 系数向量α⊤i是单位向量，即α⊤iαi=1 ,i=1,2,…,m

2. 变量yi与yj 互不相关，即cov(yi,yj)=0 (i≠j)

3. 变量y1是x的所有线性变换中方差最大的；y2是与y1不相关的x的所有线性变换中方差最大的；一般地，yi是与y1,y2,…,yi−1 (i=1,2,…,m)都不相关的x的所有线性变换中方差最大的；

   这时分别称y1,y2,…,ym (i=1,2,…,m)为x的第一主成分、第二主成分、…、第m主成分。

定义中的条件1表明线性变换是正交变换，α1,α2,…,αm是一组标准正交基，

α⊤iαj={1,i=j0,i≠j

条件2，3给出了一个求主成分的方法：

第一步：在x的所有线性变换

α⊤1x=m∑i=1αi1xi

中，在α⊤1α1=1条件下，求方差最大的，得到x的第一主成分；

第二步，在与α⊤1x不相关的x的所有线性变换：

α⊤2x=m∑i=1αi2xi

中，在 α⊤2α2=1条件下，求方差最大的，得到x的第二主成分；

第k步，在于α⊤1x,α⊤2x,…,α⊤k−1x不相关的x的所有线性变换：

α⊤kx=m∑i=1αikxi

中，在α⊤kαk=1条件下，求方差最大的，得到x的第k主成分；

如此继续下去，直到得到x的第m主成分。

定理16.1 设x是m维随机变量，Σ是x的协方差矩阵，Σ的特征值分别是λ1≥λ2≥⋯≥λm≥0，特征值对应的单位特征向量分别是α1,α2,…,αm，则x的第k主成分是:

yk=α⊤kx=α1kx1+α2kx2+⋯+αmkxm , k=1,2,…,m

关于x的第k主成分的方差是:

var(yk)=α⊤kΣαk=λk , k=1,2,…,m

即协方差矩阵Σ的第k个特征值。

若特征值有重根，对应的特征向量组成m维空间Rm的一个子空间，子空间的维数等于重根数，在子空间任取一个正交坐标系，这个坐标系的单位向量就可作为特征向量。这是坐标系取法不唯一。

证明：采用拉格朗日乘子法求出主成分。

首先，求x的第一主成分y1=α⊤1x，

即求系数向量α1。由定义16.1知，第一主成分的α1是在α⊤1α1=1条件下，x的所有线性变换中使方差：

var(α⊤1x)=α⊤1Σα1

达到最大的。

求第一主成分就是求解约束最优化问题：

maxα1α⊤1Σα1 α⊤1α1=1

定义拉格朗日函数

α⊤1Σα1−λ(α⊤1α1−1)

其中λ为拉格朗日乘子。将拉格朗日函数对α1进行求导，并令其为0，得：

Σα1−λα1=0

因此，λ是Σ得特征值，α1是对应的单位特征向量。于是，目标函数：

α⊤1Σα1=α⊤1λα1=λα⊤1α1=λ

假设α1是Σ的最大特征值λ1对应的单位特征向量，显然α1与λ1是最优化问题的解。

所以，α⊤1x构成第一主成分，其方差等于协方差矩阵的最大特征值:

var(α⊤1x)=α⊤1Σα1=λ1

接着，求x的第二主成分y2=α⊤2x，

第二主成分的α2是在α⊤2α2=1，且α⊤2x与α⊤1x不相关的条件下(也就是说α⊤2x与α⊤1x的协方差为0，协方差公式(16-4))，x的所有线性变换中使方差：

var(α⊤2x)=α⊤2Σα2

达到最大的。

求第二主成分需要求解约束最优化问题:

maxα2=α⊤2Σα2 α⊤1Σα2=0,α⊤2Σα1=0,α⊤2α2=1 α⊤1Σα2=α⊤2Σα1=α⊤2λ1α1=λ1α⊤2α1=λ1α⊤1α2 α⊤1α2=0,α⊤2α1=0

定义拉格朗日函数

α⊤2Σα2−λ(α⊤2α2−1)−ϕα⊤2α1

其中λ，ϕ是拉格朗日乘子。对α2求导，并令其为0，得:

2Σα2−2λα2−ϕα1=0

将方程左乘以α⊤1有:

2α⊤1Σα2−2λα⊤1α2−ϕα⊤1α1=0

此式前两项为0，且α⊤1α1=1，导出ϕ=0，因此式(16-10)成为:

Σα2−λα2=0

由此，λ是Σ的特征值，α2是对应的单位特征向量。于是，目标函数：

α⊤2Σα2=α⊤2λα2=λα⊤2α2=λ

假设α2是Σ的第二大特征值λ2对应的单位特征向量，显然α2与λ2是以上最优化问题的解。

于是α⊤2x构成第二主成分，其方差等于协方差矩阵的第二大特征，

var(α⊤2x)=α⊤2Σα2=λ2

一般地，x的第k主成分是α⊤kx，并且var(α⊤kx)=λk，这里λk是Σ的第k个特征值并且αk是对应的单位特征向量。

以此类推，第k个主成分的方差等于Σ的第k个特征值:

var(α⊤kx)=α⊤kΣαk=λk ,k=1,2,…,m

定理证毕。

推论16.1 m维随机变量y=(y1,y2,…,ym)⊤的分量依次是x 的第一主成分到第m主成分的充要条件是：

1. 令y=A⊤x，A为正交矩阵

   A=[α11α12…α1mα21α22…α2m⋮⋮⋮αm1αm2…αmm]

2. 其中y的协方差矩阵为对角矩阵

   cov(y)=diag(λ1,λ2,…,λm)λ1≥λ2≥⋯≥λm

   其中λk是Σ的第k个特征值，αk是对应的单位特征向量，k=1,2,…,m。

以上证明中，λk是Σ的第k个特征值，αk是对应的单位特征向量，即

Σαk=λkαk ,k=1,2,…,m

用矩阵表示即为:

ΣA=AΛ

这里A=[αij]m×m(注意：这里的A与前面的A∈m×n 不是同一个！)，Λ是对角矩阵，其第k个对角元素是λk。因为A是正交矩阵，即A⊤A=AA⊤=I，由式(16-14)得到两个公式:

A⊤ΣA=Λ Σ=AΛA⊤

下面叙述总体主成分的性质：

1. 总体主成分y的协方差矩阵是对角矩阵:

   cov(y)=Λ=diag(λ1,λ2,…,λm)

2. 总体主成分y的方差之和等于随机变量x的方差之和，即：

   m∑i=1λi=m∑i=1σii

   其中σii是随机变量xi的方差，即协方差矩阵Σ的对角元素。事实上，利用式(16-16)及矩阵的迹(trace)的性质，可知：

   m∑i=1var(xi)=tr(Σ⊤)=tr(AΛA⊤)=tr(A⊤ΛA)=tr(Λ)=m∑i=1λi=m∑i=1var(yi)

3. 第k个主成分yk与变量xi的相关系数ρ(yk,xi)称为因子负荷量(factor loading)，它表示第k个主成分yk与变量xi的相关性。计算公式：

   ρ(yk,xi)=√λkαik√σii ,k,i=1,2,…,m ρ(yk,xi)=cov(yk,xi)√var(yk)var(xi)=cov(α⊤kx,e⊤ix)√λk√σii

   其中，ei为基本单位向量，其第i个分量为1，其余分量为0.再由协方差的性质：

   cov(α⊤kx,e⊤ix)=α⊤kΣei=e⊤iΣαk=λke⊤iαk=λkαik

   故得式(16-20)。

4. 第k个主成分yk与m个变量的因子负荷量满足：

   m∑i=1σiiρ2(yk,xi)=λk

   由式(16-20)有：

   m∑i=1σiiρ2(yk,xi)=m∑i=1λkα2ij=λkα⊤kαk=λk

5. m个主成分与第i个变量xi的因子负荷量满足：

   m∑k=1ρ2(yk,xi)=1

   由于y1,y2,…,ym互不相关，如：

   ρ2(xi,(y1,y2,…,ym))=m∑i=1ρ2(yk,xi)

   又因xi可以表为y1,y2,…,ym的线性组合，所以xi 与y1,y2,…,ym的相关系数的平方为1，式(16-22)得证。

主成分的个数

主成分分析的主要目的是降维，所以一般选择k(k≪m)个主成分(线性无关变量)来代替m个原有变量(线性相关变量)，使问题得以简化，并能保留原有变量的大部分信息。

信息是指原有变量的方差。

定理16.2 对任意正整数q，1≤q≤m，考虑正交线性变换

y=B⊤x

其中y是q维向量，B⊤是qxm矩阵，令y的协方差矩阵为：

Σy=B⊤ΣB

则Σy的迹tr(Σy)在B=Aq时取得最大值，其中矩阵Aq由正交矩阵A∈m×m的前q列组成。

证明：

令βk是B的第k列，由于正交矩阵A的列构成m维空间的基，所以βk可以由A的列表示，即：

βk=m∑j=1cjkαj B=AC

其中C是mxq矩阵，其第j行第k列元素为cjk。

B⊤ΣB=C⊤A⊤ΣAC=C⊤ΛC=m∑i=1λjcjc⊤j

其中c⊤j是C的第j行。因此

tr(B⊤ΣB)=m∑j=1λjtr(cjc⊤j)=m∑j=1λjtr(c⊤jcj)=m∑j=1λjc⊤jcj=m∑j=11∑k=1λjc2jk

其次，由式(16-25)及A的正交性知:

C=A⊤B

由于A是正交的，B的列是正交的，所以:

C⊤C=B⊤AA⊤B=B⊤B=Iq

即C的列也是正交的。于是

tr(C⊤C)=tr(Iq)m∑j=1q∑k=1c2jk=q

这样，矩阵C可以认为是某个m阶正交矩阵D的前q列。正交矩阵D的行也正交，所以满足:

d⊤jdj=1 ,j=1,2,…,m

其中d⊤j是D的第j行。由于矩阵D的行包括矩阵C的行的前q个元素，所以：

c⊤jcj≤1 ,j=1,2,…,m q∑k=1c2jk≤1 ,j=1,2,…,m

式(16-26)中∑qk=1c2jk是λj的系数，由式(16-27)这些系数之和是q，且由式(16-28)知这些系数小于等于1.因为λ1≥λ2≥⋯≥λq≥⋯≥λm，显然，当能找到cjk使得

q∑k=1c2jk={1,j=1,2,…,q0,j=q+1,…,m

时，∑mj=1(∑qk=1c2jk)λj最大。而当B=Aq时，有

cjk={1,1≤j=k≤q0,其它

满足式(16-29)。所以，当B=Aq时，tr(Σy)达到最大值。

定理16.2表明，当x的线性变换y在B=Aq时，其协方差矩阵Σy的迹tr(Σy)取得最大值，也就是说，当取A的前q列(取x的前q个主成分)时，能够最大限度地保留原有变量方差的信息。

定理16.3 考虑正交变换

y=B⊤x

这里B⊤是pxm矩阵，A和Σy的定义与定理16.2相同，则tr(Σy)在B=Ap时，取得最小值，其中矩阵Ap由A的后p列组成。

定理16.3可以理解为：当舍弃A的后p列，即舍弃变量x的后p个主成分时，原有变量的方差损失最少。

以上两个定理可以作为选择k个主成分的理论依据。具体选择k的方法，通常利用方差贡献率。

定义16.2 第k主成分yk的方差贡献率定义为yk的方差与所有方差之和的比，记作ηk：

ηk=λk∑mi=1λi

k个主成分y1,y2,…,yk的累计方差贡献率定义为k个方差之和与所有方差之和的比：

k∑i=1ηi=∑ki=1λi∑mi=1λi

通常取k使得累计方差贡献率达到某个阈值，比如70%～80%以上。累计方差贡献率反映了主成分保留信息的比例。

但不能反映对某个原有变量xi的保留信息的比例，这时通常利用k个主成分y1,y2,…,yk对原有变量xi的贡献率。

定义16.3 k个主成分y1,y2,…,yk对原有变量xi的贡献率定义为xi与(y1,y2,…,yk)的相关系数的平方，记作vi

vi=ρ2(xi,(y1,y2,…,yk))

计算公式如下：

vi=ρ2(xi,(y1,y2,…,yk))=k∑j=1ρ2(xi,yj)=k∑j=1λjα2ijσii

规范化变量的总体主成分

在实际问题中，不同变量可能有不同的量纲，直接求主成分有时会产生不合理的结果。为了消除这个影响，需要对各个随机变量实施规范化，使其均值为0，方差为1.

设x=(x1,x2,…,xm)⊤为m维随机变量，xi为第i个随机变量，i=1,2,…,m，令

x∗i=xi−E(xi)√var(xi) ,i=1,2,…,m

其中，E(xi),var(xi)分别是随机变量xi的均值和方差，这时x∗i就是xi的规范化随机变量。

显然，规范化随机变量的协方差矩阵就是相关矩阵R。

对照总体主成分的性质可知，规范化随机变量的总体主成分分析有以下性质：

1. 规范化变量主成分的协方差矩阵是

   Λ∗=diag(λ∗1,λ∗2,…,λ∗m)

   其中λ∗1≥λ∗2≥⋯≥λ∗m≥0为相关矩阵R的特征值。

2. 协方差矩阵的特征值之和为m

   m∑k=1λ∗k=m

3. 规范化随机变量x∗i与主成分y∗k的相关系数(因子负荷量)为

   ρ(y∗k,x∗i)=√λ∗ke∗ik ,k,i=1,2,…,m

   其中e∗k=(e∗1k,e∗2k,…,e∗mk)⊤为矩阵R对应的特征值λ∗k对应的单位特征向量。

4. 所有规范化随机变量x∗i与主成分y∗k的相关系数的平方和等于λ∗k

   m∑i=1ρ2(y∗k,x∗i)=m∑i=1λ∗ke∗2ik=λ∗k ,k=1,2,…,m

5. 规范化随机变量x∗i与所有主成分y∗k的相关系数的平方和等于1

   m∑k=1ρ2(y∗k,x∗i)=m∑k=1λ∗ke∗2ik=1 ,k=1,2,…,m

样本主成分分析

16.1是定义在样本总体上的主成分分析。

在实际问题中，需要在观测数据上进行主成分分析，也就是样本主成分分析。

样本主成分的定义和性质

假设对m维随机变量x=(x1,x2,…,xm)⊤进行n次独立观测，x1,x2,…,xn表示观测样本，其中 xj=(x1j,x2j,…,xmj)⊤表示第j个样本，xij表示第j个观测样本的第i个变量，j=1,2,…,n。观测数据用样本矩阵X表示，记作：

X=[x1,x2,…,xn]=[x11x12…x1nx21x22…x2n⋮⋮⋮xm1xm2…xmn]

给定样本矩阵X，可以估计样本均值，以及样本协方差。样本均值向量 ˉx为

ˉx=1nn∑j=1xj

样本协方差矩阵S为：

S=[sij]m×msij=1n−1n∑k=1(xik−ˉxi)(xjk−ˉxj) ,i,j=1,2,…,m

其中ˉxi=1n∑nk=1xik为第i个变量的样本均值，ˉxj=1n∑nk=1xjk为第j个变量的样本均值。

样本相关矩阵R为:

R=[rij]m×m ,rij=sij√siisjj ,i,j=1,2,…,m

定义m维向量x=(x1,x2,…,xm)⊤到m维向量y=(y1,y2,…,ym)⊤的线性变换

y=A⊤x A=[a1,a2,…,am]=[a11a12…a1ma21a22…a2m⋮⋮⋮am1am2…amm]ai=(a1i,a2i,…,ami)⊤ ,i=1,2,…,m

考虑式(16-43)的任意一个线性变换(对照式(16-1))

yi=a⊤ix=a1ix1+a2ix2+⋯+amixm

其中yi是m维向量y的第i个变量，相应于容量为n的样本x1,x2,…,xn，yi的样本均值ˉyi为：

ˉyi=1nn∑j=1a⊤ixj=α⊤iˉx

其中ˉx是随机变量x的样本均值：

ˉx=1nn∑j=1xj

其中，yj的样本方差var(yi)为：

var(yi)=1n−1n∑j=1(a⊤ixj−a⊤iˉx)2=a⊤i[1n−1n∑j=1(xj−ˉx)(xj−ˉx)⊤]ai=c

对任意两个线性变换yi=a⊤ix,yk=a⊤kx，对应于容量为n的样本x1,x2,…,xn，yi,yk的样本协方差为:

cov(yi,yk)=a⊤iSak

定义16.4(样本主成分) 给定样本矩阵X。样本第一主成分y1=a⊤1x是在a⊤1a1=1条件下，使得a⊤1xj (j=1,2,…,n)的样本方差a⊤1Sa1最大的x的线性变换；样本第二主成分y2=a⊤2x是在a⊤2a2=1和a⊤2xj与a⊤1xj (j=1,2,…,n)的样本协方差a⊤1Sa2=0条件下，使得a⊤2xj (j=1,2,…,n)的样本方差a⊤2Sa2最大的x的线性变换；**一般地，样本在第i主成分yi=a⊤ix是在a⊤iai=1和a⊤ixj与a⊤kxj (k<i , j=1,2,…,n)的样本协方差a⊤kSai=0条件下，使得a⊤ixj (j=1,2,…,n)的样本方差a⊤iSai最大的x的线性变换。

在使用样本主成分时，一般假设样本数据是规范化的，即对样本矩阵作如下变换：

x∗ij=xij−ˉxi√sii ,i=1,2,…,m ;j=1,2,…,n ˉxi=1nn∑j=1xij ,i=1,2,…,msii=1n−1n∑j=1(xij−ˉxi)2 ,i=1,2,…,m

这时，样本协方差矩阵S就是样本相关矩阵R：

R=1n−1XX⊤

样本协方差矩阵S是总体协方差矩阵Σ的无偏估计，样本相关矩阵R是总体相关矩阵的无偏估计，S的特征值和特征向量是Σ 特征值和特征向量的极大似然估计。

相关矩阵的特征值分解算法

给定样本矩阵X，利用数据的样本协方差矩阵或者样本相关矩阵的特征值分解进行主成分分析。具体步骤如下：

1. 对观测数据按式(16-48)进行规范化处理，仍以X表示

2. 依据规范化数据矩阵，计算样本相关矩阵R

   R=[rij]m×m=1n−1XX⊤ rij=1n−1n∑l=1xilxlj , i,j=1,2,…,m

3. 求样本相关矩阵R的k个特征值和对应的k个单位特征向量。

   求解R的特征方程

   |R−λI|=0

   得R的m个特征值

   λ1≥λ2≥⋯≥λm

   求方差贡献率∑ki=1ηi达到预定值的主成分个数k。

   求前k个特征值对应的单位特征向量：

   ai=(a1i,a2i,…,ami)⊤ ,i=1,2,…,k

4. 求k个样本主成分

   以k个单位特征向量为系数进行线性变换，求出k个样本主成分：

   yi=a⊤ix ,i=1,2,…,k

5. 计算k个主成分yj与原变量xi的相关系数ρ(xi,xj)，以及k个主成分对原变量xi的贡献率vi。

6. 计算n个样本的k个主成分值

   将规范化样本数据代入k个主成分式(16-50)，得到n个样本的主成分值。第j个样本(x)j=(x1j,x2j,…,xmj)⊤的第i主成分值是：

   yij=(a1i,a2i,…,ami)(x1j,x2j,…,xmj)⊤=m∑l=1alixlj i=1,2,…,m ,j=1,2,…,n

数据矩阵的奇异值分解算法

给定样本矩阵X，利用数据矩阵奇异值分解进行主成分分析。具体过程如下。这里假设有k个主成分。

参照式(15-19)，对于mxn实矩阵A，假设其秩为r，0<k<r，则可以将矩阵A进行截断奇异值分解：

A≈UkΣkV⊤k

式中Uk是mxk矩阵，Vk是nxk矩阵，Σk是k阶对角矩阵；Uk，Vk分别由取A的完全奇异值分解的矩阵U,V的前k列，Σk由取A的完全奇异值分解的矩阵Σ的前k个对角线元素得到。

定义一个新的nxm矩阵X′

X′=1√n−1X⊤

且X′的每一个列均值为零。不难得知，

X′⊤X′=(1√n−1X⊤)⊤(1√n−1X⊤)=1√n−1XX⊤

即X′⊤X′等于X的协方差矩阵SX。

SX=X′⊤X′

主成分分析归结于求协方差矩阵SX的特征值和对应的单位特征向量，所以问题转化为求矩阵X′⊤X′的特征值和对应的单位特征向量。

假设X′的截断奇异值分解为X′=UΣV⊤，那么V的列向量就是SX=X′⊤X′的单位特征向量。因此，V的列向量就是X的主成分。

于是，求X主成分可以通过求X′的奇异值分解来实现。具体算法如下：

算法16.1(主成分分析算法)

输入：mxn样本矩阵X，其每一行元素的均值为零

输出：kxn样本主成分矩阵Y

参数：主成分个数k

1. 构造新的nxm矩阵

   X′=1√n−1X⊤

   且X′每一列的均值为零。

2. 对矩阵X′进行截断奇异值分解，得到：

   X′=UΣV⊤

   有k个奇异值、奇异向量。矩阵V 的前k列构成k个样本主成分。

3. 求kxn样本主成分矩阵