主成分分析(PCA)与K-L变换

\qquad 主成分分析 $\text{(Principal Component Analysis, PCA)}$ 基于输入数据的协方差矩阵，可实现“去除数据集的相关性” —— 在数据处理中常用于通过“高维数据的降维”来实现特征提取，在信号处理中被称为 $\text{K-L}$ 变换，常用于实现“数据压缩”。

这里的“相关性”是指数据集各维度之间的“线性”依赖关系，采用“协方差”来描述
以二维数据 $\mathbf{x}=(x_1,x_2{)}^T$ 为例，“观测数据集不相关”是指协方差 $cov(x_1,x_2)=0$，即：$x_1$ 和 $x_2$ 两个分量不是线性相关的。此时，协方差矩阵为对角阵
详细内容可参考《协方差矩阵与二维高斯分布》

\qquad

1. 主成分分析 —— 基于最大方差的描述

\qquad 考虑观测数据集 { x i } i = 1 K ∈ R N \{\boldsymbol x_i\}_{i=1}^K\in R^N {xi​}i=1K​∈RN，将每个观测数据投影到某个单位方向 $\mathbf{u}\in R^N,{\mathbf{u}}^T\mathbf{u}=1$。从数据处理的角度来看，如果观测数据集在单位方向 $\mathbf{u}$ 上的投影值具有最大方差，那么方向 $\mathbf{u}$ 就为主成分 $\text{(Principal Component)}$。

1.1 投影数据的方差

\qquad 如图 $1$ 所示，（二维）数据集 { x i } i = 1 K \{\boldsymbol x_i\}_{i=1}^K {xi​}i=1K​ （红色点）在投影之后，变成了单位方向 $\mathbf{u}$ 上的向量集 { x ~ i } i = 1 K \{\tilde{\boldsymbol x}_i\}_{i=1}^K {x~i​}i=1K​（绿色点）。

\qquad

图1 取自于《PRML》Fig 12.2
　　　 观测数据点 ${\mathbf{x}}_i$ 经过投影之后为 ${\tilde{\mathbf{x}}}_i=({\mathbf{u}}^T{\mathbf{x}}_i)\mathbf{u}$，即方向为 $\mathbf{u}$，长度为 ${\mathbf{u}}^T{\mathbf{x}}_i$

\qquad 投影后数据 { x ~ i } i = 1 K \{\tilde{\boldsymbol x}_i\}_{i=1}^K {x~i​}i=1K​ （在 $\mathbf{u}$ 方向上）的方差可定义为：

1 K ∑ i = 1 K { u T x i − u T x ˉ } 2 = u T S u \qquad\qquad\qquad\dfrac{1}{K}\displaystyle\sum_{i=1}^K\{\boldsymbol u^T\boldsymbol x_i-\boldsymbol u^T\bar{\boldsymbol x}\}^2=\boldsymbol u^TS\boldsymbol u K1​i=1∑K​{uTxi​−uTxˉ}2=uTSu

\qquad\qquad 其中，$\bar{\mathbf{x}}=\frac{1}{K}\sum _{i=1}^K{\mathbf{x}}_i$　为观测数据集的样本均值

\qquad\qquad 　　　$S=\frac{1}{K}\sum _{i=1}^K({\mathbf{x}}_i-\bar{\mathbf{x}})({\mathbf{x}}_i-\bar{\mathbf{x}}{)}^T$　为观测数据集的协方差矩阵

对于观测数据集 { x 1 , x 2 , ⋯   , x K } ,   ∀ x i ∈ R N \{\boldsymbol x_1,\boldsymbol x_2,\cdots,\boldsymbol x_K\},\ \forall \boldsymbol x_i\in R^N {x1​,x2​,⋯,xK​}, ∀xi​∈RN，记 $\bar{\mathbf{x}}$ 为数据集的样本均值
若将“去掉均值的数据集”写成矩阵形式，记为 $\mathbf{X}=[{\mathbf{x}}_1-\bar{\mathbf{x}},{\mathbf{x}}_2-\bar{\mathbf{x}},\cdots ,{\mathbf{x}}_K-\bar{\mathbf{x}}{]}_{N\times K}$
那么，数据集的（有偏）协方差矩阵为： $S=\frac{1}{K}\mathbf{X}{\mathbf{X}}^T$，【无偏样本协方差矩阵为： $S=\frac{1}{K-1}\mathbf{X}{\mathbf{X}}^T$】
可参考《PCA图像压缩的matlab实现》

