PCA数学原理及实现与实战（基于PCA的图像压缩）

目录
PCA 数学原理
PCA 接口实现
基于 PCA 的图像压缩
项目地址

主成分分析 数学原理

特征值和特征向量

特征值和特征向量由以下公式求出：

$$Av=\lambda v$$

左式表示对向量 $v$ 执行线性变换，右式表示对向量 $v$ 执行拉伸或缩放，即向量 $v$ 在经过线性变换后依旧与变换前共线，其称为特征向量。$\lambda$ 描述了特征向量被拉伸或缩放的程度，其称为特征值。

主成分分析

我们从最小均方误差角度来解释主成分分析的数学原理。

首先定义好推导用到的数学符号：

数学符号

符号	意义
$(x,y)$	$x$ 行 $y$ 列
$d$	样本原维数
$d^{\prime }$	降维后维数
$n$	样本数
$x^{(i)}$	第 $i$ 个样本，形如 $(d,1)$
$z_i$	第 $i$ 个投影到低维空间后的样本，形如 $(d^{\prime },1)$
$u_i$	第 $i$ 个正交基向量，形如 $(d,1)$
$U$	基向量组成的矩阵 $(u_1,u_2,\dots ,u_{d^{\prime }})$，形如 $(d,d^{\prime })$
$S$	${\sum }_{i=1}^n{x}^{(i)}{x^{(i)}}^T$
$C$	常数项

假设所有样本都已进行去中心化，且所有基向量为标准正交基向量：

$$\begin{array}{rl}||u_i|{|}^2& =1\\ \forall i,u_i^T{u}_i& =1\\ \forall i\ne j,u_i^T{u}_j& =0\end{array}$$

需要注意的是所有正交基向量应形如 $(d,1)$ 而非 $(d^{\prime },1)$。

为了保证降维后重建样本的质量，我们需要最小化重建样本和原样本的最小均方误差。重建样本可以写作：

$${\hat{x}}^{(i)}=U{z}_i$$

优化目标如下：

$$\begin{array}{rl}\mathcal{L}& =\sum _{i=1}^n||x^{(i)}-{\hat{x}}^{(i)}|{|}^2\\ & =\sum _{i=1}^n||x^{(i)}-U{z}_i|{|}^2\\ & =\sum _{i=1}^n||x^{(i)}|{|}^2-2\sum _{i=1}^n{z}_i^T{U}^T{x}^{(i)}+\sum _{i=1}^n{z}_i^T{U}^{T}U{z}_i\\ & =\sum _{i=1}^n{z}_i^T{z}_i+-2\sum _{i=1}^n{z}_i^T{U}^T{x}^{(i)}+C\end{array}$$

千万注意，$U$ 并非正交矩阵！！！也就是说，$z_i^T{U}^{T}U{z}_i=z_i^T{z}_i$ 并不是因为 $U^{T}U$ 为单位矩阵。

给定所有基向量 $u_1,u_2,\dots ,u_{d^{\prime }}$ 的情况下我们可以通过求导最小化以上目标 $\mathcal{L}$:

$$\begin{array}{rl}\frac{\partial \mathcal{L}}{\partial z_i}& =-2U^T{x}^{(i)}+2z_i=0\\ z_i& =U^T{x}^{(i)}\end{array}$$

将 $z_i=U^T{x}^{(i)}$ 代回优化目标并忽略常数项，得：

$$\begin{array}{rl}\mathcal{L}& =-{x^{(i)}}^{T}U{U}^T{x}^{(i)}\\ & =-\sum _{i=1}^n\sum _{j=1}^{d^{\prime }}{u}_j^T{x}^{(i)}{x^{(i)}}^T{u}_j\\ & =-\sum _{j=1}^{d^{\prime }}{u}_j^T\underset{S}{\underbrace{\sum _{i=1}^n{x}^{(i)}{x^{(i)}}^T}}{u}_j\end{array}$$

所以我们的优化问题为：

$$\begin{array}{r}\min -\sum _{j=1}^{d^{\prime }}{u}_j^{T}S{u}_j\\ s.t.u_i^T{u}_i=1,\forall i;\\ u_i^T{u}_j=0,\forall i\ne j;\end{array}$$

应用拉格朗日乘子法减少以上有约束问题的约束：

$$\begin{array}{r}\min \mathcal{W}=-\sum _{j=1}^{d^{\prime }}{u}_j^{T}S{u}_j+\sum _j^{d^{\prime }}{\lambda }_j(u_i^T{u}_i-1)\\ s.t.u_i^T{u}_j=0,\forall i\ne j;\end{array}$$

我们可以先暂时忽略约束，对优化目标求导：

$$\begin{array}{rl}\frac{\partial \mathcal{W}}{\partial u_i}& =-2S{u}_i+2\lambda u_i\\ S{u}_i& =\lambda u_i\end{array}$$

