文章目录

  • 问题设置(Problem Setting)
  • 最大化方差的角度看PCA(Maximum Variance Perspective)
    • 最大方差的方向(Direction with Maximal Variance)
    • M维子空间下的最大方差(M-dimensional Subspace with Maximal Variance)
  • 投影的角度看待PCA(Projection Perspective)
    • 问题设置和问题目标(Setting and Objective)
    • 找到最优化坐标(Finding Optimal Coordinates)
    • 找到主子空间的基向量(Finding the Basis of the Principal Subspace)
  • 特征向量计算以及低秩逼近(Eigenvector Computation and Low-Rank Approximations)
    • 用低秩逼近的PCA(PCA Using Low-Rank Matrix Approximations)
    • 实际方面(Practical Aspects)
  • 高维PCA(PCA in High Dimensions)
  • PCA在实践中的关键步骤(Key Steps of PCA in Practice)
  • 用潜变量看待PCA(Latent Variable Perspective)
    • 生成过程及概率模型(Generative Process and Probabilistic Model)
    • 似然以及联合分布(Likelihood and Joint Distribution)
    • 后验分布(Posterior Distibution)
  • 拓展阅读

对于一些高维的数据,分析难度大,而且想要对这些数据进行可视化几乎是不可能的,并且想要存储这些数据的代价也是及其昂贵的,所以我们想要找到一种能够将数据的维度降低的方法。这其中, 主成分分析法(principal component analysis (PCA))是最常用的方法之一。

问题设置(Problem Setting)

在PCA中,我们希望能够找到一个一个向量的投影向量x~n\tilde x_nx~n​,与原向量尽可能相近。
对于一个独立均匀分布的数据集X={x1,⋯,xN},xn∈RD\mathcal X=\{x_1,\cdots,x_N\}, x_n\in \mathbb R^DX={x1​,⋯,xN​},xn​∈RD,它的均值为0, 对应的数据方差矩阵为:
S=1N∑n=1Nxnx⊤S=\frac 1N \sum^N_{n=1}x_nx^\topS=N1​n=1∑N​xn​x⊤
压缩之后表示为:
zn=B⊤xn∈RMz_n = B^\top x_n\in \mathbb R^Mzn​=B⊤xn​∈RM
其中,B为投影矩阵,定义为:
B:=[b1,⋯,bM]∈RD×MB := [b_1,\cdots,b_M]\in \mathbb R^{D\times M}B:=[b1​,⋯,bM​]∈RD×M
假设bib_ibi​为正交规范基,则bi⊤bj=0,i≠j;bi⊤bi=1b_i^\top b_j=0, i\ne j;b_i^\top b_i=1bi⊤​bj​=0,i​=j;bi⊤​bi​=1.我们希望找到一个M维的子空间U⊆RD,dim⁡(U)=M<DU\subseteq \mathbb R^D,\operatorname{dim}(U)=M < DU⊆RD,dim(U)=M<D,向其中的投影的向量与原先的向量最相似,因为压缩造成的损失最小。我们将投影的数据表示为x~n∈U\tilde x_n\in Ux~n​∈U,对应的坐标为znz_nzn​(基向量为b1,⋯,bMb_1,\cdots,b_Mb1​,⋯,bM​)

PCA的目标是最小化平方重构误差(the Squared Reconstruction Error)∥xn−x~n∥2\| x_n-\tilde x_n\|^2∥xn​−x~n​∥2.
从数据压缩的角度来看,我们先是将源数据压缩到一个更低维度的空间中,对应下图的zzz,然后将压缩的信息复原,对应下图中的x~\tilde xx~; zzz控制着多少信息能够从xxx到x~\tilde xx~.在PCA中,我们考虑原始数据与低维数据之间的线性关系,所以有以下关系:
z=B⊤x;x~=Bzz = B^\top x;\tilde x = Bzz=B⊤x;x~=Bz。将PCA看成是一个数据压缩的过程,所以可以认为第一个箭头是编码器(encoder),第二个箭头是解码器(decoder)

Graphical illustration of PCA. In PCA, we find a compressed version z of original data x. The compressed data can be reconstructed into x~\tilde xx~, which lives in the original data space, but has an intrinsic lower-dimensional representation than xxx.

最大化方差的角度看PCA(Maximum Variance Perspective)

在下图中,我们丢弃了数据关于x2x_2x2​的信息,这样做能够达到降维的效果,而且使得数据的损失最小化,是源数据与降维之后的数据尽可能相似。假设忽略x1x_1x1​的信息,则得到的数据就很不相似了,也就是说这个降维操作导致了很多的信息损失。通过观察可以发现,数据在两个维度上的分散程度不一样。当数据在一个维度上越分散,说明这个维度上所包含的信息也就越多,而方差可以表示数据分散程度的大小,所以从方差的角度理解PCA就是找到低维空间中数据方差最大的维度

为了运算方便,我们对数据进行一个均值归一化(Mean Normalization),因为我们要研究的是方差,而对数据整体的几何运算并不会影响数据的方差
Vz[z]=Vx[B⊤(x−μ)]=Vx[B⊤x−B⊤μ]=Vx[B⊤x]\mathbb V_z[z]=\mathbb V_x[B^\top(x-\mu)]=\mathbb V_x[B^\top x - B^\top \mu]=\mathbb V_x[B^\top x]Vz​[z]=Vx​[B⊤(x−μ)]=Vx​[B⊤x−B⊤μ]=Vx​[B⊤x]
这时候对应的低维空间的数据的均值也是0:Ez[z]=Ex[B⊤x]=B⊤Ex[x]=0\mathbb E_z[z]=\mathbb E_x[\boldsymbol B^\top \boldsymbol x]=\boldsymbol B^\top \mathbb E_x[\boldsymbol x]=\boldsymbol 0Ez​[z]=Ex​[B⊤x]=B⊤Ex​[x]=0。

B\boldsymbol BB是投影矩阵,将源数据投影到主成分上,从而实现降维。

最大方差的方向(Direction with Maximal Variance)