\qquad
\qquad 为了求使得方差 ${\mathbf{u}}^{T}S\mathbf{u}$ 最大的 $\mathbf{u}$ 方向，采用以下步骤：

$(1)$ 采用拉格朗日乘子法构造最优化问题：

max ⁡   {   u T S u + λ ( 1 − u T u )   } \qquad\qquad\qquad\qquad\max\ \{\ \boldsymbol u^TS\boldsymbol u+\lambda(1-\boldsymbol u^T\boldsymbol u)\ \} max { uTSu+λ(1−uTu) }

$(2)$ 令其对 $\mathbf{u}$ 的偏导为 $0$，可得到：

$S\mathbf{u}=\lambda \mathbf{u}$

\qquad 　　显然，$\mathbf{u}$ 是协方差矩阵 $S$ 的特征向量
\qquad 　　　　　$\lambda$ 是 $\mathbf{u}$ 所对应的特征值

$(3)$ 将等式两端左乘 ${\mathbf{u}}^T$，可得到：

${\mathbf{u}}^{T}S\mathbf{u}={\mathbf{u}}^T\lambda \mathbf{u}=\lambda {\mathbf{u}}^T\mathbf{u}=\lambda$

\qquad 　　这就说明，观测数据集在 $S$ 的特征向量 $\mathbf{u}$ 方向上投影后，“投影数据的方差值”正好等于对应的特征值 $\lambda$。

\qquad 如图 $2$ 所示，若将协方差矩阵 $S$ 的特征值按从大到小的顺序排列，最大特征值 ${\lambda }_1$ 所对应特征向量 ${\mathbf{u}}_1$ 的方向，即为第一主成分，该观测数据集“最主要的特征”体现在 ${\mathbf{u}}_1$ 方向上，投影后的数据在该方向上分布的散布程度最大；${\lambda }_2$ 所对应特征向量 ${\mathbf{u}}_2$ 的方向，为第二主成分，“次主要的特征”体现在 ${\mathbf{u}}_2$ 方向上；$\cdots \cdots$。
\qquad

图2 从一个均值为 $0$，协方差矩阵为 $\left[\begin{array}{cc}2& -2\\ -2& 4\end{array}\right]$ 的正态分布采样了 $100$ 个数据（蓝色的 ‘+’ 号）组成观测数据集 { x i } i = 1 100 ∈ R 2 \{\boldsymbol x_i\}_{i=1}^{100}\in R^2 {xi​}i=1100​∈R2
　 　
通过对协方差矩阵进行特征分解，可求出：
$(1)$ 第一主成分为 ${\mathbf{u}}_1=\left[\begin{array}{c}0.5257\\ -0.8507\end{array}\right]$（黑色箭头），观测数据集在该方向上的投影具有更大的方差
　　对应了协方差矩阵的最大特征值 ${\lambda }_1=5.2360$
$(2)$ 第二主成分为 ${\mathbf{u}}_2=\left[\begin{array}{c}-0.8507\\ -0.5257\end{array}\right]$（红色箭头），观测数据集在该方向上的投影具有更小的方差
　　对应了协方差矩阵较小的特征值 ${\lambda }_2=0.7639$
　
由 $100$ 个观测点计算得到的“观测数据集的协方差矩阵”$S=\left[\begin{array}{cc}1.7809& -1.7640\\ -1.7640& 3.6450\end{array}\right]$，特征分解之后的两个主成分分别为：
$(1)$ 第一主成分为 ${\mathbf{u}}_1^{{}^{\prime }}=\left[\begin{array}{c}0.5162\\ -0.8565\end{array}\right]$，较大的特征值 ${\lambda }_1^{{}^{\prime }}=4.7081$（右图为投影数据）
　　可以看出，观测数据的主要变化特征（$y$ 随 $x$ 的线性变化）是沿着 ${\mathbf{u}}_1^{{}^{\prime }}$ 方向（可参考《协方差矩阵与二维高斯分布》）
$(2)$ 第二主成分为 ${\mathbf{u}}_2^{{}^{\prime }}=\left[\begin{array}{c}-0.8565\\ -0.5162\end{array}\right]$，较小的特征值 ${\lambda }_2^{{}^{\prime }}=0.7178$
　
可以看出，即使只采样了 $100$ 个数据，两组主成分的方向也是非常接近的

