文章目录
- 问题设置（Problem Setting）
- 最大化方差的角度看PCA（Maximum Variance Perspective）
- 投影的角度看待PCA(Projection Perspective)
- 特征向量计算以及低秩逼近（Eigenvector Computation and Low-Rank Approximations）
- 高维PCA（PCA in High Dimensions）
- PCA在实践中的关键步骤
- 用潜变量看待PCA
- 拓展阅读
对于一些高维的数据，分析难度大，而且想要对这些数据进行可视化几乎是不可能的，并且想要存储这些数据的代价也是及其昂贵的，所以我们想要找到一种能够将数据的维度降低的方法。这其中， 主成分分析法（principal component analysis (PCA)）是最常用的方法之一。
问题设置（Problem Setting）
在PCA中，我们希望能够找到一个一个向量的投影向量
x
~
n
\tilde x_n
x~n,与原向量尽可能相近。
对于一个独立均匀分布的数据集
X
=
{
x
1
,
⋯
,
x
N
}
,
x
n
∈
R
D
\mathcal X=\{x_1,\cdots,x_N\}, x_n\in \mathbb R^D
X={x1,⋯,xN},xn∈RD,它的均值为0， 对应的数据方差矩阵为：
S
=
1
N
∑
n
=
1
N
x
n
x
⊤
S=\frac 1N \sum^N_{n=1}x_nx^\top
S=N1n=1∑Nxnx⊤
压缩之后表示为：
z
n
=
B
⊤
x
n
∈
R
M
z_n = B^\top x_n\in \mathbb R^M
zn=B⊤xn∈RM
其中，B为投影矩阵，定义为：
B
:
=
[
b
1
,
⋯
,
b
M
]
∈
R
D
×
M
B := [b_1,\cdots,b_M]\in \mathbb R^{D\times M}
B:=[b1,⋯,bM]∈RD×M
假设
b
i
b_i
bi为正交规范基，则
b
i
⊤
b
j
=
0
,
i
≠
j
;
b
i
⊤
b
i
=
1
b_i^\top b_j=0, i\ne j;b_i^\top b_i=1
bi⊤bj=0,i=j;bi⊤bi=1.我们希望找到一个M维的子空间
U
⊆
R
D
,
dim
(
U
)
=
M
<
D
U\subseteq \mathbb R^D,\operatorname{dim}(U)=M < D
U⊆RD,dim(U)=M<D,向其中的投影的向量与原先的向量最相似，因为压缩造成的损失最小。我们将投影的数据表示为
x
~
n
∈
U
\tilde x_n\in U
x~n∈U,对应的坐标为
z
n
z_n
zn(基向量为
b
1
,
⋯
,
b
M
b_1,\cdots,b_M
b1,⋯,bM)
PCA的目标是最小化平方重构误差（the Squared Reconstruction Error）
∥
x
n
−
x
~
n
∥
2
\| x_n-\tilde x_n\|^2
∥xn−x~n∥2.
从数据压缩的角度来看，我们先是将源数据压缩到一个更低维度的空间中，对应下图的
z
z
z，然后将压缩的信息复原，对应下图中的
x
~
\tilde x
x~;
z
z
z控制着多少信息能够从
x
x
x到
x
~
\tilde x
x~.在PCA中，我们考虑原始数据与低维数据之间的线性关系，所以有以下关系：
z
=
B
⊤
x
;
x
~
=
B
z
z = B^\top x;\tilde x = Bz
z=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 x x~, which lives in the original data space, but has an intrinsic lower-dimensional representation than x x x.
最大化方差的角度看PCA（Maximum Variance Perspective）
在下图中，我们丢弃了数据关于
x
2
x_2
x2的信息，这样做能够达到降维的效果，而且使得数据的损失最小化，是源数据与降维之后的数据尽可能相似。假设忽略
x
1
x_1
x1的信息，则得到的数据就很不相似了，也就是说这个降维操作导致了很多的信息损失。通过观察可以发现，数据在两个维度上的分散程度不一样。当数据在一个维度上越分散，说明这个维度上所包含的信息也就越多，而方差可以表示数据分散程度的大小，所以从方差的角度理解PCA就是找到低维空间中数据方差最大的维度。
我们对数据进行一个均值归一化（Mean Normalization）,因为我们要研究的是方差，而对数据整体的几何运算并不会影响数据的方差
V
z
[
z
]
=
V
x
[
B
⊤
(
x
−
μ
)
]
=
V
x
[
B
⊤
x
−
B
⊤
μ
]
=
V
x
[
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:
E
z
[
z
]
=
E
x
[
B
⊤
x
]
=
B
⊤
E
x
[
x
]
=
0
\mathbb E_z[z]=\mathbb E_x[\boldsymbol B^\top \boldsymbol x]=\boldsymbol B^\top \mathbb E_x[\boldsymbol x]=\boldsymbol 0
Ez[z]=Ex[B⊤x]=B⊤Ex[x]=0
最大方差的方向（Direction with Maximal Variance）
为了找到数据在低维空间中的最大的方差，我们先找到一个向量
b
1
∈
R
D
b_1\in \mathbb R^D
b1∈RD,数据在这个向量上的投影的方差最大,也就是要最大化
z
∈
R
M
z\in \mathbb R^M
z∈RM中第一个坐标
z
1
z_1
z1的方差：
V
1
:
=
V
[
z
1
]
=
1
N
∑
n
=
1
N
z
1
n
2
V_1 := \mathbb V[z_1]=\frac1N\sum^N_{n=1}z^2_{1n}
V1:=V[z1]=N1n=1∑Nz1n2
我们将
z
1
n
z_{1n}
z1n表示为数据(
x
n
∈
R
D
x_n\in \mathbb R^D
xn∈RD)在低维空间（
z
n
∈
R
M
z_n\in \mathbb R^M
zn∈RM）的第一个坐标。
z
n
z_n
zn的第一个成分为：
z
1
n
=
b
1
⊤
x
n
z_{1n}=b_1^\top x_n
z1n=b1⊤xn
这是
x
n
x_n
xn在
b
1
b_1
b1张成的一维子空间中的正交投影，将上面二式联立：
V
1
=
1
N
∑
n
=
1
N
(
b
1
⊤
x
n
)
2
=
1
N
∑
n
=
1
N
b
1
⊤
x
n
x
n
⊤
b
1
=
b
1
⊤
(
1
N
∑
n
=
1
N
x
n
x
n
⊤
)
b
1
=
b
1
⊤
S
b
1
V_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_1
V1=N1n=1∑N(b1⊤xn)2=N1n=1∑Nb1⊤xnxn⊤b1=b1⊤(N1n=1∑Nxnxn⊤)b1=b1⊤Sb1
其中， S为数据协方差矩阵。由上式可知，正交基（
b
i
b_i
bi）会对最终的方差的结果产生影响,所以这里要求这些基向量为规范正交基（
∥
b
1
∥
2
=
1
\|b_1\|^2=1
∥b1∥2=1），这样问题就转换成一个约束问题：
max
b
1
b
1
⊤
S
b
1
,
s
.
t
.
∥
b
1
∥
2
=
1
\max_{b_1}b_1^\top Sb_1,\quad s.t.\ \|b_1\|^2 = 1
b1maxb1⊤Sb1,s.t. ∥b1∥2=1
利用拉格朗日方法：
L
(
b
1
,
λ
)
=
b
1
⊤
S
b
1
+
λ
1
(
1
−
b
1
⊤
b
1
)
\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
∂
b
1
=
2
b
1
⊤
S
−
2
λ
1
b
1
⊤
,
∂
L
∂
λ
1
=
1
−
b
1
⊤
b
1
\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λ1b1⊤,∂λ1∂L=1−b1⊤b1
令偏微分的结果为0：
S
b
1
=
λ
1
b
1
b
1
⊤
b
1
=
1
\begin{aligned}Sb_1&=\lambda_1 b_1 \\b_1^\top b_1 &= 1\end{aligned}
Sb1b1⊤b1=λ1b1=1
由上式可以知道，
λ
1
\lambda_1
λ1是方差S的一个特征值，
b
1
b_1
b1是一个特征向量，利用这个式子，我们可以将问题转化成：
V
1
=
b
1
⊤
S
b
1
=
λ
1
b
1
⊤
b
1
=
λ
1
V_1=b_1^\top Sb_1 = \lambda_1b_1^\top b_1 = \lambda_1
V1=b1⊤Sb1=λ1b1⊤b1=λ1
所以我们需要找到一个特征值最大的特征向量,这样源数据在投影之后的方差最大，这个特征向量称为主成分（Principal Component）我们可以得到投影数据点：
x
~
n
=
b
1
z
1
n
=
b
1
b
1
⊤
x
n
∈
R
D
\tilde x_n=b_1 z_{1n}=b_1b_1^\top x_n\in \mathbb R^D
x~n=b1z1n=b1b1⊤xn∈RD
注意这里的投影点上的数据是高纬度空间中的数据，但是实际上存储的时候只需要用低纬度的空间信息就可以表示了。
M维子空间下的最大方差（M-dimensional Subspace with Maximal Variance）
m
−
1
m-1
m−1个主成分对应的是
S
S
S的
m
−
1
m-1
m−1个特征向量，这些特征向量对应着最大的
m
−
1
m-1
m−1个特征值。由于
S
=
1
N
∑
n
=
1
N
x
n
x
n
⊤
S = \frac1N \sum\limits_{n=1}^Nx_nx_n^\top
S=N1n=1∑Nxnxn⊤,所以S是一个对称矩阵，所以由谱定理可以得知，这些特征向量能够形成
R
D
\mathbb R^D
RD空间下的
m
−
1
m-1
m−1维子空间的正交规范特征基。想要找到这些正交基，可以使用向量减法：
X
~
:
=
X
=
∑
i
=
1
m
−
1
b
i
b
i
⊤
X
=
X
−
B
m
−
1
X
\tilde X := X=\sum^{m-1}_{i=1}b_ib_i^\top X=X-B_{m-1}X
X~:=X=i=1∑m−1bibi⊤X=X−Bm−1X
其中，数据点的列向量
X
=
[
x
1
,
⋯
,
x
N
]
∈
R
D
×
N
X=[x_1,\cdots,x_N]\in \mathbb R^{D\times N}
X=[x1,⋯,xN]∈RD×N(这里使用列向量是为了计算方便),投影矩阵
B
m
−
1
:
=
∑
i
=
1
m
−
1
b
i
b
i
⊤
B_{m-1}:=\sum\limits^{m-1}_{i=1}b_ib_i^\top
Bm−1:=i=1∑m−1bibi⊤
所以想要找到第m个主成分，我们需要最大化方差;
V
m
=
V
[
z
m
]
=
1
N
∑
n
=
1
N
(
b
m
⊤
x
^
n
)
2
=
b
m
⊤
S
^
b
m
,
s
.
t
.
∥
b
m
∥
2
=
1
V_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=1
Vm=V[zm]=N1n=1∑N(bm⊤x^n)2=bm⊤S^bm,s.t. ∥bm∥2=1
其中，
S
^
\hat S
S^表示为数据集在正交变换之后（
X
^
\hat\mathcal X
X^）的方差.
假设我们已经知道了
S
^
\hat S
S^的特征向量，设
b
i
b_i
bi为S的特征向量：
S
^
b
i
=
1
N
X
^
X
^
⊤
b
i
=
1
N
(
X
−
B
m
−
1
X
)
(
X
−
B
m
−
1
X
)
⊤
b
i
=
(
S
−
S
B
m
−
1
−
B
m
−
1
S
+
B
m
−
1
S
B
m
−
1
)
b
i
,
(
∗
)
\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=N1X^X^⊤bi=N1(X−Bm−1X)(X−Bm−1X)⊤bi=(S−SBm−1−Bm−1S+Bm−1SBm−1)bi,(∗)
由于
b
i
b_i
bi都是这个子空间下的规范正交基（ONB），所以：
B
m
−
1
b
i
=
{
b
i
,
i
<
m
0
,
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−1bi={bi,i<m0,i≥m
当
i
<
m
i<m
i<m时，说明
b
i
b_i
bi是子空间下的一个正交基，由于是规范正交基，所以与其他基向量的乘积为0，与自身相乘仍为自身。当
i
≥
m
i\ge m
i≥m时，说明
b
i
b_i
bi不是子空间下的正交基，这时候，这
b
i
b_i
bi与其他的所有的正交基相互垂直，所以与他们的乘积也就为0.
由上面的关系可以得到：
S
^
b
i
=
(
S
−
B
m
−
1
S
)
b
i
=
S
b
i
=
λ
i
b
i
S
^
b
m
=
S
b
m
=
λ
m
b
m
\hat S b_i=(S-B_{m-1}S)b_i=Sb_i=\lambda_ib_i\\\hat Sb_m = Sb_m=\lambda_mb_m
S^bi=(S−Bm−1S)bi=Sbi=λibiS^bm=Sbm=λmbm
这可以知道正交投影之后的向量的特征向量的是一致的。
当
i
<
m
i<m
i<m时，
B
m
−
1
B_{m-1}
Bm−1的关系式带入到（*）中：
S
^
b
i
=
(
S
−
S
B
m
−
1
−
B
m
−
1
S
+
B
m
−
1
S
B
m
−
1
)
b
i
=
0
=
0
b
i
\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−1S+Bm−1SBm−1)bi=0=0bi
所以可以发现
b
1
,
⋯
,
b
m
−
1
b_1,\cdots,b_{m-1}
b1,⋯,bm−1张成于
S
^
\hat S
S^的零空间
由
S
^
b
m
=
S
b
m
=
λ
m
b
m
\hat Sb_m = Sb_m=\lambda_mb_m
S^bm=Sbm=λmbm和
b
m
⊤
b
m
=
1
b^\top_mb_m=1
bm⊤bm=1,可以得到数据在m维上的正交投影的方差为：
V
m
=
b
m
⊤
S
b
m
=
λ
m
b
m
⊤
b
m
=
λ
m
V_m=b_m^\top S b_m=\lambda_mb^\top_mb_m=\lambda_m
Vm=bm⊤Sbm=λmbm⊤bm=λm
由上式可以看到数据方差于对应的特征值之间的关系。
由上表可知，在200 个特征值中，仅有少数的特征值是显著大于0的，所以方差只存在于少数的主成分之中。
为了评估PCA造成的信息损失，我们有以下标准：
M个主成分所能包含的最大方差：
V
m
=
∑
m
=
1
M
λ
m
V_m=\sum^M_{m=1}\lambda_m
Vm=m=1∑Mλm
其中的
λ
m
\lambda_m
λm是前M个最大的特征值
因数据压缩导致的方差损失：
J
m
:
=
∑
j
=
M
+
1
D
λ
i
=
V
D
−
V
m
J_m:=\sum^D_{j=M+1}\lambda_i=V_D-V_m
Jm:=j=M+1∑Dλi=VD−Vm
或者使用相对方差捕获率（the relative variance captured）
V
M
V
D
\frac{V_M}{V_D}
VDVM，或者是压缩方差损失
1
−
V
M
V
D
1-\frac{V_M}{V_D}
1−VDVM
投影的角度看待PCA(Projection Perspective)
我们可以将PCA理解为找到一个子空间，源数据在上面的正交投影与源数据最为相似，也就是正交投影的数据与源数据的欧几里得距离最小。
问题设置和问题目标（Setting and Objective）
假设一个规范正交基
B
=
(
b
1
,
⋯
,
b
D
)
∈
R
D
B=(b_1,\cdots,b_D)\in \mathbb R^D
B=(b1,⋯,bD)∈RD,所以在这个空间中的所有的向量都可以看成是这些正交基的线性组合：
x
=
∑
d
=
1
D
ζ
d
b
d
=
∑
m
=
1
M
ζ
m
b
m
+
∑
j
=
M
+
1
D
ζ
j
b
j
,
ζ
∈
R
x=\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 R
x=d=1∑Dζdbd=m=1∑Mζmbm+j=M+1∑Dζjbj,ζ∈R
在一个低维的子空间中(
U
⊆
R
D
,
dim
(
U
)
=
M
U\subseteq \mathbb R^D, \operatorname{dim}(U)=M
U⊆RD,dim(U)=M):
x
~
=
∑
m
=
1
M
z
m
b
m
∈
U
∈
R
D
\tilde x = \sum^M_{m=1}z_mb_m\in U\in\mathbb R^D
x~=m=1∑Mzmbm∈U∈RD
我们的目标就是最小化两种向量之间的欧几里得距离
∥
x
−
x
~
∥
\|x-\tilde x\|
∥x−x~∥,这个最小化的向量所在的空间被称为主子空间（Principal Subspace）,标记为：
x
~
n
:
=
∑
m
=
1
M
z
m
n
b
m
=
B
z
n
∈
R
D
,
z
n
:
=
[
z
1
n
,
⋯
,
z
M
n
]
⊤
∈
R
M
\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^M
x~n:=m=1∑Mzmnbm=Bzn∈RD,zn:=[z1n,⋯,zMn]⊤∈RM
z
n
z_n
zn为投影矩阵的坐标。
描述PCA之后的损失的量度为重构误差（Reconstruction Error）:
J
m
:
=
1
N
∑
n
=
1
N
∥
x
n
−
x
~
n
∥
2
J_m:=\frac 1N \sum^N_{n=1}\| x_n-\tilde x_n\|^2
Jm:=N1n=1∑N∥xn−x~n∥2
找到最优化坐标（Finding Optimal Coordinates）
想要找到最优化的坐标，需要找到原向量在基向量空间中的正交映射.如下图所示，我的目标也可以理解为找到最小的
x
~
−
x
\tilde x-x
x~−x,由图中可以知道最小的时候是向量正交投影到基向量上的时候。
接下来我们从数学的角度理解这个结论。
对于一个规范正交基
b
=
(
b
1
,
⋯
,
b
M
)
,
U
⊆
R
D
b=(b_1,\cdots,b_M),U \subseteq \mathbb R^D
b=(b1,⋯,bM),U⊆RD,假设最优的坐标为
z
1
n
,
⋯
,
z
M
n
z_{1n},\cdots,z_{Mn}
z1n,⋯,zMn,对应投影
x
~
n
,
n
=
1
,
⋯
,
N
\tilde x_n,n=1,\cdots,N
x~n,n=1,⋯,N。为了找到各个维度（坐标）下的最佳的坐标，我们需要将目标函数对坐标进行求导
∂
J
M
∂
z
i
n
=
∂
J
M
∂
x
~
n
∂
x
~
n
∂
z
i
n
∂
J
M
∂
x
~
n
−
2
N
(
x
n
−
x
~
n
)
⊤
∈
R
1
×
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
=
1
M
z
m
n
b
m
=
B
z
n
∈
R
D
\tilde x_n:=\sum^M_{m=1}z_{mn}b_m=Bz_n\in \mathbb R^D
x~n:=m=1∑Mzmnbm=Bzn∈RD
所以有：
∂
J
M
∂
z
i
n
=
−
2
N
(
x
n
−
x
~
n
)
⊤
b
i
=
−
2
N
(
x
n
−
∑
m
=
1
M
z
m
n
b
m
)
⊤
b
i
=
b
i
b
j
=
0
−
2
N
(
x
n
⊤
b
i
−
z
i
n
b
i
⊤
b
i
)
=
−
2
N
(
x
n
⊤
b
i
−
z
i
n
)
\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∑Mzmnbm)⊤bi=bibj=0−N2(xn⊤bi−zinbi⊤bi)=−N2(xn⊤bi−zin)
将上面的偏微分设为0，可以得到最优情况下的坐标：
z
i
n
=
x
n
⊤
b
i
=
b
i
⊤
x
n
,
i
=
1
⋯
M
,
n
=
1
,
⋯
,
N
z_{in}=x_n^\top b_i=b_i^\top x_n,\quad i=1\cdots M,n=1,\cdots ,N
zin=xn⊤bi=bi⊤xn,i=1⋯M,n=1,⋯,N
这就说明最优坐标就是将原始数据做正交投影到目标向量空间中的坐标。
现在我们稍微复习一下向基向量的正交投影：
一个向量向正交基 ( b 1 , ⋯ , b D ) ∈ R D (b_1,\cdots,b_D)\in \mathbb R^D (b1,⋯,bD)∈RD进行正交投影
x ~ = b j ( b j ⊤ b j ⏟ O N B , = I ) − 1 b j ⊤ x = b j b j ⊤ x ∈ R D \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^D x~=bj(ONB,=I bj⊤bj)−1bj⊤x=bjbj⊤x∈RD
其中， b j ⊤ x b_j^\top x bj⊤x是正交投影之后的坐标
在新的坐标系下，虽然 x ~ ∈ R D \tilde x\in \mathbb R^D x~∈RD，但是我们只需要用前M个坐标，因为在这个坐标系下剩下的坐标都是0.
找到主子空间的基向量（Finding the Basis of the Principal Subspace）
为了找到主子空间的基向量，我们需要对原先的代价函数的形式进行一些改造。：
x
~
n
=
∑
m
=
1
M
z
m
n
b
m
=
∑
m
=
1
M
(
x
n
⊤
b
m
)
b
m
\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∑Mzmnbm=m=1∑M(xn⊤bm)bm
根据点积的对称性：
x
~
n
=
(
∑
m
=
1
M
b
m
b
m
⊤
)
x
n
\tilde x _n=(\sum^M_{m=1}b_mb_m^\top)x_n
x~n=(m=1∑Mbmbm⊤)xn
补充（原因）
原先提到原始数据可以用基向量线性组合表示，所以(这里可以理解为将原向量分解为投影向量和位移向量)
x
n
=
∑
d
=
1
D
z
d
n
b
d
=
(
10.32
)
∑
d
=
1
D
(
x
n
⊤
b
d
)
b
d
=
(
∑
d
=
1
D
b
d
b
d
⊤
)
x
n
=
(
∑
m
=
1
M
b
m
b
m
⊤
)
x
n
+
(
∑
j
=
M
+
1
D
b
j
b
j
⊤
)
x
n
\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∑Dzdnbd=(10.32)d=1∑D(xn⊤bd)bd=(d=1∑Dbdbd⊤)xn=(m=1∑Mbmbm⊤)xn+⎝⎛j=M+1∑Dbjbj⊤⎠⎞xn
所以位移向量（displacement vector）为：
x
n
−
x
~
n
=
(
∑
j
=
M
+
1
D
b
j
b
j
⊤
)
x
n
=
∑
j
=
M
+
1
D
(
x
n
⊤
b
j
)
b
j
\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∑Dbjbj⊤)xn=j=M+1∑D(xn⊤bj)bj
其中，
∑
j
=
M
+
1
D
b
j
b
j
⊤
\sum_{j=M+1}^Db_jb_j^\top
∑j=M+1Dbjbj⊤为投影矩阵。
这里可以看出，位移矩阵是在垂直于主子空间的空间中。
低秩近似（(Low-Rank Approximation）：
由之前的讨论得知，投影矩阵为:
∑ m = 1 M b m b m ⊤ = B B ⊤ \sum_{m=1}^Mb_mb^\top_m=BB^\top m=1∑Mbmbm⊤=BB⊤
由此，原先的平均平方重构误差可以写为：
1 N ∑ n = 1 N ∥ x n − x ~ ∥ 2 = 1 N ∑ n = 1 N ∥ x n − B B ⊤ x n ∥ 2 = 1 N ∑ n = 1 N ∥ ( I − B B ⊤ ) x n ∥ 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\|^2 N1n=1∑N∥xn−x~∥2=N1n=1∑N∥xn−BB⊤xn∥2=N1n=1∑N∥(I−BB⊤)xn∥2
所以可以将PCA理解为找到与单位矩阵最接近的 B B ⊤ BB^\top BB⊤的 M M M秩逼近。
现在我们能够重构损失函数：
J
M
=
1
N
∑
n
=
1
N
∥
x
n
−
x
~
n
∥
2
=
1
N
∑
n
=
1
N
∥
∑
j
=
M
+
1
D
(
b
j
⊤
x
n
)
b
j
∥
2
J_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\|^2
JM=N1n=1∑N∥xn−x~n∥2=N1n=1∑N∥j=M+1∑D(bj⊤xn)bj∥2
我们将平方范数展开，并结合
b
j
b_j
bj是来源于规范正交基，可以得到下式：
J
M
=
1
N
∑
n
=
1
N
∑
j
=
M
+
1
D
(
b
j
⊤
x
n
)
2
=
1
N
∑
n
=
1
N
∑
j
=
M
+
1
D
b
j
⊤
x
n
b
j
⊤
x
n
=
1
N
∑
n
=
1
N
∑
j
=
M
+
1
D
b
j
⊤
x
n
x
n
⊤
b
j
J_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_j
JM=N1n=1∑Nj=M+1∑D(bj⊤xn)2=N1n=1∑Nj=M+1∑Dbj⊤xnbj⊤xn=N1n=1∑Nj=M+1∑Dbj⊤xnxn⊤bj
补充推导过程
由于点乘的对称性，我们可以知道
b
j
⊤
x
n
=
x
n
⊤
b
j
b^\top_jx_n=x^\top_nb_j
bj⊤xn=xn⊤bj,带入上式：
J
M
=
∑
j
=
M
+
1
D
b
j
⊤
(
1
N
∑
n
=
1
N
x
n
x
n
⊤
)
⏟
=
:
S
b
j
=
∑
j
=
M
+
1
D
b
j
⊤
S
b
j
=
∑
j
=
M
+
1
D
tr
(
b
j
⊤
S
b
j
)
=
∑
j
=
M
+
1
D
tr
(
S
b
j
b
j
⊤
)
=
tr
(
(
∑
j
=
M
+
1
D
b
j
b
j
⊤
)
⏟
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∑Dbj⊤=:S
(N1n=1∑Nxnxn⊤)bj=j=M+1∑Dbj⊤Sbj=j=M+1∑Dtr(bj⊤Sbj)=j=M+1∑Dtr(Sbjbj⊤)=tr(projection matrix
⎝⎛j=M+1∑Dbjbj⊤⎠⎞S)
由上可知，损失函数可以被理解为源数据在主子空间的正交补上的方差。这也对应这主成分分析是在最小化我们忽略的维度上的误差。等价的来说也就是我们需要保留方差最大的那几个维度。所以当我们投影到M维主子空间的时候，所对应的重构误差为：
J
M
=
∑
j
=
M
+
1
D
λ
j
J_M=\sum^D_{j=M+1}\lambda_j
JM=j=M+1∑Dλj
为什么是这个？
其中的 λ \lambda λ数据协方差的奇异值。所以想要最小化这个重构误差，就需要选择 D − M D-M D−M个最小的特征值，这些特征值对应的是主子空间的正交基的特征向量。这也就是说，主子空间所对应的特征向量的特征值是协方差矩阵中的最大的M个特征值。
这一节有很多问题，待补充。。。
特征向量计算以及低秩逼近（Eigenvector Computation and Low-Rank Approximations）
为了计算方差矩阵的特征值，我们可以采用特征值分解或者是奇异值分解，前者可以直接计算出矩阵的特征值和特征向量。而使用SVD的可行性，是因为方差矩阵是对称并且能够分解为
X
X
⊤
XX^\top
XX⊤所以，方差矩阵的特征值就是
X
X
X的奇异值的平方。
S
=
1
N
∑
n
=
1
N
x
n
x
n
⊤
=
1
N
X
X
⊤
,
X
=
[
x
1
,
⋯
,
x
N
]
∈
R
D
×
N
S=\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=N1n=1∑Nxnxn⊤=N1XX⊤,X=[x1,⋯,xN]∈RD×N
矩阵
X
X
X对应的SVD为：
X
⏟
D
×
N
=
U
⏟
D
×
D
Σ
⏟
D
×
N
V
⊤
⏟
N
×
N
\underbrace X_{D\times N}=\underbrace U_{D\times D}\underbrace\Sigma_{D\times N}\underbrace {V^\top}_{N\times N}
D×N
X=D×D
UD×N
ΣN×N
V⊤
其中U和V都是正交矩阵，
Σ
\Sigma
Σ为对角矩阵，主对角线上的元素为奇异值
σ
i
i
≥
0
\sigma_{ii}\ge 0
σii≥0.将这个式子带入到方差矩阵中：
S
=
1
N
X
X
⊤
=
1
N
U
Σ
V
⊤
V
⏟
=
I
N
Σ
⊤
U
⊤
=
1
N
U
Σ
Σ
⊤
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^\top
S=N1XX⊤=N1UΣ=IN
V⊤VΣ⊤U⊤=N1UΣΣ⊤U⊤
SVD分解之后的两端的矩阵是酉矩阵（ V ⊤ = V − 1 V^\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 × m m\times m m×m complex unitary matrix, Σ \mathbf{\Sigma} Σ is an m × n m\times n m×n rectangular diagonal matrix with non-negative real numbers on the diagonal, and V is an n × n n\times n n×n complex unitary matrix(酉矩阵).
所以U的列向量是
X
X
⊤
XX^\top
XX⊤的特征向量，也是方差矩阵的特征向量。其中特征值与奇异值的关系为：
λ
d
=
σ
d
2
N
\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
∈
R
D
×
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
=
U
M
⏟
D
×
M
Σ
M
⏟
M
×
M
V
M
⊤
⏟
M
×
N
∈
R
D
×
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×M
UMM×M
ΣMM×N
VM⊤∈RD×N
其中，
Σ
\Sigma
Σ包含X的前M个最大的奇异值。
实际方面（Practical Aspects）
我们可以直接采用特征多项式求解出特征值和特征多项式，但是由Abel-Ruffini theorem，五阶或者五阶以上的多项式方程没有几何解。所以在解决大于
4
×
4
4\times 4
4×4的矩阵的时候会遇到这个问题。
由于在主成分分析中，我们会只需要前M大的特征向量和特征多项式，所以计算出所有的特征向量和特征值然后再舍弃一些特征值是很没有必要的。在一些极端的情况下，我们只需要第一个特征向量，这时候使用幂迭代（Power iteration）效率会非常高。
幂迭代
首先先随机选取一个不在S的零空间的向量 x 0 x_0 x0,然后按照下式进行迭代：
x k + 1 = S x k ∥ S x k ∥ , k = 0 , 1 , ⋯ x_{k+1}=\frac{Sx_k}{\|Sx_k\|},\quad k=0,1,\cdots xk+1=∥Sxk∥Sxk,k=0,1,⋯
这个式子总是有 ∥ x k ∥ = 1 \|x_k\|=1 ∥xk∥=1,最终这个式子会收敛于最大的特征值所对应的特征向量。当S为不可逆的时候，应该保证 x 0 ≠ 0 x_0 \ne 0 x0=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, A v = λ v Av=\lambda v Av=λ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
(
D
3
)
O(D^3)
O(D3).所以需要找到一种更加高效的方法解决这个问题。
下面我们讨论数据的个数远小于数据维度的情况，即
N
≪
D
N\ll D
N≪D
假设一组中心化（均值为0）的数据集
x
1
,
⋯
,
x
N
,
x
n
∈
R
D
×
D
x_1,\cdots,x_N,\ \ x_n\in\mathbb R^{D\times D}
x1,⋯,xN, xn∈RD×D，对应的协方差矩阵为：
S
=
1
N
X
X
⊤
∈
R
D
×
D
,
X
=
[
x
1
,
⋯
,
x
N
]
∈
R
D
×
N
S=\frac 1N XX^\top\in\mathbb R^{D\times D},\quad X=[x_1,\cdots,x_N]\in\mathbb R^{D\times N}
S=N1XX⊤∈RD×D,X=[x1,⋯,xN]∈RD×N
由于我们假设
N
≪
D
N\ll D
N≪D所以数据点的数量远小于数据的维度，也就是说数据的秩为N，则有
D
−
N
+
1
D-N+1
D−N+1个特征值为0，接下来我们探究将D维协方差矩阵转换成N维，且对应的特征值都是正数。所以有特征向量的等式：
S
b
m
=
λ
m
b
m
,
m
=
1
,
⋯
M
Sb_m=\lambda_m b_m,\quad m=1,\cdots M
Sbm=λmbm,m=1,⋯M
其中b是主子空间的基向量，现在将S的定义带入：
S
b
m
=
1
N
X
X
⊤
b
m
=
λ
m
b
m
Sb_m =\frac 1N XX^\top b_m=\lambda_m b_m
Sbm=N1XX⊤bm=λmbm
现在等式两边同时乘以
X
⊤
∈
R
N
×
D
X^\top\in \mathbb R^{N\times D}
X⊤∈RN×D
1
N
X
⊤
X
⏟
N
×
N
X
⊤
b
m
⏟
=
:
c
m
=
λ
m
X
⊤
b
m
⇔
1
N
X
⊤
X
c
m
=
λ
m
c
m
\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_m
N1N×N
X⊤X=:cm
X⊤bm=λmX⊤bm⇔N1X⊤Xcm=λmcm
所以可以发现协方差矩阵的特征值为
λ
m
\lambda_m
λm对应的特征向量为
c
m
c_m
cm
印证原先提到的： X X ⊤ XX^\top XX⊤的非零特征值等于 X ⊤ X X^\top X X⊤X的非零特征值
现在我们得到了映射之后的特征值和特征向量，现在我们需要找到源数据的特征值和特征向量。现在对上式两边同时左乘
X
X
X：
1
N
X
X
⊤
⏟
S
X
c
m
=
λ
m
X
c
m
\underbrace{\frac 1NXX^\top}_SXc_m=\lambda_mXc_m
S
N1XX⊤Xcm=λmXcm
这样我们得到了源数据(X是酉矩阵？)，这仍旧是S的特征向量。
PCA在实践中的关键步骤
- 减去数据均值（Key Steps of PCA in Practice）
- 这一步将所有数据减去数据的均值，使得处理后的数据的均值为0，这一步不是必须的，但是减小遇到数值问题的风险。
- 规范化（Standardization）：
- 将数据除以数据的标准偏差 σ d \sigma_d σd
- 协方差矩阵的特征值分解（Eigendecomposition of the covariance matrix）
- 由于协方差是对称的，根据谱定理，我们能够找到特征向量的规范正交基
- 投影（Projection）
- 将数据点 x ∗ ∈ R D x_*\in \mathbb R^D x∗∈RD投影到主子空间中：
-
x
(
d
)
←
x
∗
(
d
)
−
μ
d
σ
d
d
=
1
,
⋯
,
D
x^{(d)}\leftarrow\frac{x_*^{(d)}-\mu_d}{\sigma_d}\quad d = 1,\cdots,D
x(d)←σdx∗(d)−μdd=1,⋯,D
10.其中， x ∗ ( d ) x^{(d)}_* x∗(d)代表的是 x ∗ x_* x∗的第d个成分，所以对应的投影为：
x ~ ∗ = B B ⊤ 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, D x~∗(d)←x~∗(d)σd+μd,d=1,⋯,D
MNIST数字：重构（MNIST Digits: Reconstruction）
由下图可知，当主成分为一的时候，图像就是一个可以识别的数字了，随着主成分的增加，图像变得清晰了些
下图中展示了图像信息损失和主成分数量之间的关系：
1 N ∑ n = 1 N ∥ x n − x ~ n ∥ 2 = ∑ i = M + 1 D λ i \frac 1N \sum^N_{n=1}\|x_n-\tilde x_n\|^2=\sum^D_{i=M+1}\lambda_i N1n=1∑N∥xn−x~n∥2=i=M+1∑Dλi
这也印证之前提到的，大多数的信息只存在于少量的主成分之中。
用潜变量看待PCA
原先讨论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 ∈ R M z\in \mathbb R^M z∈RM, 我们能够将PCA用概率潜变量模型(probabilistic latent-variable model)来描述，这被称为概率主成分分析（probabilistic PCA ， PPCA）
生成过程及概率模型（Generative Process and Probabilistic Model）
我们考虑一个线性降维，对于一个连续随机变量
z
∈
R
M
z\in \mathbb R^M
z∈RM以及一个标准正态先验
p
(
z
)
=
N
(
0
,
I
)
p(z)=\mathcal N(0, I)
p(z)=N(0,I), 潜变量以及观测值之间的关系为：
x
=
B
z
+
μ
+
ϵ
∈
R
D
x=Bz+\mu+\epsilon\in \mathbb R^D
x=Bz+μ+ϵ∈RD
其中
ϵ
∼
N
(
0
,
σ
2
I
)
\epsilon \sim \mathcal N(0,\sigma^2I)
ϵ∼N(0,σ2I)为高斯观测噪音，而
B
∈
R
D
×
M
,
μ
∈
R
D
B\in \mathbb R^{D\times M},\quad \mu \in \mathbb R^D
B∈RD×M,μ∈RD是潜变量到观测变量的线性/仿射映射。所以，潜变量于观测值之间的联系方式为：
p
(
x
∣
z
,
B
,
μ
,
σ
2
)
=
N
(
x
∣
B
z
+
μ
,
σ
2
I
)
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
)
x
n
∣
z
n
∼
N
(
x
∣
B
z
n
+
μ
,
σ
2
I
)
\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 实际上就是通过采样解决条件概率问题。
在这里，先采样得到潜变量
z
z
z,然后再从潜变量中采样得到预测数据。于是，上面的生辰过程可以写成:
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 x n x_n xn explicitly depend on corresponding latent variables z n ∼ N ( 0 , I ) z_n \sim \mathcal N(0,I) zn∼N(0,I) The model parameters B ; μ B;\mu B;μ 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
)
d
z
=
∫
N
(
x
∣
B
z
+
μ
,
σ
2
I
)
N
(
z
∣
0
,
I
)
d
z
p(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)dz
p(x∣B,μ,σ2)=∫p(x∣z,B,μ,σ2)p(z)dz=∫N(x∣Bz+μ,σ2I)N(z∣0,I)dz
还是由原先的知识，我们可以知道这个积分的结果是高斯分布，其均值及方差为：
E
x
[
x
]
=
E
z
[
B
z
+
μ
]
+
E
ϵ
[
ϵ
]
=
μ
V
[
x
]
=
V
z
[
B
z
+
μ
]
+
V
ϵ
[
ϵ
]
=
V
z
[
B
z
]
+
σ
2
I
=
B
V
z
[
z
]
B
⊤
+
σ
2
I
=
B
B
⊤
+
σ
2
I
\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
先前我们不适用条件概率分布的原因是极大似然估计以及极大似然后验估计需要的似然函数可以是数据以及模型参数的函数，但是不能是潜变量的函数，这里是用积分消去潜变量之后才用的。
潜变量以及模型参数的关系（需要补充）
因为潜变量
z
z
z的线性/仿射变换
x
=
B
z
x=Bz
x=Bz是联合高斯分布，现在已知一些边际概率分布:
p
(
z
)
=
N
(
z
∣
0
,
I
)
;
p
(
x
)
=
N
(
x
∣
μ
,
B
B
⊤
+
σ
2
I
)
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
[
B
z
+
μ
]
=
B
Cov
z
[
z
,
z
]
=
B
\operatorname{Cov}[x,z]=\operatorname{Cov}_z[Bz+\mu]=B\operatorname{Cov}_z[z,z]=B
Cov[x,z]=Covz[Bz+μ]=BCovz[z,z]=B
所以潜变量以及观测到的随机变量之间的联合分布为：
p
(
x
,
z
∣
B
,
μ
,
σ
2
)
=
N
(
[
x
z
]
∣
[
μ
0
]
,
[
B
B
⊤
+
σ
2
I
B
B
⊤
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
+
M
D+M
D+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
⊤
(
B
B
⊤
+
σ
2
I
)
−
1
(
x
−
μ
)
C
=
I
−
B
⊤
(
B
B
⊤
+
σ
2
I
)
−
1
B
\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
∈
R
D
x\in \mathbb R^D
x∈RD用矩阵
B
⊤
B^\top
B⊤转换成用低维表示的数据
z
∈
R
M
z\in \mathbb R^M
z∈RM,其中B由协方差矩阵的最大的特征值所对应的特征向量组成。得到低阶矩阵之后，我们可以利用投影矩阵
B
B
⊤
BB^\top
BB⊤将数据复原到源数据的维度:
x
≈
x
~
=
B
z
=
B
B
⊤
x
∈
R
D
x\approx\tilde x=Bz=BB^\top x\in\mathbb R^D
x≈x~=Bz=BB⊤x∈RD.
当然我们还将PCA看成一个线性自动编码机(Linear Auto-encoder)
由此可以得到重构误差：
1
N
∑
n
=
1
N
∥
x
n
−
x
~
n
∥
2
=
1
N
∑
n
=
1
N
∥
x
n
−
B
B
⊤
x
n
∥
\frac 1N \sum^N_{n=1}\|x_n-\tilde x_n\|^2=\frac 1N\sum^N_{n=1}\|x_n-BB^\top x_n\|
N1n=1∑N∥xn−x~n∥2=N1n=1∑N∥xn−BB⊤xn∥
如果我们将线性映射转换成非线性映射，我们就会得到非线性自动编码机。当编码器是神经网络的时候，这个被称为认知网络或推理网络（recognition network or inference network），编码器称为生成器(Generator)。
还有一种对PCA的理解涉及到信息论（information theory），就是将编码当成原始数据的压缩版本。当我们将压缩的信息还原，这并不能得到与原始一摸一样的数据，我们称这个压缩过程为有损失的。所以我们的目标就是尽可能将原始数据与压缩数据之间的相关性最大化。这种关系称为交互信息（the mutual information）。
在讨论PPCA的时候，我们默认模型的参数（
B
,
μ
B,\mu
B,μ,似然参数
σ
2
\sigma^2
σ2）都是已知的，这些参数为：（我们将
D
D
D维数据投影到
M
M
M维子空间中）
μ
M
l
=
1
N
∑
n
=
1
N
x
n
B
M
L
=
T
(
Λ
−
σ
2
I
)
1
2
R
σ
M
L
2
=
1
D
−
M
∑
j
=
M
+
1
D
λ
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}
μMlBMLσML2=N1n=1∑Nxn=T(Λ−σ2I)21R=D−M1j=M+1∑Dλj
其中，
T
∈
R
D
×
M
T\in \mathbb R^{D\times M}
T∈RD×M包含协方差矩阵的M个特征向量，
Λ
=
diag
(
λ
1
,
⋯
,
λ
M
)
∈
R
M
×
M
\Lambda=\operatorname{diag}(\lambda_1,\cdots,\lambda_M)\in \mathbb R^{M\times M}
Λ=diag(λ1,⋯,λM)∈RM×M是一个对角矩阵，包含主子空间所对应的特征向量所对应的特征值。
R
∈
R
M
×
M
R\in \mathbb R^{M\times M}
R∈RM×M是一个随意的正交矩阵。
B
M
L
B_{ML}
BML是极大似然的解。
σ
M
L
2
\sigma_{ML}^2
σML2是主子空间的正交补上的平均方差，可以认为是正交映射之后造成的损失。
当处理一个无噪音的数据的时候，也就是
σ
→
0
\sigma \rightarrow 0
σ→0,这时候PPCA与PCA得到的结构是一致的。由于协方差矩阵是对称的，所以可以被正交化，所以存在一个矩阵T包含S的特征向量：
S
=
T
Λ
T
−
1
S=T\Lambda T^{-1}
S=TΛT−1
数据的协方差矩阵就是高斯似然函数（
p
(
x
∣
B
,
μ
,
σ
2
)
p(x|B,\mu,\sigma^2)
p(x∣B,μ,σ2)）的协方差矩阵,也就是
B
B
⊤
+
σ
2
I
BB^\top+\sigma^2I
BB⊤+σ2I。当
σ
→
0
\sigma\rightarrow 0
σ→0时，两种PCA的数据方差相等，所以有：
Cov
[
X
]
=
T
Λ
T
−
1
=
B
B
⊤
⇔
B
=
T
Λ
1
2
R
\operatorname{Cov}[\mathcal X]=T\Lambda T^{-1}=BB^\top\Leftrightarrow B=T\Lambda^{\frac 12}R
Cov[X]=TΛT−1=BB⊤⇔B=TΛ21R
所以实际上，这些PCA都是在对数据的协方差矩阵进行分解。
接触下来的内容难度较大，理解不够透彻，后续补充
- iterative expectation maximization (EM) algorithm
- Bayesian PCA
- Markov chain Monte Carlo /variational inference.
- independent component analysis
- blind-source separation
- deep auto-encoder
- Gaussian process latent-variable model (GP-LVM)