PCA降维提取图像特征

PCA算法应用背景

在机器学习和数据处理中，较多的场景会处理存在多个字段的关联数据，例如某个视频网站平台的用户画像，每一次观看视频记录可能会存在用户性别、用户年龄、用户设备登陆地理位置、所观看视频时长、视频分类等。

其中某些字段为标识量而非具体的量化值，再经过业务层的数据统一与清理后可以得到能够量化与计算的向量，如：

(1,23,340102,23,003)^T

即使用0、1、2表示性别为不明确、男、女，23表示23岁，340102表示设备登录地的合肥市瑶海区的行政区划代码，23表示视频时长23分钟，003表示平台视频分区及子分类。

但是在实际的业务场景中字段数与记录数可能要比以上所示例的五维数据更多，从而带来的问题便是算法复杂度呈现指数级的增加，对于机器学习可能会需要更长时间才能够训练出符合业务需求的模型，大大降低了业务的灵活程度。

这时需要通过在不影响最终模型效果的情况下对数据进行降维，在数据降维的过程中，利用数据本身存在的关联性，使得信息损失尽量降低。

例如在上面的例子中，视频的分区与子类型和视频的长度存在一定的关联性，那么在一个有m条记录n个字段的数据集，即m×n的矩阵中，去掉视频时长的字段即矩阵中所在列，实际的信息损失对最终训练的模型没有较大影响。

PCA算法简述

基与投影值的引入

在数学中，若要准确的描述向量，那么首先要确定一组基，向量的各个分量即该向量在基的各个直线上的投影值。

在常见的二维平面直角坐标系中，我们经常默认基为(0,1) ,(1,0)，又因为线性空间中的基可以是任意两个线性无关的向量，可以变换基来达到坐标变换的目的。

在新的基下，可通过原坐标点乘新基来获得新坐标，如图所示，在原始基(1,0) ,(0,1)下，基变换为(\frac{1}{\sqrt2},\frac{1}{\sqrt2}) ,(-\frac{1}{\sqrt2},\frac{1}{\sqrt2})，坐标从原始的(2,3)变换为(\frac{5}{\sqrt2},-\frac{1}{\sqrt2}).

PCA_SVD-20200327-1.jpg

这种计算过程也可以使用矩阵表示，上图的变换过程用矩阵表示为：

\left\lgroup \begin{matrix} \frac{1}{\sqrt2} & \frac{1}{\sqrt2} \\ -\frac{1}{\sqrt2} & \frac{1}{\sqrt2} \end{matrix} \right\rgroup \times \left\lgroup \begin{matrix} 3 & \frac{1}{\sqrt2} \\ 2 & \frac{1}{\sqrt2} \end{matrix} \right\rgroup

其中第一个矩阵的两行分别为新基的两个向量，乘上坐标向量后即可得到变换后的坐标.

PCA_SVD-20200327-1.jpg

推广到一般情况下可得，如果有m个n维的向量，变换为由r个n维向量的空间中，先将新空间中的r个基按行组成矩阵A，将待变换的向量按列组成B，变换结果为AB.

\left\lgroup \begin{matrix} p_1 \\ p_2 \\ \vdots \\ p_r \end{matrix} \right\rgroup \left\lgroup \begin{matrix} a_1 & a_2 & \cdots & a_m \end{matrix} \right\rgroup = \left\lgroup \begin{matrix} p_1a_1 & p_1a_2 & \cdots & p_1a_m \\ p_2a_1 & p_2a_2 & \cdots & p_2a_m \\ \vdots & \vdots & \ddots & \vdots \\ p_ra_1 & p_ra_2 & \cdots & p_ra_m \end{matrix} \right\rgroup

其中p_1,p_2,\cdots,p_r为新空间的各个基，a_1,a_2,\cdots,a_m为待变换的向量，这里的r可以小于m，即高维的向量可以变换到低维的空间中.

PCA过程分析

通过基变换可以将一个高维的向量变换为低维的向量，那么主要问题便是如何确定新的低维空间的基才能够使得在变换过程中损失信息最低，最大的保留原有信息，下面将通过一个实例来分解这个问题.

假设某视频平台的某个数据库中现有5条记录，每个记录有两个字段，数值为方便计算与实际情况不同.