\qquad

1.2 高维数据的降维

\qquad 考虑 $N$ 维观测数据集 { x i } i = 1 K \{\boldsymbol x_i\}_{i=1}^K {xi​}i=1K​，当 $N$ 值比较大时，会极大增加运算量。例如，一幅 $100\times 100$ 大小的图像实际上是 $R^{10000}$ 空间中的一个元素（$N=10000$）。

\qquad 然而，数据通常都是包含一些冗余的，“高维数据中的有效成分”实际上只包含在一个比较小的 $M$ 维子空间（$R^M\subset R^N,M<N$）中。

\qquad 采用 $\text{PCA}$ 提取出最重要的 $M$ 个成分（ 最大的 $M$ 个特征值所对应的成分）组成 ${\tilde{\mathbf{x}}}_i\in R^M$ 来表示完整的 $N$ 维数据 ${\mathbf{x}}_i\in R^N$，从而实现了高维数据的降维。

\qquad

2. K-L变换 —— 基于最小误差的描述

\qquad 主成分分析，在信号处理领域也被称为 $\text{K-L}$ 变换。也就是说，一个信号可以在 $\text{K-L}$ 基上展开。

2.1 内积空间中的正交变换

\qquad 由线性代数理论可知，实对称矩阵 $S_{N\times N}$ 可以通过一个正交矩阵 $A$ 实现对角化（实对称矩阵的特征分解，可参考本文）：

$AS{A}^{-1}=AS{A}^T=\left[\begin{array}{cccc}{\lambda }_0& & & \\ & {\lambda }_1& & \\ & & \ddots & \\ & & & {\lambda }_{N-1}\end{array}\right]$

$S{p}_i={\lambda }_i{p}_i(i=0,1,\cdots ,N-1)$　　【 $p_i$ 为特征值 ${\lambda }_i$ 所对应的特征向量】

\qquad\qquad 其中，正交矩阵 $A=[p_0,p_1,\cdots ,p_{N-1}]$ 满足 $A^{-1}=A^T$

\qquad 正交矩阵 $A$ 的主要特点在于：$A$ 中各列是两两正交的。