为了找到数据在低维空间中的最大的方差,我们先找到一个向量b1∈RDb_1\in \mathbb R^Db1​∈RD,数据在这个向量上的投影的方差最大,也就是要最大化z∈RMz\in \mathbb R^Mz∈RM中第一个坐标z1z_1z1​的方差:
V1:=V[z1]=1N∑n=1Nz1n2V_1 := \mathbb V[z_1]=\frac1N\sum^N_{n=1}z^2_{1n}V1​:=V[z1​]=N1​n=1∑N​z1n2​
我们将z1nz_{1n}z1n​表示为数据(xn∈RDx_n\in \mathbb R^Dxn​∈RD)在低维空间(zn∈RMz_n\in \mathbb R^Mzn​∈RM)的第一个坐标。znz_nzn​的第一个成分为:
z1n=b1⊤xnz_{1n}=b_1^\top x_nz1n​=b1⊤​xn​
这是xnx_nxn​在b1b_1b1​张成的一维子空间中的正交投影,将上面二式联立:
V1=1N∑n=1N(b1⊤xn)2=1N∑n=1Nb1⊤xnxn⊤b1=b1⊤(1N∑n=1Nxnxn⊤)b1=b1⊤Sb1V_1=\frac1N\sum^N_{n=1}(b_1^\top x_n)^2=\frac 1N \sum^N_{n=1}b_1^\top x_n x_n^\top b_1=b_1^\top\begin{pmatrix}\frac 1N\sum\limits^N_{n=1}x_nx_n^\top \end{pmatrix}b_1=b_1^\top Sb_1V1​=N1​n=1∑N​(b1⊤​xn​)2=N1​n=1∑N​b1⊤​xn​xn⊤​b1​=b1⊤​(N1​n=1∑N​xn​xn⊤​​)b1​=b1⊤​Sb1​
其中, S为数据协方差矩阵。由上式可知,正交基(bib_ibi​)会对最终的方差的结果产生影响,所以这里要求这些基向量为规范正交基(∥b1∥2=1\|b_1\|^2=1∥b1​∥2=1),这样问题就转换成一个约束问题:
max⁡b1b1⊤Sb1,s.t.∥b1∥2=1\max_{b_1}b_1^\top Sb_1,\quad s.t.\ \|b_1\|^2 = 1b1​max​b1⊤​Sb1​,s.t. ∥b1​∥2=1
利用拉格朗日方法:
L(b1,λ)=b1⊤Sb1+λ1(1−b1⊤b1)\mathfrak L(b_1,\lambda)=b_1^\top Sb_1+\lambda_1(1-b_1^\top b_1)L(b1​,λ)=b1⊤​Sb1​+λ1​(1−b1⊤​b1​)
对上式分别求偏导:
∂L∂b1=2b1⊤S−2λ1b1⊤,∂L∂λ1=1−b1⊤b1\frac{\partial \mathfrak L}{\partial b_1}=2b_1^\top S-2\lambda_1 b_1^\top,\quad \frac{\partial \mathfrak L}{\partial \lambda_1}=1-b_1^\top b_1∂b1​∂L​=2b1⊤​S−2λ1​b1⊤​,∂λ1​∂L​=1−b1⊤​b1​
令偏微分的结果为0:
Sb1=λ1b1b1⊤b1=1\begin{aligned}Sb_1&=\lambda_1 b_1 \\b_1^\top b_1 &= 1\end{aligned}Sb1​b1⊤​b1​​=λ1​b1​=1​
由上式可以知道,λ1\lambda_1λ1​是方差S的一个特征值,b1b_1b1​是一个特征向量,利用这个式子,我们可以将问题转化成:
V1=b1⊤Sb1=λ1b1⊤b1=λ1V_1=b_1^\top Sb_1 = \lambda_1b_1^\top b_1 = \lambda_1V1​=b1⊤​Sb1​=λ1​b1⊤​b1​=λ1​
所以我们需要找到一个特征值最大的特征向量,这样源数据在投影之后的方差最大,这个特征向量称为主成分(Principal Component)我们可以得到投影数据点:
x~n=b1z1n=b1b1⊤xn∈RD\tilde x_n=b_1 z_{1n}=b_1b_1^\top x_n\in \mathbb R^Dx~n​=b1​z1n​=b1​b1⊤​xn​∈RD
注意这里的投影点上的数据是高纬度空间中的数据,但是实际上存储的时候只需要用低纬度的空间信息就可以表示了。

M维子空间下的最大方差(M-dimensional Subspace with Maximal Variance)

m−1m-1m−1个主成分对应的是SSS的m−1m-1m−1个特征向量,这些特征向量对应着最大的m−1m-1m−1个特征值。由于S=1N∑n=1Nxnxn⊤S = \frac1N \sum\limits_{n=1}^Nx_nx_n^\topS=N1​n=1∑N​xn​xn⊤​,所以S是一个对称矩阵,所以由谱定理可以得知,这些特征向量能够形成RD\mathbb R^DRD空间下的m−1m-1m−1维子空间的正交规范特征基。想要找到这些正交基,可以使用向量减法:
X~:=X=∑i=1m−1bibi⊤X=X−Bm−1X\tilde X := X=\sum^{m-1}_{i=1}b_ib_i^\top X=X-B_{m-1}XX~:=X=i=1∑m−1​bi​bi⊤​X=X−Bm−1​X
其中,数据点的列向量X=[x1,⋯,xN]∈RD×NX=[x_1,\cdots,x_N]\in \mathbb R^{D\times N}X=[x1​,⋯,xN​]∈RD×N(这里使用列向量是为了计算方便),投影矩阵Bm−1:=∑i=1m−1bibi⊤B_{m-1}:=\sum\limits^{m-1}_{i=1}b_ib_i^\topBm−1​:=i=1∑m−1​bi​bi⊤​
所以想要找到第m个主成分,我们需要最大化方差;
Vm=V[zm]=1N∑n=1N(bm⊤x^n)2=bm⊤S^bm,s.t.∥bm∥2=1V_m=\mathbb V[z_m]=\frac1N \sum^N_{n=1}(b^\top_m \hat x_n)^2=b^\top_m \hat Sb_m,\quad s.t. \ \|b_m\|^2=1Vm​=V[zm​]=N1​n=1∑N​(bm⊤​x^n​)2=bm⊤​S^bm​,s.t. ∥bm​∥2=1
其中,S^\hat SS^表示为数据集在正交变换之后(X^\hat\mathcal XX^)的方差.
假设我们已经知道了S^\hat SS^的特征向量,设bib_ibi​为S的特征向量:
S^bi=1NX^X^⊤bi=1N(X−Bm−1X)(X−Bm−1X)⊤bi=(S−SBm−1−Bm−1S+Bm−1SBm−1)bi,(∗)\begin{aligned}\hat Sb_i &= \frac 1N \hat X \hat X^\top b_i=\frac1N(X-B_{m-1}X)(X-B_{m-1}X)^\top b_i\\&=(S-SB_{m-1}-B_{m-1}S+B_{m-1}SB_{m-1})b_i,\end{aligned}\quad (*)S^bi​​=N1​X^X^⊤bi​=N1​(X−Bm−1​X)(X−Bm−1​X)⊤bi​=(S−SBm−1​−Bm−1​S+Bm−1​SBm−1​)bi​,​(∗)
由于bib_ibi​都是这个子空间下的规范正交基(ONB),所以:

Bm−1bi={bi,i<m0,i≥m\boldsymbol B_{m-1}\boldsymbol b_i=\left\{ \begin{aligned} \boldsymbol b_i, \quad i< m \\ \boldsymbol0, \quad i\ge m \end{aligned} \right.\\ Bm−1​bi​={bi​,i<m0,i≥m​
当i<mi<mi<m时,说明bib_ibi​是子空间下的一个正交基,由于是规范正交基,所以与其他基向量的乘积为0,与自身相乘仍为自身。当i≥mi\ge mi≥m时,说明bib_ibi​不是子空间下的正交基,这时候,这bib_ibi​与其他的所有的正交基相互垂直,所以与他们的乘积也就为0.
由上面的关系可以得到:
S^bi=(S−Bm−1S)bi=Sbi=λibiS^bm=Sbm=λmbm\hat S b_i=(S-B_{m-1}S)b_i=Sb_i=\lambda_ib_i\\\hat Sb_m = Sb_m=\lambda_mb_mS^bi​=(S−Bm−1​S)bi​=Sbi​=λi​bi​S^bm​=Sbm​=λm​bm​
这可以知道正交投影之后的向量的特征向量的是一致的。
当i<mi<mi<m时,Bm−1B_{m-1}Bm−1​的关系式带入到(*)中:
S^bi=(S−SBm−1−Bm−1S+Bm−1SBm−1)bi=0=0bi\hat{\boldsymbol{S}} b_{i}=\left(\boldsymbol{S}-\boldsymbol{S} \boldsymbol{B}_{m-1}-\boldsymbol{B}_{m-1} \boldsymbol{S}+\boldsymbol{B}_{m-1} \boldsymbol{S} \boldsymbol{B}_{m-1}\right) \boldsymbol{b}_{i}=\mathbf{0}=0 \boldsymbol{b}_{i} S^bi​=(S−SBm−1​−Bm−1​S+Bm−1​SBm−1​)bi​=0=0bi​
所以可以发现b1,⋯,bm−1b_1,\cdots,b_{m-1}b1​,⋯,bm−1​张成于S^\hat SS^的零空间
由S^bm=Sbm=λmbm\hat Sb_m = Sb_m=\lambda_mb_mS^bm​=Sbm​=λm​bm​和bm⊤bm=1b^\top_mb_m=1bm⊤​bm​=1,可以得到数据在m维上的正交投影的方差为:
Vm=bm⊤Sbm=λmbm⊤bm=λmV_m=b_m^\top S b_m=\lambda_mb^\top_mb_m=\lambda_mVm​=bm⊤​Sbm​=λm​bm⊤​bm​=λm​
由上式可以看到数据方差于对应的特征值之间的关系。

由上表可知,在200 个特征值中,仅有少数的特征值是显著大于0的,所以方差只存在于少数的主成分之中。
为了评估PCA造成的信息损失,我们有以下标准:
M个主成分所能包含的最大方差:
Vm=∑m=1MλmV_m=\sum^M_{m=1}\lambda_mVm​=m=1∑M​λm​
其中的λm\lambda_mλm​是前M个最大的特征值
因数据压缩导致的方差损失:
Jm:=∑j=M+1Dλi=VD−VmJ_m:=\sum^D_{j=M+1}\lambda_i=V_D-V_mJm​:=j=M+1∑D​λi​=VD​−Vm​
或者使用相对方差捕获率(the relative variance captured)VMVD\frac{V_M}{V_D}VD​VM​​,或者是压缩方差损失1−VMVD1-\frac{V_M}{V_D}1−VD​VM​​

投影的角度看待PCA(Projection Perspective)

我们可以将PCA理解为找到一个子空间,源数据在上面的正交投影与源数据最为相似,也就是正交投影的数据与源数据的欧几里得距离最小。

问题设置和问题目标(Setting and Objective)

假设一个规范正交基B=(b1,⋯,bD)∈RDB=(b_1,\cdots,b_D)\in \mathbb R^DB=(b1​,⋯,bD​)∈RD,所以在这个空间中的所有的向量都可以看成是这些正交基的线性组合:
x=∑d=1Dζdbd=∑m=1Mζmbm+∑j=M+1Dζjbj,ζ∈Rx=\sum_{d=1}^D\zeta_db_d=\sum^M_{m=1}\zeta_mb_m+\sum^D_{j=M+1}\zeta_jb_j,\quad \zeta \in \mathbb Rx=d=1∑D​ζd​bd​=m=1∑M​ζm​bm​+j=M+1∑D​ζj​bj​,ζ∈R
在一个低维的子空间中(U⊆RD,dim⁡(U)=MU\subseteq \mathbb R^D, \operatorname{dim}(U)=MU⊆RD,dim(U)=M):
x~=∑m=1Mzmbm∈U∈RD\tilde x = \sum^M_{m=1}z_mb_m\in U\in\mathbb R^Dx~=m=1∑M​zm​bm​∈U∈RD
我们的目标就是最小化两种向量之间的欧几里得距离∥x−x~∥\|x-\tilde x\|∥x−x~∥,这个最小化的向量所在的空间被称为主子空间(Principal Subspace),标记为:
x~n:=∑m=1Mzmnbm=Bzn∈RD,zn:=[z1n,⋯,zMn]⊤∈RM\tilde x_n:=\sum^M_{m=1}z_{mn}b_m=Bz_n\in \mathbb R^D,\quad z_n := [z_{1n},\cdots,z_{Mn}]^\top\in \mathbb R^Mx~n​:=m=1∑M​zmn​bm​=Bzn​∈RD,zn​:=[z1n​,⋯,zMn​]⊤∈RM
znz_nzn​为投影矩阵的坐标。
描述PCA之后的损失的量度为重构误差(Reconstruction Error):
Jm:=1N∑n=1N∥xn−x~n∥2J_m:=\frac 1N \sum^N_{n=1}\| x_n-\tilde x_n\|^2Jm​:=N1​n=1∑N​∥xn​−x~n​∥2

找到最优化坐标(Finding Optimal Coordinates)

想要找到最优化的坐标,需要找到原向量在基向量空间中的正交映射.如下图所示,我的目标也可以理解为找到最小的x~−x\tilde x-xx~−x,由图中可以知道最小的时候是向量正交投影到基向量上的时候。

接下来我们从数学的角度理解这个结论。
对于一个规范正交基b=(b1,⋯,bM),U⊆RDb=(b_1,\cdots,b_M),U \subseteq \mathbb R^Db=(b1​,⋯,bM​),U⊆RD,假设最优的坐标为z1n,⋯,zMnz_{1n},\cdots,z_{Mn}z1n​,⋯,zMn​,对应投影x~n,n=1,⋯,N\tilde x_n,n=1,\cdots,Nx~n​,n=1,⋯,N。为了找到各个维度(坐标)下的最佳的坐标,我们需要将目标函数对坐标进行求导
∂JM∂zin=∂JM∂x~n∂x~n∂zin∂JM∂x~n−2N(xn−x~n)⊤∈R1×D\begin{aligned}&\frac {\partial J_M}{\partial z_{in}}=\frac{\partial J_M}{\partial\tilde x_n}\frac{\partial \tilde x_n}{\partial z_{in}} \\ &\frac{\partial J_M}{\partial \tilde x_n}-\frac{2}{N}(x_n-\tilde x_n)^\top\in \mathbb R^{1\times D}\end {aligned}​∂zin​∂JM​​=∂x~n​∂JM​​∂zin​∂x~n​​∂x~n​∂JM​​−N2​(xn​−x~n​)⊤∈R1×D​
因为:
x~n:=∑m=1Mzmnbm=Bzn∈RD\tilde x_n:=\sum^M_{m=1}z_{mn}b_m=Bz_n\in \mathbb R^Dx~n​:=m=1∑M​zmn​bm​=Bzn​∈RD
所以有:
∂JM∂zin=−2N(xn−x~n)⊤bi=−2N(xn−∑m=1Mzmnbm)⊤bi=bibj=0−2N(xn⊤bi−zinbi⊤bi)=−2N(xn⊤bi−zin)\frac {\partial J_M}{\partial z_{in}}=-\frac 2N(x_n-\tilde x_n)^\top b_i=-\frac 2N (x_n-\sum_{m=1}^Mz_{mn}b_m)^\top b_i\overset{b_ib_j=0}{=}-\frac{2}{N}(x_n^\top b_i-z_{in}b^\top_ib_i)=-\frac2N(x_n^\top b_i-z_{in})∂zin​∂JM​​=−N2​(xn​−x~n​)⊤bi​=−N2​(xn​−m=1∑M​zmn​bm​)⊤bi​=bi​bj​=0−N2​(xn⊤​bi​−zin​bi⊤​bi​)=−N2​(xn⊤​bi​−zin​)
将上面的偏微分设为0,可以得到最优情况下的坐标:
zin=xn⊤bi=bi⊤xn,i=1⋯M,n=1,⋯,Nz_{in}=x_n^\top b_i=b_i^\top x_n,\quad i=1\cdots M,n=1,\cdots ,Nzin​=xn⊤​bi​=bi⊤​xn​,i=1⋯M,n=1,⋯,N
这就说明最优坐标就是将原始数据做正交投影到目标向量空间中的坐标。

现在我们稍微复习一下向基向量的正交投影:
一个向量向正交基(b1,⋯,bD)∈RD(b_1,\cdots,b_D)\in \mathbb R^D(b1​,⋯,bD​)∈RD进行正交投影
x~=bj(bj⊤bj⏟ONB,=I)−1bj⊤x=bjbj⊤x∈RD\tilde x=b_j(\underbrace{ b_j^\top b_j}_{ONB,=I})^{-1}b_j^\top x=b_j b_j^\top x \in \mathbb R^Dx~=bj​(ONB,=Ibj⊤​bj​​​)−1bj⊤​x=bj​bj⊤​x∈RD
其中,bj⊤xb_j^\top xbj⊤​x是正交投影之后的坐标

在新的坐标系下,虽然x~∈RD\tilde x\in \mathbb R^Dx~∈RD,但是我们只需要用前M个坐标,因为在这个坐标系下剩下的坐标都是0.

找到主子空间的基向量(Finding the Basis of the Principal Subspace)

为了找到主子空间的基向量,我们需要对原先的代价函数的形式进行一些改造。:
x~n=∑m=1Mzmnbm=∑m=1M(xn⊤bm)bm\tilde x _n = \sum_{m=1}^Mz_{mn}b_m=\sum _{m=1}^M(x_n^\top b_m)b_m x~n​=m=1∑M​zmn​bm​=m=1∑M​(xn⊤​bm​)bm​
根据点积的对称性:
x~n=(∑m=1Mbmbm⊤)xn\tilde x _n=(\sum^M_{m=1}b_mb_m^\top)x_nx~n​=(m=1∑M​bm​bm⊤​)xn​

补充(原因)

原先提到原始数据可以用基向量线性组合表示,所以(这里可以理解为将原向量分解为投影向量和位移向量)
xn=∑d=1Dzdnbd=(10.32)∑d=1D(xn⊤bd)bd=(∑d=1Dbdbd⊤)xn=(∑m=1Mbmbm⊤)xn+(∑j=M+1Dbjbj⊤)xn\begin{aligned} \boldsymbol{x}_{n} &=\sum_{d=1}^{D} z_{d n} \boldsymbol{b}_{d} \stackrel{(10.32)}{=} \sum_{d=1}^{D}\left(\boldsymbol{x}_{n}^{\top} \boldsymbol{b}_{d}\right) \boldsymbol{b}_{d}=\left(\sum_{d=1}^{D} b_{d} \boldsymbol{b}_{d}^{\top}\right) \boldsymbol{x}_{n} \\ &=\left(\sum_{m=1}^{M} \boldsymbol{b}_{m} \boldsymbol{b}_{m}^{\top}\right) \boldsymbol{x}_{n}+\left(\sum_{j=M+1}^{D} \boldsymbol{b}_{j} \boldsymbol{b}_{j}^{\top}\right) \boldsymbol{x}_{n} \end{aligned} xn​​=d=1∑D​zdn​bd​=(10.32)d=1∑D​(xn⊤​bd​)bd​=(d=1∑D​bd​bd⊤​)xn​=(m=1∑M​bm​bm⊤​)xn​+⎝⎛​j=M+1∑D​bj​bj⊤​⎠⎞​xn​​
所以位移向量(displacement vector)为:
xn−x~n=(∑j=M+1Dbjbj⊤)xn=∑j=M+1D(xn⊤bj)bj\begin{aligned} x_n-\tilde x_n&=(\sum_{j=M+1}^Db_jb_j^\top)x_n\\&=\sum^D_{j=M+1}(x_n^\top b_j)b_j\end{aligned}xn​−x~n​​=(j=M+1∑D​bj​bj⊤​)xn​=j=M+1∑D​(xn⊤​bj​)bj​​
其中,∑j=M+1Dbjbj⊤\sum_{j=M+1}^Db_jb_j^\top∑j=M+1D​bj​bj⊤​为投影矩阵。

这里可以看出,位移矩阵是在垂直于主子空间的空间中。

低秩近似((Low-Rank Approximation):
由之前的讨论得知,投影矩阵为:
∑m=1Mbmbm⊤=BB⊤\sum_{m=1}^Mb_mb^\top_m=BB^\topm=1∑M​bm​bm⊤​=BB⊤
由此,原先的平均平方重构误差可以写为:
1N∑n=1N∥xn−x~∥2=1N∑n=1N∥xn−BB⊤xn∥2=1N∑n=1N∥(I−BB⊤)xn∥2\frac 1N\sum^N_{n=1}\|x_n-\tilde x\|^2=\frac 1N \sum_{n=1}^N\|x_n-BB^\top x_n\|^2=\frac 1N \sum^N_{n=1}\|(I-BB^\top)x_n\|^2N1​n=1∑N​∥xn​−x~∥2=N1​n=1∑N​∥xn​−BB⊤xn​∥2=N1​n=1∑N​∥(I−BB⊤)xn​∥2
所以可以将PCA理解为找到与单位矩阵最接近的BB⊤BB^\topBB⊤的MMM秩逼近。

现在我们能够重构损失函数:
JM=1N∑n=1N∥xn−x~n∥2=1N∑n=1N∥∑j=M+1D(bj⊤xn)bj∥2J_M=\frac 1N\sum^N_{n=1}\|x_n-\tilde x_n\|^2=\frac 1N \sum^N_{n=1}\Vert\sum_{j=M+1}^D(b_j^\top x_n)b_j\|^2JM​=N1​n=1∑N​∥xn​−x~n​∥2=N1​n=1∑N​∥j=M+1∑D​(bj⊤​xn​)bj​∥2
我们将平方范数展开,并结合bjb_jbj​是来源于规范正交基,可以得到下式:
JM=1N∑n=1N∑j=M+1D(bj⊤xn)2=1N∑n=1N∑j=M+1Dbj⊤xnbj⊤xn=1N∑n=1N∑j=M+1Dbj⊤xnxn⊤bjJ_M=\frac 1N \sum^N_{n=1}\sum^D_{j=M+1}(b_j^\top x_n)^2=\frac {1}{N}\sum^N_{n=1}\sum_{j=M+1}^D b_j^\top x_nb^\top_j x_n=\frac 1N \sum^N_{n=1}\sum^D_{j=M+1}b_j^\top x_n x_n^\top b_jJM​=N1​n=1∑N​j=M+1∑D​(bj⊤​xn​)2=N1​n=1∑N​j=M+1∑D​bj⊤​xn​bj⊤​xn​=N1​n=1∑N​j=M+1∑D​bj⊤​xn​xn⊤​bj​

补充推导过程

由于点乘的对称性,我们可以知道bj⊤xn=xn⊤bjb^\top_jx_n=x^\top_nb_jbj⊤​xn​=xn⊤​bj​,带入上式:
JM=∑j=M+1Dbj⊤(1N∑n=1Nxnxn⊤)⏟=:Sbj=∑j=M+1Dbj⊤Sbj=∑j=M+1Dtr⁡(bj⊤Sbj)=∑j=M+1Dtr⁡(Sbjbj⊤)=tr⁡((∑j=M+1Dbjbj⊤)⏟projection matrix S)\begin{aligned} J_{M} &=\sum_{j=M+1}^{D} \boldsymbol{b}_{j}^{\top} \underbrace{\left(\frac{1}{N} \sum_{n=1}^{N} \boldsymbol{x}_{n} \boldsymbol{x}_{n}^{\top}\right)}_{=: \boldsymbol{S}} \boldsymbol{b}_{j}=\sum_{j=M+1}^{D} \boldsymbol{b}_{j}^{\top} \boldsymbol{S} \boldsymbol{b}_{j} \\ &=\sum_{j=M+1}^{D} \operatorname{tr}\left(\boldsymbol{b}_{j}^{\top} \boldsymbol{S} \boldsymbol{b}_{j}\right)=\sum_{j=M+1}^{D} \operatorname{tr}\left(\boldsymbol{S} \boldsymbol{b}_{j} \boldsymbol{b}_{j}^{\top}\right)=\operatorname{tr}(\underbrace{\left(\sum_{j=M+1}^{D} \boldsymbol{b}_{j} \boldsymbol{b}_{j}^{\top}\right)}_{\text {projection matrix }} \boldsymbol{S}) \end{aligned} JM​​=j=M+1∑D​bj⊤​=:S(N1​n=1∑N​xn​xn⊤​)​​bj​=j=M+1∑D​bj⊤​Sbj​=j=M+1∑D​tr(bj⊤​Sbj​)=j=M+1∑D​tr(Sbj​bj⊤​)=tr(projection matrix ⎝⎛​j=M+1∑D​bj​bj⊤​⎠⎞​​​S)​
由上可知,损失函数可以被理解为源数据在主子空间的正交补上的方差。这也对应这主成分分析是在最小化我们忽略的维度上的误差。等价的来说也就是我们需要保留方差最大的那几个维度。所以当我们投影到M维主子空间的时候,所对应的重构误差为:
JM=∑j=M+1DλjJ_M=\sum^D_{j=M+1}\lambda_jJM​=j=M+1∑D​λj​

为什么是这个?

其中的λ\lambdaλ数据协方差的奇异值。所以想要最小化这个重构误差,就需要选择D−MD-MD−M个最小的特征值,这些特征值对应的是主子空间的正交基的特征向量。这也就是说,主子空间所对应的特征向量的特征值是协方差矩阵中的最大的M个特征值。

这一节有很多问题,待补充。。。

特征向量计算以及低秩逼近(Eigenvector Computation and Low-Rank Approximations)

为了计算方差矩阵的特征值,我们可以采用特征值分解或者是奇异值分解,前者可以直接计算出矩阵的特征值和特征向量。而使用SVD的可行性,是因为方差矩阵是对称并且能够分解为XX⊤XX^\topXX⊤所以,方差矩阵的特征值就是XXX的奇异值的平方。
S=1N∑n=1Nxnxn⊤=1NXX⊤,X=[x1,⋯,xN]∈RD×NS=\frac 1N \sum^N_{n=1}x_nx_n^\top = \frac 1N XX^\top,\quad X=[x_1,\cdots , x_N]\in \mathbb R^{D\times N}S=N1​n=1∑N​xn​xn⊤​=N1​XX⊤,X=[x1​,⋯,xN​]∈RD×N
矩阵XXX对应的SVD为:
X⏟D×N=U⏟D×DΣ⏟D×NV⊤⏟N×N\underbrace X_{D\times N}=\underbrace U_{D\times D}\underbrace\Sigma_{D\times N}\underbrace {V^\top}_{N\times N}D×NX​​=D×DU​​D×NΣ​​N×NV⊤​​
其中U和V都是正交矩阵,Σ\SigmaΣ为对角矩阵,主对角线上的元素为奇异值σii≥0\sigma_{ii}\ge 0σii​≥0.将这个式子带入到方差矩阵中:
S=1NXX⊤=1NUΣV⊤V⏟=INΣ⊤U⊤=1NUΣΣ⊤U⊤S=\frac 1N XX^\top=\frac 1NU\Sigma\underbrace{V^\top V}_{=I_N}\Sigma^\top U^\top=\frac 1N U\Sigma\Sigma^\top U^\topS=N1​XX⊤=N1​UΣ=IN​V⊤V​​Σ⊤U⊤=N1​UΣΣ⊤U⊤

SVD分解之后的两端的矩阵是酉矩阵(V⊤=V−1V^\top=V^{-1}V⊤=V−1):
Specifically, the singular value decomposition of an m\times n complex matrix M is a factorization of the form UΣV∗\mathbf {U\Sigma V^{*}}UΣV∗, where U is an m×mm\times mm×m complex unitary matrix, Σ\mathbf{\Sigma}Σ is an m×nm\times nm×n rectangular diagonal matrix with non-negative real numbers on the diagonal, and V is an n×nn\times nn×n complex unitary matrix(酉矩阵).

所以U的列向量是XX⊤XX^\topXX⊤的特征向量,也是方差矩阵的特征向量。其中特征值与奇异值的关系为:
λd=σd2N\lambda_d=\frac{\sigma^2_d}{N}λd​=Nσd2​​
S的特征值和X的奇异值的关系对应的是原先的最大方差视角和奇异值分解之间的关系。

如何理解?

用低秩逼近的PCA(PCA Using Low-Rank Matrix Approximations)

PCA需要找出前N个最大特征值所对应的特征向量,实现这个目标可以采用低秩逼近的方式。

Eckart-Young theorem:就是评估低秩逼近之后造成的损失

由Eckart-Young theorem:
X~M:=argmin⁡rk(A)≤M⁡∥X−A∥2∈RD×N\tilde X_M:=\operatorname{argmin}_{\operatorname{rk(A)\le M}}\|X-A\|_2\in \mathbb R^{D\times N}X~M​:=argminrk(A)≤M​∥X−A∥2​∈RD×N
所以,对应的低秩逼近就是找出前M大的奇异值:
X~M=UM⏟D×MΣM⏟M×MVM⊤⏟M×N∈RD×N\tilde X_M=\underbrace {U_M}_{D\times M}\underbrace{\Sigma_M}_{M\times M}\underbrace{V_M^\top}_{M\times N}\in \mathbb R^{D\times N}X~M​=D×MUM​​​M×MΣM​​​M×NVM⊤​​​∈RD×N
其中,Σ\SigmaΣ包含X的前M个最大的奇异值。

实际方面(Practical Aspects)

我们可以直接采用特征多项式求解出特征值和特征多项式,但是由Abel-Ruffini theorem,五阶或者五阶以上的多项式方程没有几何解。所以在解决大于4×44\times 44×4的矩阵的时候会遇到这个问题。
由于在主成分分析中,我们会只需要前M大的特征向量和特征多项式,所以计算出所有的特征向量和特征值然后再舍弃一些特征值是很没有必要的。在一些极端的情况下,我们只需要第一个特征向量,这时候使用幂迭代(Power iteration)效率会非常高。

幂迭代
首先先随机选取一个不在S的零空间的向量x0x_0x0​,然后按照下式进行迭代:
xk+1=Sxk∥Sxk∥,k=0,1,⋯x_{k+1}=\frac{Sx_k}{\|Sx_k\|},\quad k=0,1,\cdotsxk+1​=∥Sxk​∥Sxk​​,k=0,1,⋯
这个式子总是有∥xk∥=1\|x_k\|=1∥xk​∥=1,最终这个式子会收敛于最大的特征值所对应的特征向量。当S为不可逆的时候,应该保证x0≠0x_0 \ne 0x0​​=0

In mathematics, power iteration (also known as the power method) is an eigenvalue algorithm: given a diagonalizable matrix A, the algorithm will produce a number \lambda , which is the greatest (in absolute value) eigenvalue of A, and a nonzero vector v, which is a corresponding eigenvector of λ\lambdaλ , that is, Av=λvAv=\lambda vAv=λv. The algorithm is also known as the Von Mises iteration.

Power iteration is a very simple algorithm, but it may converge slowly. The most time-consuming operation of the algorithm is the multiplication of matrix A by a vector, so it is effective for a very large sparse matrix with appropriate implementation.

高维PCA(PCA in High Dimensions)

想要对数据使用PCA,需要求解出数据的协方差矩阵,对于一个D维的数据,如果使用特征多项式(∣λE−A∣=0|\lambda E-A|=0∣λE−A∣=0)的时间复杂度为O(D3)O(D^3)O(D3).所以需要找到一种更加高效的方法解决这个问题。
下面我们讨论数据的个数远小于数据维度的情况,即N≪DN\ll DN≪D
假设一组中心化(均值为0)的数据集x1,⋯,xN,xn∈RD×Dx_1,\cdots,x_N,\ \ x_n\in\mathbb R^{D\times D}x1​,⋯,xN​,  xn​∈RD×D,对应的协方差矩阵为:
S=1NXX⊤∈RD×D,X=[x1,⋯,xN]∈RD×NS=\frac 1N XX^\top\in\mathbb R^{D\times D},\quad X=[x_1,\cdots,x_N]\in\mathbb R^{D\times N}S=N1​XX⊤∈RD×D,X=[x1​,⋯,xN​]∈RD×N
由于我们假设N≪DN\ll DN≪D所以数据点的数量远小于数据的维度,也就是说数据的秩为N,则有D−N+1D-N+1D−N+1个特征值为0,接下来我们探究将D维协方差矩阵转换成N维,且对应的特征值都是正数。所以有特征向量的等式:
Sbm=λmbm,m=1,⋯MSb_m=\lambda_m b_m,\quad m=1,\cdots MSbm​=λm​bm​,m=1,⋯M
其中b是主子空间的基向量,现在将S的定义带入:
Sbm=1NXX⊤bm=λmbmSb_m =\frac 1N XX^\top b_m=\lambda_m b_mSbm​=N1​XX⊤bm​=λm​bm​
现在等式两边同时乘以X⊤∈RN×DX^\top\in \mathbb R^{N\times D}X⊤∈RN×D
1NX⊤X⏟N×NX⊤bm⏟=:cm=λmX⊤bm⇔1NX⊤Xcm=λmcm\frac 1N \underbrace {X^\top X}_{N\times N}\underbrace{X^\top b_m}_{=:c_m}=\lambda_m X^\top b_m\Leftrightarrow\frac 1N X^\top Xc_m=\lambda_mc_mN1​N×NX⊤X​​=:cm​X⊤bm​​​=λm​X⊤bm​⇔N1​X⊤Xcm​=λm​cm​
所以可以发现协方差矩阵的特征值为λm\lambda_mλm​对应的特征向量为cmc_mcm​

印证原先提到的:XX⊤XX^\topXX⊤的非零特征值等于X⊤XX^\top XX⊤X的非零特征值

现在我们得到了映射之后的特征值和特征向量,现在我们需要找到源数据的特征值和特征向量。现在对上式两边同时左乘XXX:
1NXX⊤⏟SXcm=λmXcm\underbrace{\frac 1NXX^\top}_SXc_m=\lambda_mXc_mSN1​XX⊤​​Xcm​=λm​Xcm​
这样我们得到了源数据(X是酉矩阵?),这仍旧是S的特征向量。

PCA在实践中的关键步骤(Key Steps of PCA in Practice)

  1. 减去数据均值(Key Steps of PCA in Practice)
  2. 这一步将所有数据减去数据的均值,使得处理后的数据的均值为0,这一步不是必须的,但是减小遇到数值问题的风险。
  3. 规范化(Standardization):
  4. 将数据除以数据的标准偏差σd\sigma_dσd​
  5. 协方差矩阵的特征值分解(Eigendecomposition of the covariance matrix)
  6. 由于协方差是对称的,根据谱定理,我们能够找到特征向量的规范正交基
  7. 投影(Projection)
  8. 将数据点x∗∈RDx_*\in \mathbb R^Dx∗​∈RD投影到主子空间中:
  9. x(d)←x∗(d)−μdσdd=1,⋯,Dx^{(d)}\leftarrow\frac{x_*^{(d)}-\mu_d}{\sigma_d}\quad d = 1,\cdots,Dx(d)←σd​x∗(d)​−μd​​d=1,⋯,D
    10.其中,x∗(d)x^{(d)}_*x∗(d)​代表的是x∗x_*x∗​的第d个成分,所以对应的投影为:
    x~∗=BB⊤x∗\tilde x_*=BB^\top x_*x~∗​=BB⊤x∗​
    对应的坐标为:
    z∗=B⊤x∗z_*=B^\top x_*z∗​=B⊤x∗​
    其中B由数据协方差矩阵最大的几个特征值所对应的特征向量组成。注意PCA返回的是坐标,而不是投影向量。
    要得到原始数据的投影,我们需要将投影之后的数据进行“反规范化”:
    x~∗(d)←x~∗(d)σd+μd,d=1,⋯,D\tilde x^{(d)}_*\leftarrow \tilde x^{(d)}_*\sigma_d+\mu_d,\quad d= 1,\cdots, Dx~∗(d)​←x~∗(d)​σd​+μd​,d=1,⋯,D

MNIST数字:重构(MNIST Digits: Reconstruction)
由下图可知,当主成分为一的时候,图像就是一个可以识别的数字了,随着主成分的增加,图像变得清晰了些

下图中展示了图像信息损失和主成分数量之间的关系:
1N∑n=1N∥xn−x~n∥2=∑i=M+1Dλi\frac 1N \sum^N_{n=1}\|x_n-\tilde x_n\|^2=\sum^D_{i=M+1}\lambda_iN1​n=1∑N​∥xn​−x~n​∥2=i=M+1∑D​λi​

这也印证之前提到的,大多数的信息只存在于少量的主成分之中。

用潜变量看待PCA(Latent Variable Perspective)

原先讨论PCA的时候没有使用概率方面的理论,这样能够帮助我们避开一些由概率论引起的数学上的困难,但是用概率论的能够帮助我们更好地理解PCA,而且在处理带有噪音的数据的时候,概率论中的似然函数提供了分析方式。

Observational Noise. The error between the true value in a system and its observed value due to imprecision in measurement. Also called Measurement Noise.
观测噪音(Observation Noise)实际上就是我们所说的测量误差,由仪器等因素导致的于真实值的偏差。

通过介绍连续的潜变量z∈RMz\in \mathbb R^Mz∈RM, 我们能够将PCA用概率潜变量模型(probabilistic latent-variable model)来描述,这被称为概率主成分分析(probabilistic PCA , PPCA)

生成过程及概率模型(Generative Process and Probabilistic Model)

我们考虑一个线性降维,对于一个连续随机变量z∈RMz\in \mathbb R^Mz∈RM以及一个标准正态先验p(z)=N(0,I)p(z)=\mathcal N(0, I)p(z)=N(0,I), 潜变量以及观测值之间的关系为:
x=Bz+μ+ϵ∈RDx=Bz+\mu+\epsilon\in \mathbb R^Dx=Bz+μ+ϵ∈RD
其中ϵ∼N(0,σ2I)\epsilon \sim \mathcal N(0,\sigma^2I)ϵ∼N(0,σ2I)为高斯观测噪音,而B∈RD×M,μ∈RDB\in \mathbb R^{D\times M},\quad \mu \in \mathbb R^DB∈RD×M,μ∈RD是潜变量到观测变量的线性/仿射映射。所以,潜变量于观测值之间的联系方式为:
p(x∣z,B,μ,σ2)=N(x∣Bz+μ,σ2I)p(x|z,B,\mu,\sigma^2)=\mathcal N(x|Bz+\mu, \sigma^2I)p(x∣z,B,μ,σ2)=N(x∣Bz+μ,σ2I)
整体来说,PPCA的生成过程为:
z∼N(z∣0,I)xn∣zn∼N(x∣Bzn+μ,σ2I)\begin{aligned} z&\sim \mathcal N(z|0,I)\\ x_n|z_n&\sim\mathcal N(x|Bz_n+\mu,\sigma^2I)\end{aligned}zxn​∣zn​​∼N(z∣0,I)∼N(x∣Bzn​+μ,σ2I)​
想要得到获得这些参数,需要一些典型数据,想要得到这样的数据可以使用祖先抽样(Ancestral sampling)


Ancestral sampling 实际上就是通过采样解决条件概率问题。

在这里,先采样得到潜变量zzz,然后再从潜变量中采样得到预测数据。于是,上面的生辰过程可以写成:
p(x,z∣B,μ,σ2)=p(x∣z,B,μ,σ2)p(x)p(x,z|B,\mu,\sigma^2)=p(x|z,B,\mu,\sigma^2)p(x)p(x,z∣B,μ,σ2)=p(x∣z,B,μ,σ2)p(x)
对应的图模型:

Graphical model for probabilistic PCA. The observations xnx_nxn​ explicitly depend on corresponding latent variables zn∼N(0,I)z_n \sim \mathcal N(0,I)zn​∼N(0,I) The model parameters B;μB;\muB;μ and the likelihood parameter σ\sigmaσ are shared across the dataset.

可以将潜变量用于生成新的数据(补充理解

似然以及联合分布(Likelihood and Joint Distribution)

由原先的概率论部分,我们知道可以采用积分将潜变量消掉:
p(x∣B,μ,σ2)=∫p(x∣z,B,μ,σ2)p(z)dz=∫N(x∣Bz+μ,σ2I)N(z∣0,I)dzp(x|B,\mu,\sigma^2)=\int p(x|z,B,\mu,\sigma^2)p(z)dz=\int \mathcal N(x|Bz+\mu,\sigma^2I)\mathcal N(z|0,I)dzp(x∣B,μ,σ2)=∫p(x∣z,B,μ,σ2)p(z)dz=∫N(x∣Bz+μ,σ2I)N(z∣0,I)dz
还是由原先的知识,我们可以知道这个积分的结果是高斯分布,其均值及方差为:
Ex[x]=Ez[Bz+μ]+Eϵ[ϵ]=μV[x]=Vz[Bz+μ]+Vϵ[ϵ]=Vz[Bz]+σ2I=BVz[z]B⊤+σ2I=BB⊤+σ2I\begin{aligned}\mathbb E_x[x]&=\mathbb E_z[Bz+\mu]+\mathbb E_\epsilon[\epsilon]=\mu \\ \mathbb V[x]&=\mathbb V_z[Bz+\mu]+\mathbb V_\epsilon[\epsilon]=\mathbb V_z[Bz]+\sigma^2I\\&=B\mathbb V_z[z]B^\top+\sigma^2I=BB^\top+\sigma^2I\end{aligned}Ex​[x]V[x]​=Ez​[Bz+μ]+Eϵ​[ϵ]=μ=Vz​[Bz+μ]+Vϵ​[ϵ]=Vz​[Bz]+σ2I=BVz​[z]B⊤+σ2I=BB⊤+σ2I​

先前我们不适用条件概率分布的原因是极大似然估计以及极大似然后验估计需要的似然函数可以是数据以及模型参数的函数,但是不能是潜变量的函数,这里是用积分消去潜变量之后才用的。

潜变量以及模型参数的关系(需要补充)

因为潜变量zzz的线性/仿射变换x=Bzx=Bzx=Bz是联合高斯分布,现在已知一些边际概率分布:p(z)=N(z∣0,I);p(x)=N(x∣μ,BB⊤+σ2I)p(z)=\mathcal N(z|0,I);p(x)=\mathcal N(x|\mu,BB^\top +\sigma^2I)p(z)=N(z∣0,I);p(x)=N(x∣μ,BB⊤+σ2I).所以对应交叉协方差(cross-covariance)为:
Cov⁡[x,z]=Cov⁡z[Bz+μ]=BCov⁡z[z,z]=B\operatorname{Cov}[x,z]=\operatorname{Cov}_z[Bz+\mu]=B\operatorname{Cov}_z[z,z]=BCov[x,z]=Covz​[Bz+μ]=BCovz​[z,z]=B
所以潜变量以及观测到的随机变量之间的联合分布为:
p(x,z∣B,μ,σ2)=N([xz]∣[μ0],[BB⊤+σ2IBB⊤I])p\left(\boldsymbol{x}, \boldsymbol{z} \mid \boldsymbol{B}, \boldsymbol{\mu}, \sigma^{2}\right)=\mathcal{N}\left(\left[\begin{array}{l} \boldsymbol{x} \\ \boldsymbol{z} \end{array}\right] \mid\left[\begin{array}{l} \boldsymbol{\mu} \\ \mathbf{0} \end{array}\right],\left[\begin{array}{cc} \boldsymbol{B} \boldsymbol{B}^{\top}+\sigma^{2} \boldsymbol{I} & \boldsymbol{B} \\ \boldsymbol{B}^{\top} & \boldsymbol{I} \end{array}\right]\right) p(x,z∣B,μ,σ2)=N([xz​]∣[μ0​],[BB⊤+σ2IB⊤​BI​])
其中均值向量的长度为D+MD+MD+M,协方差矩阵的大小为(D+M)×(D+M)(D+M)\times (D+M)(D+M)×(D+M)

后验分布(Posterior Distibution)

由前面提到的联合概率分布p(x,z∣B,μ,σ2)p(x,z|B,\mu,\sigma^2)p(x,z∣B,μ,σ2)可以求得后验分布p(z∣x)p(z|x)p(z∣x)(参数求解方式在概率论那一章有提及)
p(z∣x)=N(z∣m,C)m=B⊤(BB⊤+σ2I)−1(x−μ)C=I−B⊤(BB⊤+σ2I)−1B\begin{aligned} p(\boldsymbol{z} \mid \boldsymbol{x}) &=\mathcal{N}(\boldsymbol{z} \mid \boldsymbol{m}, \boldsymbol{C}) \\ \boldsymbol{m} &=\boldsymbol{B}^{\top}\left(\boldsymbol{B} \boldsymbol{B}^{\top}+\sigma^{2} \boldsymbol{I}\right)^{-1}(\boldsymbol{x}-\boldsymbol{\mu}) \\ \boldsymbol{C} &=\boldsymbol{I}-\boldsymbol{B}^{\top}\left(\boldsymbol{B} \boldsymbol{B}^{\top}+\sigma^{2} \boldsymbol{I}\right)^{-1} \boldsymbol{B} \end{aligned} p(z∣x)mC​=N(z∣m,C)=B⊤(BB⊤+σ2I)−1(x−μ)=I−B⊤(BB⊤+σ2I)−1B​
注意后验协方差与数据无关,协方差矩阵C告诉我们(?)嵌入的可信度(?p343)
我们可以利用这个后验分布得到数据对应的潜变量,然后再利用潜变量得到重构向量x~∗∼p(x∣z∗,B,μ,σ2)\tilde x_*\sim p(x|z_*,B,\mu,\sigma^2)x~∗​∼p(x∣z∗​,B,μ,σ2).将这个过程重复多次,我们能够得到潜变量的后验分布以及其暗含的观测数据

拓展阅读

现在我们想想之前做了什么。我们使用两个角度看待PCA,一个是投影的角度(最小化重构误差),一个是最大化方差的角度,除此之外还有其他的角度。我们先将高维数据x∈RDx\in \mathbb R^Dx∈RD用矩阵B⊤B^\topB⊤转换成用低维表示的数据z∈RMz\in \mathbb R^Mz∈RM,其中B由协方差矩阵的最大的特征值所对应的特征向量组成。得到低阶矩阵之后,我们可以利用投影矩阵BB⊤BB^\topBB⊤将数据复原到源数据的维度:x≈x~=Bz=BB⊤x∈RDx\approx\tilde x=Bz=BB^\top x\in\mathbb R^Dx≈x~=Bz=BB⊤x∈RD.
当然我们还将PCA看成一个线性自动编码机(Linear Auto-encoder)

由此可以得到重构误差:
1N∑n=1N∥xn−x~n∥2=1N∑n=1N∥xn−BB⊤xn∥\frac 1N \sum^N_{n=1}\|x_n-\tilde x_n\|^2=\frac 1N\sum^N_{n=1}\|x_n-BB^\top x_n\|N1​n=1∑N​∥xn​−x~n​∥2=N1​n=1∑N​∥xn​−BB⊤xn​∥
如果我们将线性映射转换成非线性映射,我们就会得到非线性自动编码机。当编码器是神经网络的时候,这个被称为认知网络或推理网络(recognition network or inference network),编码器称为生成器(Generator)。
还有一种对PCA的理解涉及到信息论(information theory),就是将编码当成原始数据的压缩版本。当我们将压缩的信息还原,这并不能得到与原始一摸一样的数据,我们称这个压缩过程为有损失的。所以我们的目标就是尽可能将原始数据与压缩数据之间的相关性最大化。这种关系称为交互信息(the mutual information)。
在讨论PPCA的时候,我们默认模型的参数(B,μB,\muB,μ,似然参数σ2\sigma^2σ2)都是已知的,这些参数为:(我们将DDD维数据投影到MMM维子空间中)
μMl=1N∑n=1NxnBML=T(Λ−σ2I)12RσML2=1D−M∑j=M+1Dλj\begin{aligned}\mu_{Ml} &=\frac 1N \sum^N_{n=1}x_n\\ B_{ML}&=T(\Lambda-\sigma^2I)^{\frac 12}R\\ \sigma^2_{ML}&=\frac{1}{D-M}\sum^D_{j=M+1}\lambda_j\end{aligned}μMl​BML​σML2​​=N1​n=1∑N​xn​=T(Λ−σ2I)21​R=D−M1​j=M+1∑D​λj​​
其中,T∈RD×MT\in \mathbb R^{D\times M}T∈RD×M包含协方差矩阵的M个特征向量,Λ=diag⁡(λ1,⋯,λM)∈RM×M\Lambda=\operatorname{diag}(\lambda_1,\cdots,\lambda_M)\in \mathbb R^{M\times M}Λ=diag(λ1​,⋯,λM​)∈RM×M是一个对角矩阵,包含主子空间所对应的特征向量所对应的特征值。R∈RM×MR\in \mathbb R^{M\times M}R∈RM×M是一个随意的正交矩阵。BMLB_{ML}BML​是极大似然的解。σML2\sigma_{ML}^2σML2​是主子空间的正交补上的平均方差,可以认为是正交映射之后造成的损失。
当处理一个无噪音的数据的时候,也就是σ→0\sigma \rightarrow 0σ→0,这时候PPCA与PCA得到的结构是一致的。由于协方差矩阵是对称的,所以可以被正交化,所以存在一个矩阵T包含S的特征向量:
S=TΛT−1S=T\Lambda T^{-1}S=TΛT−1
数据的协方差矩阵就是高斯似然函数(p(x∣B,μ,σ2)p(x|B,\mu,\sigma^2)p(x∣B,μ,σ2))的协方差矩阵,也就是BB⊤+σ2IBB^\top+\sigma^2IBB⊤+σ2I。当σ→0\sigma\rightarrow 0σ→0时,两种PCA的数据方差相等,所以有:
Cov⁡[X]=TΛT−1=BB⊤⇔B=TΛ12R\operatorname{Cov}[\mathcal X]=T\Lambda T^{-1}=BB^\top\Leftrightarrow B=T\Lambda^{\frac 12}RCov[X]=TΛT−1=BB⊤⇔B=TΛ21​R
所以实际上,这些PCA都是在对数据的协方差矩阵进行分解。

接触下来的内容难度较大,理解不够透彻,后续补充

  1. iterative expectation maximization (EM) algorithm
  2. Bayesian PCA
  3. Markov chain Monte Carlo /variational inference.
  4. independent component analysis
  5. blind-source separation
  6. deep auto-encoder
  7. Gaussian process latent-variable model (GP-LVM)

MML ch 10 主成分分析降维(Dimensionality Reduction with Principal Component Analysis)相关推荐

  1. 机器学习与高维信息检索 - Note 4 - 主成分分析及其现代解释(Principal Component Analysis, PCA)及相关实例

    主成分分析及其现代解释 4. 主成分分析及其现代解释 Principal Component Analysis and Its Modern Interpretations 4.1 几何学解释 The ...

  2. 主成分分析(PCA)(principal component analysis)

    本文主要讲PCA的相关数学推导.PCA的数学推导用线性代数的知识就可以完成. 参考deeplearningbook.org一书2.12 Example: Principal Components An ...

  3. 通过主成分分析实现三维模型对齐【Principal Component Analysis】

    The PCA is based on the statistical representation of a random variable. Suppose we have a finite se ...

  4. 机器学习降维算法一:PCA (Principal Component Analysis)

    引言: 机器学习领域中所谓的降维就是指采用某种映射方法,将原高维空间中的数据点映射到低维度的空间中.降维的本质是学习一个映射函数 f : x->y,其中x是原始数据点的表达,目前最多使用向量表达 ...

  5. pca降维python实例_主成分分析(Principal component analysis, PCA)例子–Python | 文艺数学君...

    摘要这一篇是关于PCA的实战, 我们会举一个例子, 看一下PCA具体在实战中是如何来进行的. 同时我们会比较同一个数据下, 使用PCA(主成分分析)和FA(因子分析)得到结果的不同. 简介 这一篇文章 ...

  6. 【机器学习sklearn】主成分分析PCA(Principal Component Analysis)

    主成分分析方法PCA 前言 一.PCA是什么? 二.代码实践 使用MNIST数据集实现sklearn库里的主成分分析方法 不同主成分个数对应的可解释方差分析(Explained Variance) 总 ...

  7. 机器学习实战-65:主成因分析降维算法(Principal Component Analysis)

    机器学习实战-65:主成因分析降维算法(PCA) 深度学习原理与实践(开源图书)-总目录,建议收藏,告别碎片阅读! 机器学习分为监督学习.无监督学习和半监督学习(强化学习).无监督学习最常应用的场景是 ...

  8. 主成分分析(principal component analysis, PCA)公式

    主成分分析(principal component analysis, PCA)公式 主成分分析 摘要 什么是主成分 求解 PCA 的公式 数学证明 程序验证 参考文献 主成分分析 摘要 主成分分析作 ...

  9. 【ML】主成分分析 PCA(Principal Component Analysis)原理 + 实践 (基于sklearn)

    [ML]主成分分析 PCA(Principal Component Analysis)原理 + 实践 (基于sklearn) 原理简介 实践 数据集 数据处理 使用KNN模型进行分类预测(为了和PCA ...

  10. 基于spss的主成分分析法(Principal Component Analysis,PCA)

    主成分分析(Principal Component Analysis,PCA), 将多个变量通过线性变换以选出较少个数重要变量的一种多元统计分析方法. 在实际课题中,为了全面分析问题,往往提出很多与此 ...

最新文章

  1. LINUX下c语言调用math.h库函数的注意事项
  2. 富士康第三季度净利润10.9亿美元 同比下滑8.7%
  3. reverse()反转字符串的正确使用方式
  4. Event handling in Angular
  5. 电脑pin重置_如果忘记了如何重置Windows PIN
  6. 基于Accord.Audio和百度语言识别
  7. Git:查看所有远程分支以及同步远程代码
  8. SpringBoot实现前后端数据交互、json数据交互、Controller接收参数的几种常用方式...
  9. cmd运行sql文件
  10. String类常用方法
  11. ST语言和C语言的区别 STC
  12. C语言中的EOF是什么?
  13. %02x与%2x 之间的区别
  14. 【详解】模型优化技巧之优化器和学习率调整
  15. 声散射 matlab,一种基于声波散射的高强度聚焦超声声场测量方法与流程
  16. 无线充电组别国一队:浙江工业大学
  17. 走进音视频的世界——音视频的基本概念
  18. PHP使用socks5代理发送邮件
  19. 使用idea打包web项目为war
  20. MVC框架的学习总结

热门文章

  1. 以为精通Java 线程池,看到这些误区,还是年轻了
  2. 计算机信息技术结束语,新学期初二年级计算机信息技术课第四节结尾
  3. UVM中starting_phase
  4. 【Fusion】Conic Quadratic Optimization
  5. 分析XBrowser地址栏使用案例
  6. 通过几道CTF题学习yii2框架
  7. groovy 字符串截取最后一个_认识python之字符串的下标和切片(17)
  8. 计算机集成牌照,车牌识别+证件识别嵌入式识别系统集成
  9. android中实现图片圆形效果
  10. 安庆集团-冲刺日志(第五天)