字段名	记录1	记录2	记录3	记录4	记录5
视频时长	1	1	2	4	2
用户年龄	1	3	3	4	4

将以上的记录抽象为一个矩阵并均值化

\left\lgroup \begin{matrix} 1 & 1 & 2 & 4 & 2 \\ 1 & 3 & 3 & 4 & 4 \end{matrix} \right\rgroup \xrightarrow{均值化} \left\lgroup \begin{matrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{matrix} \right\rgroup

明显可得以上5个数据为二维数据，若要将为一维的数据并且最大程度的保留信息，需要选择合适的一维空间的基，即选定合适的一维空间的方向.

首先需要满足在新的基上投影后的投影值尽可能分散，如果在投影后出现投影值重合，那么将会严重的丢失信息；

其次在上面的例子中是从二维将为一维，只需要选定一维空间中的一个方向，如果是三维降维二维或更高维降维，需要选定低维空间中的每个方向，若仍然按照在新的基上投影后的投影值尽可能分散，那么低维空间中的各个维度方向几乎是重合的，所以需要其它的约束条件.直观的说，让各个字段尽可能多的保留信息，希望它们之间不存在相关性，即线性无关.

通过以上分析可以得到两个约束条件：

• 新的基上投影后的投影值尽可能分散.
• 各个字段相关性尽可能低.

在数学上可以通过方差来表示分散程度，协方差来表示字段相关性，当协方差为0时，低维空间各个维度的方向即低维空间的基是正交的.且每个字段已经经过均质化处理，其均值为0.

投影值分散程度:Var(a)=\frac{1}{m}\Sigma^{m}_{i=1}a^2_i \\ 字段相关性:Cov(a,b)=\frac{1}{m}\Sigma^{m}_{i=1}a_ib_i

从而得到了可以使用数学表示的约束条件及降维变换目标：将一组n维向量降维至k维(0 < k < n)，需要选定k个单位正交基，使得原n维向量变换到选定的基上后，各个字段间的协方差为0，字段的方差尽可能大.

最终的变换目标与字段内的方差和字段间的协方差有关系，为了使实际计算可行，需要将二者统一表示.

假设现有的数据仍为两个字段，按行排列为矩阵X

X= \left\lgroup \begin{matrix} a_1 & a_2 & \cdots & a_m \\ b_1 & b_2 & \cdots & b_m \end{matrix} \right\rgroup \frac{1}{m}XX^T= \left\lgroup \begin{matrix} \displaystyle\frac{1}{m}\sum^{m}_{i=1}a^2_i & \displaystyle\frac{1}{m}\sum^{m}_{i=1}a_ib_i\\ \displaystyle\frac{1}{m}\sum^{m}_{i=1}a_ib_i & \displaystyle\frac{1}{m}\sum^{m}_{i=1}b^2_i \end{matrix} \right\rgroup

此时每个字段内的方差与字段间的协方差已经被统一表示在了一个矩阵内，可以称矩阵\frac{1}{m} XX^T为协方差矩阵，推广到一般情况可得：若有m个n维的记录，将这些数据按列组成n×m矩阵X，C=\frac{1}{m} XX^T，C为对称矩阵，对角线为每个字段的方差，其余第i行第j列、第j行第i列的元素为i字段与j字段的协方差.

根据约束条件及降维变换目标，可以发现将协方差矩阵C对角化即可，以下为推导过程.

原始数据矩阵为X，X的协方差矩阵\frac{1}{m}XX^T为C，新空间的基按行组成矩阵P，降维变换到新空间后的数据矩阵为Y，Y的协方差矩阵\frac{1}{m}YY^T为D

C=\frac{1}{m}XX^T Y=PX D=\frac{1}{m}YY^T

降维变换后的数据矩阵Y已经满足约束条件及降维变换目标，所以Y的协方差矩阵D为对角矩阵，且D与原始矩阵X存在关系，可做如下推导

D=\frac{1}{m}YY^T=\frac{1}{m}(PX)(PX)^T =\frac{1}{m}PXX^TP^T=P(\frac{1}{m}XX^T)P^T=PCP^T

那么问题简化为寻找一个矩阵P使得PCP^T为对角阵，所求的P前k行即为所要寻找的基.

又因为C为实对称矩阵，实对称矩阵具有以下性质

• 实对称矩阵各个特征值对应的特征向量正交
• 若特征向量λ有r重，存在r个线性无关的的特征向量对应λ，这r个线性无关的特征向量正交

那么n×n的实对称矩阵一定有n个单位正交特征向量，这些特征向量设为e_1,e_2,\cdots,e_n，按列组成矩阵

E=(e_1,e_2,\cdots,e_n)

原始数据矩阵的协方差矩阵C有

E^TCE= \left\lgroup \begin{matrix} \lambda_1 & & & \\ & \lambda_1 & & \\ & & \ddots & \\ & & & \lambda_1 \end{matrix} \right\rgroup

即可得到P=E^T，P为协方差矩阵特征向量单位化后按行列组成的矩阵，每一行为C的特征向量

PCA算法的数学实例

通过上一节的分析，得到如下PCA算法计算步骤

• 原始数据按列组成n×m的矩阵X
• 对X每一行数据均值化
• 求协方差矩阵C=\frac{1}{m} XX^T
• 求协方差矩阵C的特征值及对应的特征向量
• 将特征向量对应特征值从大到小按行排列，取前k行组成矩阵P
• 降维后的数据即为Y=PX

仍使用上一节的原始数据

X= \left\lgroup \begin{matrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{matrix} \right\rgroup

此时数据已经均值化，求取协方差矩阵

C=\frac{1}{5} \left\lgroup \begin{matrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{matrix} \right\rgroup \left\lgroup \begin{matrix} -1 & -2 \\ -1 & 0 \\ 0 & 0 \\ 2 & 1 \\ 0 & 1 \end{matrix} \right\rgroup = \left\lgroup \begin{matrix} \displaystyle\frac{6}{5} & \displaystyle\frac{4}{5} \\ \\ \displaystyle\frac{4}{5} & \displaystyle\frac{6}{5} \end{matrix} \right\rgroup

求C的特征值

\lambda_1=2,\lambda_2=\frac{2}{5}

对应的特征向量为，其中c_1,c_2可取任意实数

c_1\left\lgroup \begin{matrix} 1 \\ 1 \end{matrix} \right\rgroup , c_2\left\lgroup \begin{matrix} -1 \\ 1 \end{matrix} \right\rgroup

特征向量标准化后，按行排列为P

\left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} \\ \displaystyle\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup , \left\lgroup \begin{matrix} \displaystyle-\frac{1}{\sqrt{2}} \\ \displaystyle\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup , P= \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} \\ \displaystyle-\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup

验证协方差矩阵C对角化

PCP^T= \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} \\ \displaystyle-\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup \left\lgroup \begin{matrix} \displaystyle\frac{6}{5} & \displaystyle\frac{4}{5} \\ \\ \displaystyle\frac{4}{5} & \displaystyle\frac{6}{5} \end{matrix} \right\rgroup \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle-\frac{1}{\sqrt{2}} \\ \displaystyle\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup = \left\lgroup \begin{matrix} 2 & 0 \\ 0 & \displaystyle\frac{2}{5} \end{matrix} \right\rgroup