结果发现任意基向量 $u_i$ 都是矩阵 $S={\sum }_{i=1}^n{x}^{(i)}{x^{(i)}}^T$ 的特征向量！！！又因为半正定矩阵 $S$ 的特征向量满足约束 $u_i^T{u}_j=0,\forall i\ne j$，所以我们可以大胆地将这一结论作为最终优化问题的解。注意到我们取最大的几个特征值 ${\lambda }_i$ 的时候能使优化目标 $\mathcal{W}$ 最小——而这就是主成分分析步骤的由来。

接下来我们根据以上推导可以给出主成分分析的步骤：

1. 去中心化($x^{(i)}=x^{(i)}-\frac{1}{n}{\sum }_{i=1}^n{x}^{(i)}$)，得到样本矩阵 $X=(x^{(1)},x^{(2)},\dots ,x^{(n)})$ 进行；
2. 计算基向量矩阵 $U$，即计算 $S=X{X}^T$ 的特征值，将特征向量按照对应的特征值降序排列，得到 $U=(u_1,u_2,\dots ,u_{d^{\prime }})$；
3. 计算降维后样本矩阵 $Z=(z_1,z_2,\dots ,z_n)$，显然有 $Z=U^{T}X$。
4. 计算重建样本矩阵 $\hat{X}=({\hat{x}}_1,{\hat{x}}_2,\dots ,{\hat{x}}_n)$，显然有 $\hat{X}=UZ$。
5. 莫忘了逆去中心化！！！

实际实验中，为了迎合 pytorch 库的接口要求，我们排列样本的形式会有所不同，简单来说就是涉及的所有矩阵都转置了（大概）。此外，代码用到的符号也未必和上边推导的一致。

---

再随手附上矩阵分解的一些数学原理：

特征值分解

若矩阵 $A$ 是一个 $n\ast n$ 的方阵，矩阵 $A$ 可以如下分解：

$$A=W\Sigma W^{-1}$$

其中 $W$ 为 $n$ 个特征向量张成的 $n\ast n$ 的方阵，$\Sigma$ 为 $n$ 个特征值位于对角线上的对角阵。我们一般把特征向量和特征值按特征值降序排列。如果我们把所有特征向量进行标准化，即 $||w_i|{|}_2=1$，此时有 $W^T=W$ ，特征分解可以写作：

$$A=W\Sigma W^T$$

注意到特征值分解要求矩阵必须为方阵，此外的情况就需要我们的奇异值分解登场了！

---

奇异值分解

奇异值分解如下：

$$M=U\Sigma V^T$$

其中中 $M$ 为 $m\ast n$ 的矩阵，$U$ 为 $m\ast m$ 的矩阵，$V$ 为 $n\ast n$ 的矩阵。我们可以推出：

$$\begin{array}{rl}& \\ M^{T}M& =(U\Sigma V^T{)}^{T}U\Sigma V^T\\ & =V\Sigma U^{T}U\Sigma V^T\\ & =V\Sigma \Sigma V^T\\ & =VL{V}^T\end{array}$$

其中 $L$ 为特征值对角矩阵。类似地可以推出：

$$\begin{array}{rl}M{M}^T& =UL{U}^T\end{array}$$
可见 $U$ 为 $M{M}^T$ 的特征向量阵，$V$ 为 $M^{T}M$ 的特征向量阵，奇异值矩阵 $\Sigma$ 为 $M{M}^T$ 或 $M^{T}M$ 开方得到的对角阵。

---

基于 PCA 的图像压缩

1. 熟悉并掌握主成分分析的基本原理
2. 学会应用主成分分析实现数据降维，并应用到图像压缩

一、PCA 函数接口实现

本次实验我们使用 pytorch 库来实现 PCA 接口。代码非常简单，除去注释和形状检查剩不下几行。PCA 的实现步骤我们上文已提过了，实验中唯一的雷点在于 torch.linalg.eig 函数会修改实参的值（找这个bug找了半天）。

class PCA():
    def __init__(self, n_component: int = 16, device='cpu') -> None:
        """主成分分析

        Args:
            n_component (int): 保留的主成分数
        """
        super().__init__()
        self.n_component = n_component
        self.device = device

    def CHECK_SHAPE(self, shape: torch.Size) -> None:
        assert len(shape) >= 2, 'Shape of input is expected bigger than 2!!!'
        limit = 1
        for i in range(1, len(shape)):
            limit *= shape[i]
        assert limit >= self.n_component, f'n_component = {self.n_component}, expected <= {limit}'

    @torch.no_grad()
    def fit(self, X: torch.Tensor) -> None:
        """提取主成分

        Args:
            X (torch.Tensor): 待进行主成分分析的输入张量，形状应当为 (batch_size, ...)
        """
        self.CHECK_SHAPE(X.shape)
        Y = X.reshape(X.shape[0], -1).to(self.device)
        self.mean = Y.mean(0)
        Z = Y - self.mean
        
        covariance = Z.T @ Z
        _, eig_vec = eig(covariance)

        self.components = eig_vec[:, :self.n_component]

    @torch.no_grad()
    def transform(self, X: torch.Tensor) -> torch.Tensor:
        """数据降维

        Args:
            X (torch.Tensor): 待降维数据，形状应当为 (batch_size, ...)

        Returns:
            torch.Tensor: 降维后数据
        """
        self.CHECK_SHAPE(X.shape)
        Z = X.reshape(X.shape[0], -1).to(self.device)

        return (Z - self.mean) @ self.components.real
    
    @torch.no_grad()
    def reconstruct(self, X: torch.Tensor) -> torch.Tensor:
        """高维数据重建

        Args:
            X (torch.Tensor): 待重建数据，形状应当为 (batch_size, ...)

        Returns:
            torch.Tensor: 重建后数据
        """
        assert len(X.shape) == 2, 'Shape of input is expected to equal to 2!!!'

        return (X @ self.components.real.T) + self.mean