$A^{-1}=A^T⟹A^{T}A=\mathbf{I}⟹\{\begin{array}{rlr}{p}_i^T{p}_j& =0& ,i\ne j\\ & & \\ p_i^T{p}_j& =1& ,i=j\end{array}$

\qquad 若将正交矩阵 $A$ 看成两个内积空间 $X,Y$ 之间的线性变换，那么

$A:X⟶Y$
$\mathbf{x}⟶\mathbf{y}=A\mathbf{x}$　　（正变换）

\qquad 正交变换可以保证信号在变换前后的能量保持不变，即：$\Vert \mathbf{x}{\Vert }_X=\Vert \mathbf{y}{\Vert }_Y=\Vert A\mathbf{x}{\Vert }_Y$
\qquad

2.2 信号的正交分解

\qquad 由于正交矩阵 $A$ 中的各列两两正交，$A$ 中的列构成了空间的一组“正交规范基”，若认为 $\mathbf{y}=A\mathbf{x}$ 为正变换，反变换也必然存在，即：

$\mathbf{x}=A^{-1}\mathbf{y}=A^T\mathbf{y}$　　（逆变换）

\qquad 上述逆变换意味着，一个信号 $\mathbf{x}=[x_0,x_1,\cdots ,x_{N-1}{]}^T$ 可以在一组正交基 { φ ( ⋅ , n ) } \{\varphi(\cdot,n)\} {φ(⋅,n)} 上展开，展开系数为 $\mathbf{y}=[y_0,y_1,\cdots ,y_{N-1}{]}^T$，也就是：

$x_i=\sum _{n=0}^{N-1}{y}_n\phi (i,n)$，　　其中 $\phi (i,n)$ 为变换核

逆变换矩阵 $A^T=\left[\begin{array}{ccccc}\phi (0,0)& \cdots & \phi (0,n)& \cdots & \phi (0,N-1)\\ ⋮& & ⋮& & ⋮\\ \phi (i,0)& \cdots & \phi (i,n)& \cdots & \phi (i,N-1)\\ ⋮& & ⋮& & ⋮\\ \phi (N-1,0)& \cdots & \phi (N-1,n)& \cdots & \phi (N-1,N-1)\end{array}\right]$

\qquad 例如，一维离散傅里叶变换的逆变换核为 $\phi (i,n)=e^{j\frac{2\pi in}{N}}$，那么

$\mathbf{y}=A\mathbf{x}⟶y_n=\sum _{i=0}^{N-1}{x}_i{e}^{-j\frac{2\pi in}{N}},n=0,1,\cdots ,N-1$
$\mathbf{x}=A^T\mathbf{y}⟶x_i=\frac{1}{N}\sum _{n=0}^{N-1}{y}_n{e}^{j\frac{2\pi in}{N}},i=0,1,\cdots ,N-1$

可以验证一维离散傅里叶逆变换 $A^T=[e^{j\frac{2\pi in}{N}}]$ 为正交矩阵

\qquad

2.3 K-L变换

\qquad 对于一个宽平稳的随机向量 $\mathbf{x}=[x_0,x_1,\cdots ,x_{N-1}{]}^T$，其协方差矩阵为 $S_{\mathbf{x}}=E[(\mathbf{x}-{\mathbf{\mu }}_{\mathbf{x}})(\mathbf{x}-{\mathbf{\mu }}_{\mathbf{x}}{)}^T]$，其中 ${\mathbf{\mu }}_{\mathbf{x}}=E[\mathbf{x}]$。

$\text{K-L}$ 变换记为 $\mathbf{y}=A\mathbf{x}$ —— 寻找正交矩阵 $A$ 将输入信号 $\mathbf{x}$ 变换为输出信号 $\mathbf{y}$（变换域），同时使得随机变量 $\mathbf{y}$ 的协方差矩阵 $S_{\mathbf{y}}$ 为对角阵，即：

$S_{\mathbf{y}}=A{S}_{\mathbf{x}}{A}^T=\left[\begin{array}{cccc}{\lambda }_0& & & \\ & {\lambda }_1& & \\ & & \ddots & \\ & & & {\lambda }_{N-1}\end{array}\right]$

\qquad 由于随机变量 $\mathbf{y}$ 的协方差矩阵 $S_{\mathbf{y}}$ 为对角阵，各分量之间的协方差均为零，原始信号 $\mathbf{x}$ 经过 $\mathbf{y}=A\mathbf{x}$ 变换后，完全去除了相关性。
\qquad

以图像块为例：
(1) 图像像素在空间域中相关性很强（局部区域像素通常比较接近），能量分布比较均匀
(2) 图像块经过正交变换（坐标旋转或变换）后，能量集中在少数坐标轴上，变换系数 { y n } \{y_n\} {yn​}之间的相关性近似统计独立
(3) 图像块的能量在变换域中的分布相对集中（集中在直流和少数低频系数），可采用较少的编码比特表示，达到压缩编码的目的

\qquad 在 $\text{K-L}$ 基 { ϕ 0 , ϕ 1 , ⋯   , ϕ N − 1 } \{\phi_0,\phi_1,\cdots,\phi_{N-1}\} {ϕ0​,ϕ1​,⋯,ϕN−1​} 上可以将输入信号 $\mathbf{x}$ 在变换域展开：

$\mathbf{x}=A^T\mathbf{y}=[{\varphi }_0,{\varphi }_1,\cdots ,{\varphi }_{N-1}]\mathbf{y}=\sum _{n=0}^{N-1}{y}_n{\varphi }_n$

\qquad\qquad 其中，$y_n=⟨\mathbf{x},{\varphi }_n⟩={\mathbf{x}}^T{\varphi }_n$ 是 $\mathbf{x}$ 在 ${\varphi }_n$ 上的投影

\qquad 因此
$\mathbf{x}=A^T\mathbf{y}=\sum _{n=0}^{N-1}⟨\mathbf{x},{\varphi }_n⟩{\varphi }_n$

由 $A^T=A^{-1}⟹A^{T}A=\mathbf{I}$，说明 $A$ 中各列 { p i } \{p_i\} {pi​} 两两正交
又由于 $A^T=A^{-1}⟹A{A}^T=\mathbf{I}⟹(A^T{)}^T(A^T)=\mathbf{I}$，说明 $A^T$ 中各列 { ϕ n } \{\phi_n\} {ϕn​} 也两两正交

\qquad

2.4 基于K-L变换的数据压缩

\qquad 要对信号 $\mathbf{x}$ 做数据压缩，只要舍去变换系数 $\mathbf{y}$ 的一些系数。假设只保留 $M<N$ 个系数，也就是：

$\hat{\mathbf{x}}=\sum _{n=0}^{M-1}{y}_n{\varphi }_n=\sum _{n=0}^{M-1}⟨\mathbf{x},{\varphi }_n⟩{\varphi }_n$

\qquad 压缩后数据 $\hat{\mathbf{x}}$ 对原始数据 $\mathbf{x}$ 的均方误差为：$\epsilon =E[(\hat{\mathbf{x}}-\mathbf{x}{)}^2]$

\qquad 因此，$\text{K-L}$ 变换也可以理解为：在保留 $M$ 个系数的前提下，为了获得最大压缩率，选择一组标准正交基 { ϕ 0 , ϕ 1 , ⋯   , ϕ N − 1 } \{\phi_0,\phi_1,\cdots,\phi_{N-1}\} {ϕ0​,ϕ1​,⋯,ϕN−1​}，使得上述均方误差最小。

$\begin{array}{rl}\epsilon & =E[(\hat{\mathbf{x}}-\mathbf{x}{)}^2]=E\left\{{\left[\sum _{n=M}^{N-1}{y}_n{\varphi }_n\right]}^2\right\}\\ & =E\left\{⟨\sum _{n=M}^{N-1}{y}_n{\varphi }_n,\sum _{n=M}^{N-1}{y}_n{\varphi }_n⟩\right\},⟨{\varphi }_i,{\varphi }_j⟩=0(i\ne j)\\ & =E\left\{\sum _{n=M}^{N-1}{y}_n^2\right\},⟨{\varphi }_n,{\varphi }_n⟩=1\\ & =E\left\{\sum _{n=M}^{N-1}[{\mathbf{x}}^T{\varphi }_n{]}^2\right\},[{\mathbf{x}}^T{\varphi }_n{]}^2=({\mathbf{x}}^T{\varphi }_n{)}^T({\mathbf{x}}^T{\varphi }_n)\\ & =E\left\{\sum _{n=M}^{N-1}[{\varphi }_n^T\mathbf{x}{\mathbf{x}}^T{\varphi }_n]\right\}\\ & =\sum _{n=M}^{N-1}{\varphi }_n^{T}E\left\{\mathbf{x}{\mathbf{x}}^T\right\}{\varphi }_n\end{array}$

\qquad 假设 $\mathbf{x}$ 已去除均值，那么 $E\left\{\mathbf{x}{\mathbf{x}}^T\right\}=S_{\mathbf{x}}$，那么最小均方误差：$\epsilon =\sum _{n=M}^{N-1}{\varphi }_n^T{S}_{\mathbf{x}}{\varphi }_n$

\qquad 同样采用拉格朗日乘子法构造出最优化问题：

min ⁡   {   ε + λ ( 1 − ϕ T ϕ )   } \qquad\qquad\qquad\min\ \{\ \varepsilon+\lambda(1-\phi^T\phi)\ \} min { ε+λ(1−ϕTϕ) }

\qquad 对${\varphi }_n$求偏导：

∂ ∂ ϕ n {   ε + λ ( 1 − ϕ T ϕ )   } = 0 \qquad\qquad\qquad\dfrac{\partial}{\partial \phi_n}\{\ \varepsilon+\lambda(1-\phi^T\phi)\ \}=0 ∂ϕn​∂​{ ε+λ(1−ϕTϕ) }=0

\qquad 可求得：

$S_{\mathbf{x}}{\varphi }_n={\lambda }_n{\varphi }_n,n=M,\cdots ,N-1$

\qquad
\qquad 因此，在变换域截短数据后的均方误差为：

$\epsilon =\sum _{n=M}^{N-1}{\varphi }_n^T{S}_{\mathbf{x}}{\varphi }_n=\sum _{n=M}^{N-1}{\varphi }_n^T{\lambda }_n{\varphi }_n=\sum _{n=M}^{N-1}{\lambda }_n$

\qquad 为了使均方误差最小，将协方差矩阵 $S_{\mathbf{x}}$ 的特征值按照从大到小的顺序排列，即 ${\lambda }_0\ge {\lambda }_1\ge \cdots \ge {\lambda }_{N-1}$，保留最大的 $M$ 个特征值，相当于在变换域（保留$M$个变换系数时）保留了信号的最大能量。

\qquad
\qquad 由于 $\text{K-L}$ 变换能够完全去除原始信号 $\mathbf{x}$ 的相关性，通过保留最大的 $M$ 个特征值来进行数据压缩时可以实现截短后数据 $\hat{\mathbf{x}}$ 的均方误差最小，$\text{K-L}$ 变换也被称为“最佳变换”。然而，由于 $\text{K-L}$ 变换的基取决于协方差矩阵 $S_{\mathbf{x}}$，协方差矩阵的特征值和特征向量的计算没有像 $\text{FFT}$ 之类的快速算法，限制了其应用场景。离散余弦变换 $\text{DCT}$ 在满足一阶马尔可夫过程假设的前提下，可以极好地近似 $\text{K-L}$ 变换，又有快速算法，因而在压缩编码中得到了广泛应用。

\qquad
代码实现：PCA图像压缩
实现原理：
（1）将 $M\times N$ 大小的图像进行分块，若每个图像块大小为 $B_s\times B_s$，整幅图像包含了 $N_s$ 个图像块，也就是满足 $MN=N_s{B}_s^2$
（2）相当于构建了一个 $B_s^2\times 1$ 的随机向量，整幅图像包含了该随机向量的 $N_s$ 次实现
（3）$\text{PCA}$ 针对该随机向量的 $N_s$ 次实现来展开，寻找该随机向量的 $B_s^2$ 个分量中最重要的 $p$ 个主分量
　　【从K-L变换的角度来看】：一个 $B_s^2\times 1$ 输入图像数据，经过K-L变换后，在新的 $B_s^2$ 维空间中仅保留 $p$ 个维度（其实就是降维），若采用 $\text{PCA}$ 选中的这 $p$ 个维度来表示图像，可以使得“变换后图像数据”的能量达到最大。
（4）从 $\text{K-L}$ 变换域回到原始图像的空间域，只需要在 $B_s^2$ 维变换域空间中保留 $p$ 个维度的分量值，将剩余的 $B_s^2-p$ 个分量置 $0$。

上述过程仅仅考虑了图像压缩在变换域的表现，并未考虑其他的量化、编码等过程。