用矩阵P的第一行乘以原始数据矩阵X，降维数据为

Y= \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup \left\lgroup \begin{matrix} -1 & -2 \\ -1 & 0 \\ 0 & 0 \\ 2 & 1 \\ 0 & 1 \end{matrix} \right\rgroup = \left\lgroup \begin{matrix} \displaystyle-\frac{3}{\sqrt{2}} & \displaystyle-\frac{1}{\sqrt{2}} & 0 & \displaystyle\frac{3}{\sqrt{2}} & \displaystyle-\frac{1}{\sqrt{2}} & \end{matrix} \right\rgroup

如图所示，数据已经从二维降维为一维

PCA_SVD-20200327-1.jpg

PCA算法图像处理中的实际应用

原始数据的数据集使用英国剑桥大学的Olivetti研究实验室公开数据集ORL_faces，此数据集为该实验室在1992年4月至1994年4月期间拍摄的40个人在不同角度和不同光照条件下的共400张人脸图像，并已经做了分辨率统一与灰度化处理，以下为部分图像展示

PCA_SVD-20200327-4.png

下面将使用PCA算法提取这400张图像中每个图像的特征，首先需要将原始数据集进行输入处理

PCA_SVD-20200327-4.png

原始数据集中的图片文件为PGM格式，相对于主流的图片格式，此格式不会对数据进行压缩，并且图片文件中的每个像素均使用数字按照像素实际行列位置进行保存，数字大小代表该像素的灰度.