二、(a)：对应代码中的 main_pr1.py

我们对 5000 张 32 * 32 人脸灰度图 (./data/faces.mat) 进行主成分分析，分别降维至 (10, 50, 100, 150) 维，并截取部分人脸的重建图片用于对比，如下所示：

原始图片 (5000, 32, 32)=(5000, 1024)

(5000, 150)

(5000, 100)

(5000, 50)

(5000, 10)

可以看到，随着降维的维度越来越低，图片越来越糊。

PSNR(Peak Signal-to-Noise Ratio) 对比

给定干净图像 $I$ 和噪声图像 $K$，PSNR 定义为：

$$\begin{gathered} MSE=\frac{1}{mn}\sum _{i=0}^{m-1}\sum _{j=0}^{n-1}||I(i,j)-K(i,j)|{|}^2 \\ PSNR=10{\log }_{10}\left(\frac{MA{X}_I^2}{MSE}\right) \end{gathered}$$
其中 $MA{X}_I^2$ 为图像 $I$ 的最大灰度值。实验中我们直接调用 skimage.metrics.peak_signal_noise_ratio 函数。
对比结果如下：

avg_psnr_top_10: 21.333251165491017
avg_psnr_top_50: 25.304461936942744
avg_psnr_top_100: 28.19060614744335
avg_psnr_top_150: 30.383015243548055

PSNR 值越大，代表图片失真越小。可以看到随着降维维数的减小，5000 张图片的平均 PSNR 值越来越小。

二、(b)：对应代码中的 main_pr2.py

我们用 PCA 对 data/lena.jpg 进行降维并重建，同样是分别降维至 (10, 50, 100, 150) 维。直接上结果：

原始图片 (512, 512, 3) = (512, 1536)

(512, 150)

(512, 100)

(512, 50)

(512, 10)

同样可以看到，随着降维程度增大，图片越来越糊，10 维时已没了人样。

PSNR 定量对比

avg_psnr_top_10: 21.957091434918926
avg_psnr_top_50: 27.010780484659534
avg_psnr_top_100: 31.253271060749288
avg_psnr_top_150: 34.0504621750583

同样地，随着降维维数的减小，图片的 PSNR 值越来越小，图片失真程度越来越严重。

二、©：对应代码中的 main_pr3.py

我们将 PCA 降维应用到手写数字识别任务上，对图像进行压缩降维，再训练聚类模型。我们实验中的聚类模型为 KMEANS，直接拿上次实验用 pytorch 实现的 KMEANS 来调用，测试集则选用 torchvision 中手写数字集的测试集，包含 10000 张图片。图片降维后，KMEANS 的训练速度要显著提升——这是显而易见的，降维后需要处理的数据量大大减少了。

聚类正确率分析

我们将五次运行的正确率取平均：

降维前	降维到150维后
0.5296	0.53415

降维前	降维到100维后
0.5553	0.58132

降维前	降维到50维后
0.537	0.576

降维前	降维到10维后
0.4856	0.539

可以看到，大多时候聚类正确率有了一定的提升。这可能是因为 PCA 去除无用的了噪声与冗余，通过消除特征间的相关性使得距离度量能更好地反映样本间的实际相似度。然而实际上很多时候降维前后聚类的正确率似乎没有大的变动，PCA 对于聚集不同数据分布的效果十分有限。

KMEANS 聚类结果降维可视化

PCA 可用于高维数据可视化。我们将原始数据维度的 KMEANS 聚类结果降维至二维并呈现在一张图上。

K = 3

K = 4

K = 5

K = 10

可以看到，同类的数据点降维至二维后大多汇聚在一起，这说明我们的聚类和 PCA 都确乎起到了作用。

---

环境依赖

实验中使用 python 3.10.0，但更低版本的 python 应该也没问题。

torch==1.13.1
torchvision==0.14.1
numpy==1.26.3
tqdm==4.66.1
scipy==1.12.0
matplotlib==3.8.2
skimage==0.23.2

运行命令

确保环境依赖安装完成后，前俩小实验都可以直接跑，最后一个实验需要指定一些参数（前俩实验也可以指定一些参数，但需要直接修改代码）。详见代码。

python main_pr1.py
python main_pr2.py
python main_pr1.py --eval --draw --num_clust 10 --dataset './ur_path_to_MNIST' --n_components 100

项目地址

https://gitee.com/laozibabac/pca-pytorch.git