假设原始数据集为10张照片，每张图片为100×100的分辨率，即有10000个数值用来保存该图像信息，那么该图像即为一个100×100的矩阵.将每一个图片的所有行首尾相连并按列排列，串成一个一维的10000×1的向量，将所有图片的列向量组成一个10000×10的矩阵X.该矩阵即可表示10张图片，其中每一列为一张图片.

原始数据集转换为矩阵X后，对X的每一行求均值，然后将X减去均值，即完成矩阵的均值化.

计算X的协方差矩阵C=XX^T，但在实际计算中，X有10000行，其协方差矩阵为10000×10000大小，计算量过大，将会导致计算机极大的消耗内存，这里可以通过计算C^T=X^T X，这里的C^T为10×10大小，设E'为C^T的特征值，那么C的特征值即为XE'.

通过以上方法计算得到矩阵C的特征值与特征向量，将特征向量还原为二维图像后看起来与人脸相似，所以C的每个特征向量即为一个特征脸.下图为从原始400张图片中通过PCA算法求得的部分特征脸.

PCA_SVD-20200327-4.png

所求得特征脸在实际的业务逻辑中无需还原为二维图像用来展示，而是保留部分特征脸用于代替大多数脸的近似组合，人脸可以通过一系列的向量进行保存，相较于原始的图像文件，能够节省大量的存储空间，也便于后续的模型训练.

SVD分解图像压缩

SVD算法简述

假设现有一实数域上的任意m×n矩阵A，这个矩阵可以使用SVD分解来表示

PCA_SVD-20200327-4.png

其中U是一个m×m的矩阵，其向量为正交的，称为左奇异向量；Σ是一个m×n的矩阵，其对角线上的元素为奇异值，其余元素均为0；V是一个n×n的矩阵，其向量为正交的，称为右奇异向量.

直接计算U,V矩阵比较困难，可以利用如下性质

AA^T=UΣV^T VΣ^T U^T=UΣΣ^T U^T \\ A^T A=VΣ^T U^T UΣV^T=VΣ^T ΣV^T

那么在计算矩阵U时，对矩阵AA^T进行特征分解，即求AA^T的特征值与其对应的特征向量，同理在计算矩阵V时，求A^T A的特征值与其对应的特征向量.计算矩阵Σ是将在计算U或V时的非零特征值从大到小排列后开根号作为对角线的元素.

在大多数情况下，计算得出的奇异值中（经过从大到小排列），前10%的奇异值的和占据了全部奇异值的和90%以上.如下图，在一个600×1200的矩阵中，绘制出奇异值数值与奇异值个数的关系，可以发现奇异值数值下降的非常迅速，能够通过前几个奇异值来表示原始矩阵的信息.

PCA_SVD-20200327-4.png

因此可以通过取前r个大的奇异值来近似等于原始矩阵，即

A_{m×n}≈U_{m×r} Σ_{r×r} V_{r×n}^T

PCA_SVD-20200327-4.png

SVD算法数学实例

通过上一节的分析，可得SVD分解矩阵并压缩的步骤

• 计算矩阵AA^T的特征值与其对应的特征向量，特征向量组成U
• 计算矩阵A^T A的特征值与其对应的特征向量，特征向量组成V
• 将第一步或第二部中求得的特征值从小到大排列后开根号，取前r个作为Σ的对角线元素，其余元素为0
• U矩阵取前r列，V^T矩阵取前r行，完成压缩

假设现有原始矩阵A，使用SVD分解对A进行压缩

A= \left\lgroup \begin{matrix} 1 & 1 \\ 1 & 1 \\ 0 & 0 \end{matrix} \right\rgroup

计算AA^T,A^TA

AA^T= \left\lgroup \begin{matrix} 2 & 2 & 0 \\ 2 & 2 & 0 \\ 0 & 0 & 0 \end{matrix} \right\rgroup , A^TA= \left\lgroup \begin{matrix} 2 & 2 \\ 2 & 2 \end{matrix} \right\rgroup

计算AA^T的特征值和对应的特征向量为

\lambda_1=4 , \lambda_2=0 , \lambda_3=0 \\ \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} & 0 \end{matrix} \right\rgroup ^T , \left\lgroup \begin{matrix} \displaystyle-\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} & 0 \end{matrix} \right\rgroup ^T , \left\lgroup \begin{matrix} 0 & 0 & 1 \end{matrix} \right\rgroup ^T \\ \therefore U = \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle-\frac{1}{\sqrt{2}} & 0 \\ \displaystyle\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} & 0 \\ 0 & 0 & 1 \end{matrix} \right\rgroup

计算A^TA的特征值和对应的特征向量为

\lambda_1=4 , \lambda_2=0 \\ \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup ^T , \left\lgroup \begin{matrix} \displaystyle-\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup ^T \\ \therefore V = \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle-\frac{1}{\sqrt{2}} \\ \displaystyle\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup

将第一步求得的非零特征值从大到小排列并开根号，作为Σ的对角线元素，其余元素为0，即

\Sigma= \left\lgroup \begin{matrix} 2 & 0 \\ 0 & 0 \\ 0 & 0 \end{matrix} \right\rgroup

原始矩阵A的SVD分解如下

A=U\Sigma V^T= \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle-\frac{1}{\sqrt{2}} & 0 \\ \displaystyle\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} & 0 \\ 0 & 0 & 1 \end{matrix} \right\rgroup \left\lgroup \begin{matrix} 2 & 0 \\ 0 & 0 \\ 0 & 0 \end{matrix} \right\rgroup \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle-\frac{1}{\sqrt{2}} \\ \displaystyle\frac{1}{\sqrt{2}} & \displaystyle\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup ^T = \left\lgroup \begin{matrix} 1 & 1 \\ 1 & 1 \\ 0 & 0 \end{matrix} \right\rgroup

对矩阵进行压缩，取前1个奇异值，可得

A\approx \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}}\\ \displaystyle\frac{1}{\sqrt{2}}\\ 0 \end{matrix} \right\rgroup (2) \left\lgroup \begin{matrix} \displaystyle\frac{1}{\sqrt{2}} & \displaystyle-\frac{1}{\sqrt{2}} \end{matrix} \right\rgroup ^T

SVD算法图像处理中的实际应用

下面将使用SVD分解算法对图像进行压缩，在图像转化为矩阵时无法将每个像素的颜色信息表示出来，只能够表示一个维度，即颜色的深浅.首先将原始图片拆分成红、绿、蓝（RGB）通道的灰度图，如下图所示.

PCA_SVD-20200327-10.png

在分离为3个单色通道图片后，即可使用SVD分解算法对单色图像所转化的矩阵进行SVD分解，并取前r个的奇异值来达到压缩的目的，r值的不同得到的压缩程度也就不同，以下为不同压缩程度的对比，其中origin标签为原始图像，不同百分比代表文件大小为原始文件大小的比值

PCA_SVD-20200327-10.png

参考文献

[1] 李昕怡.基于改进Kernel_PCA的图像去噪算法[D].西安:西安电子科技大学,2017

[2] 项晓丽,武圣,许一菲.二阶分块PCA的人脸特征提取方法研究[J].测控技术,2016,35(9):25-32.

[3] 徐竟泽.融合PCA_LDA和SVM算法的人脸识别[J].计算机工程与应用,2019,55(18):34-37.

[4] 王萼芳,石生明.高等代数[M].高等教育出版社:北京,2013:26-124.

[5] Jan Erik Solem.Python计算机视觉编程[M].人民邮电出版社:北京,2014:4-29.

[6] Google.Convolutional Neural Network (CNN)|TensorFlow Core[EB/OL].https://tensorflow.google.cn/tutorials/images/cnn,2019.

[7] Nvidia.Deep Learning SDK Documentation[EB/OL].https://docs.nvidia.com/deeplearning/sdk/index.html,2019.