线性降维：主成分分析原理及仿真

莺嘴啄花红溜，燕尾点波绿皱。
指冷玉笙寒，吹彻小梅春透。
依旧，依旧，人与绿杨俱瘦。
——《如梦令·春景》 秦观

更多精彩内容请关注微信公众号 “优化与算法”

一、背景

随着信息技术的发展，数据量呈现爆照式增加，高维海量数据给传统的数据处理方法带来了严峻的挑战，所以，开发高效的数据处理技术是很是必要的。数据降维是解决维度灾难的一种有效方法，之因此对数据进行降维是由于：

1. 在原始的高维空间中，包含有冗余信息以及噪音信息，在实际应用例如图像识别中形成了偏差，下降了准确率；而经过降维,咱们但愿减小冗余信息所形成的偏差，提升识别（或其余应用）的精度。
2. 但愿经过降维算法来寻找数据内部的本质结构特征。

主成分分析（Principal Component Analysis，PCA）是一种常见的线性降维方法，普遍应用于图像处理、人脸识别、数据压缩、信号去噪等领域。

二、算法思想与原理

设原数据大小为 \(N\times M\) ，通过PCA降维后的数据大小为 \(N \times K\)，其中 \(K<M\)。
PCA的中心思想是：将高维数据投影到具备最大方差的低维空间中。

2.1 从不一样角度理解PCA

PCA原理能够从不一样角度来理解：

1. 从空间变换的角度来看，这里咱们把 \(N\) 称为数据维度，\(M\) 称为原子个数，所以，PCA也能够直白的理解为将 \(N\) 维空间中相关性较大的 \(M\) 个原子组成的空间压缩成相关性最小的 \(K\) 个原子组成的空间。
2. 从分类的角度来看，这里咱们把 \(N\) 称为样本数量，\(M\) 称为样本类数，所以，PCA也能够普通的理解为将 \(M\) 个样本类别从新整理成具备最大区分度的 \(K\) 个类别。
3. 从信号表示的角度来看，这里咱们把 \(N\) 称为单个信号长度，\(M\) 称为信号数量，所以，PCA也能够粗略的理解为将 \(M\) 个很类似的信号经过适当的投影（或者说变换\映射）用 \(K\) 个信号来简单的表示原来的。就好比说一只猫的叫声为“喵喵喵~”，另外一只猫的叫声也为“喵喵喵~”，虽然他们的叫声会存在必定的差别，可是对于人类来讲几乎是同样的，这两只猫发出的声音信号能够直接用一个信号来表示。
4. 从特征提取的角度来讲，这里咱们把 \(N\) 称为数据维度，\(M\) 称为特征数量，所以，PCA也能够近似的理解为将数据集中 \(M\) 个特征用 \(K\) 个特征来表示，其中去掉的 \(M-K\) 个特征是可有可无的。

2.2 PCA直观思想

为了方便理解，下面以一个具体实例来讲明PCA的这种思想。
假设某地的房子面积和总价在理想状况下符合以下正比关系：

表1. 理想房子面积与总价关系

这些数据在二维平面中能够表示为：

图1. 理想房子面积与总价关系

然而在实际中并不是如此，因为各类缘由总会产生一些波动，有可能会是如表2所示的关系。

表2. 实际房子面积与总价关系

一样的，这些数据在二维平面中能够表示为：

图2. 实际房子面积与总价关系

若把表一中的数据看成一个矩阵来看，就变成：

表3. 表1中数据的矩阵形式

从表1中的数据来看，在理想状况下，房子面积和总价符合正比关系，也就是说房子面积和总价具备很大的相关性，若把它当成一个 \(2\times 5\) 的矩阵的话，它的秩为1，这说明用表1中的10个数据来描述房子面积和总价的这种关系显得有些冗余，咱们只须要一行数据（\(5\)个）就能够描述清楚这种关系。在实际状况中，表2中的两行数据虽然不知足彻底的正比关系，可是它们也很是接近正比关系，所以，他们之间的相关性也是很是大的。那么咱们如何将 \(\bf{X}\) 中的数据用一行数（也就是从二维降到一维）来表示呢？由前面相关性分析可知，应该是要进行“去相关性”操做，也就是要使两行数据经过一个变换，使其相关性变小。

数学里一般用方差和协方差来描述数据的相关性，方差定义以下：

\[Var({\bf{x}}) = \frac{1}{N}\sum\limits_{i = 1}^N {{{\left( {{x_i} - E({\bf{x}})} \right)}^2}} \]

其中 \(E({\bf{x}}) = \frac{1}{N}\sum\limits_{i = 1}^N {\left( {{x_i}} \right)}\)为向量 \(\bf{x}\) 的指望。

很明显，方差是描述一个向量（维度）内数据的分散程度（相关性）的度量，若若是数据是像表3同样的矩阵形式呢？那么就用协方差来描述数据之间的相关性：

\[Cov( \bf{x},\bf{y}) \ =\ \frac{1}{N}\sum ^{N}_{i=1}( x_{i} \ -E(\bf{x}))( y_{i} -E(\bf{y})) =E[( \bf{x}-E(\bf{x}))( \bf{y}-E(\bf{y}))\]

数量高于2个的 \(N\) 维数据一般用协方差矩阵来表示:

\[ \bf{C}\ =\begin{pmatrix} cov(\bf{x},\bf{x}) & cov( \bf{x},\bf{y})\\ cov( \bf{y},\bf{x}) & cov( \bf{y},\bf{y}) \end{pmatrix}\]

假设数据已经通过中心化处理，数据矩阵 \(\bf{X}\) 的协方差矩阵就能够写成：

\[ \bf{C_{\bf{X}}}\ =\begin{pmatrix} cov(\bf{X}_{1},\bf{X}_{1}) & cov(\bf{X}_{1},\bf{X}_{2})\\ cov( \bf{X}_{2},\bf{X}_{1}) & cov( \bf{X}_{2},\bf{X}_{2}) \end{pmatrix}\]

其中，\(\bf{X}_{1}\) 和 \(\bf{X}_{2}\) 分别为 \(\bf{X}\) 的第一行和第二行。
在这个例子中，协方差具体表示什么意义呢？

图3. 协方差的意义

如图3所示，将各个数据点分别对横轴和纵轴作投影，其中绿色线与横轴的交点为数据在横轴上的投影，蓝线与纵轴的交点为数据在纵轴上的投影。那么横轴上交点的分散程度（方差）就对应了 \(\bf{C}\) 中的 \(cov(\bf{x},\bf{x})\)，而纵轴上的交点的分散程度就对应了 \(\bf{C}\) 中的 \(cov(\bf{y},\bf{y})\)。\(cov( \bf{y},\bf{x}\) 则表示 \(\bf{x}\) 和 \(\bf{y}\) 之间的相关程度。

那么由PCA的中心思想：将高维数据投影到具备最大方差的低维空间中，咱们就能够知道只要将数据投影到方差最大的低维（\(K\) 维）空间中就行，那么如何找到方差最大的低维空间呢？如图4所示，数据不只能够投影到横轴和纵轴，还能够360度旋转，投影到任何直线（1维空间）上，对于图4中的 \(pc1\) 方向和 \(pc2\) 方向，能够直观的看出来确定要投影到 \(pc1\) 方向要好一点，由于投影以后的数据点更分散。

数据投影到一个维度越分散在必定程度上就说明两类数据的相关性越小，这就指导咱们一条思路：只须要将 \(\bf{C_{\bf{X}}}\) 中的两个协方差 \(cov(\bf{X}_{1},\bf{X}_{2})\) 和 \(cov(\bf{X}_{2},\bf{X}_{1})\) 变成0就能够了，这就意味着将 \(\bf{C_{\bf{X}}}\) 变成对角矩阵就好了。

图4. 如何选择投影方向？

咱们发现 \(\bf{C_{\bf{X}}}\) 是一个对称矩阵，天然就会想到特征值分解。

\(\bf{C_{\bf{X}}}\) 是 \(\bf{M}\) 阶对称矩阵，则必有正交矩阵 \(\bf{P}\)使得：
\({{\bf{P}}^{ - 1}}{{\bf{C}}_{\bf{X}}}{\bf{P}} = {{\bf{P}}^{\rm T}}{{\bf{C}}_{\bf{X}}}{\bf{P}} = {\bf{\Lambda }}\)
其中 \(\bf{\Lambda }\) 是以 \(\bf{C_{\bf{X}}}\) 的 \(\bf{M}\) 个特征值为对角元

所以咱们只要找到这个正交矩阵 \(\bf{P}\) 和它对应的特征值，而后取前 \(K\) 个最大的特征值及其对应的特征向量就能够了，这实际上就是进行特征值分解：

\[{{\bf{C}}_{\bf{X}}} = {{\bf{P}}^{\rm T}}{\bf{\Lambda P}} \]

固然也能够经过奇异值分解来解决问题：

\[{{\bf{C}}_{\bf{X}}} = {\bf{U\Sigma V}} \]

通过奇异值分解后获得的 \(\bf{\Sigma}\) 中最大的 \(K\) 个奇异值及其对应的 \(\bf{U}\) 中的 \(K\) 个列向量组成的矩阵就是变换矩阵。若是用右奇异矩阵 \(\bf{U}\) 的话就是对 \(\bf{X}\) 的列进行降维了。

求协方差矩阵 \({{\bf{C}}_{\bf{X}}} = {{\bf{X}}}{\bf{X}^{\rm T}}\) (这里其实是散度矩阵，是协方差矩阵的 \(N-1\) 倍)的复杂度为 \(O\left( {M{N^2}} \right)\)，当数据维度过大时，仍是很耗时的，此时用奇异值分解有必定的优点，好比有一些SVD的实现算法能够先不求出协方差矩阵也能求出右奇异矩阵 \(\bf{V}\)，也就是说，PCA算法能够不用作特征分解而是经过SVD来完成，这个方法在样本量很大的时候颇有效。

至于为何协方差矩阵通过特征值（奇异值）分解后前 \(K\) 个最大特征值对应的份量就是方差最大的 \(K\) 维空间呢？这里给出一个简单的理解方法，数学证实就不贴出来了。\(K\) 个特征值与\(K\) 个特征向量构成了一个 \(K\) 维空间，而特征值正是衡量每一维空间所占的份量的，特征值越大，说明这一维度所包含的信息量越大，权重就越大。前 \(K\) 个最大特征值组成的空间固然是权重最大的，从而矩阵 \(\bf{X}\) 投影到这 \(K\) 维空间中足够分散，也就是说方差足够大。

2.3 PCA的理论推导

这里引用一下

这里是引用主成分分析（PCA）原理总结这篇博客的证实。

2.4 PCA实例求解

仍是以某地房子面积与总价的问题做为例子来总结PCA算法步骤。

步骤1. 数据去中心化：

步骤2：计算散度矩阵(协方差矩阵)

步骤3：对协方差矩阵进行特征值（奇异值分解）：

步骤4：取最大的 \(K\) 个特征值与特征向量构成一个 变换矩阵 \(\bf{P}_K\)（这里 \(K=1\)）：

步骤5：将原始数据投影到 \(K\) 维空间：

\[{\bf{Y}} = {\bf{P}}_K^{\rm T}{\bf{X}} \]

咱们能够看到最后结果，那条投影线跟咱们直观感觉的同样：

3. PCA的应用

PCA在特征工程、计算机视觉、人脸识别、去噪等领域有很普遍的应用。这里仅仅测试一下PCA在图像压缩上的效果。matlab代码以下：

function [P,r] = fun_pca(X,K)
% 输入：
% X：通过中心化处理的数据
% K: 须要下降到的维度
% 输出：
% P: 投影矩阵
% r:前K个奇异值占全部奇异值的比

Cov = X*X';
[u,s] = svd(Cov) ;
P = u(:,1:K) ;
s = diag(s) ;
r = sum(abs(s(1:K)))/sum(abs(s)) ;
end

PCA for image compression

clear
img = imread('xinhuan.jpg');
img = im2double(img);
 
K=50;
% 保留主成分个数
img_matrix = [img(:,:,1),img(:,:,2),img(:,:,3)] ;
% 将RGB三通道图像合并成一个矩阵
img_matrix_mean = img_matrix - mean(img_matrix) ;
% 求均值
Proj = fun_pca(img_matrix_mean,K) ;
% 用PCA算法得到投影矩阵P
img_K = Proj'*img_matrix_mean ;
% 压缩到K维后的图像
    
%% 图像重建
img_rec = Proj*img_K ;
% 图像解压缩
img_rec = img_rec + mean(img_matrix) ;
% 重建后的图像加上均值
c_n = size(img_matrix,2)/3 ;
img_rec = cat(3,img_rec(:,1:c_n),img_rec(:,c_n+1:2*c_n),img_rec(:,2*c_n+1:end)) ;

figure(1)
subplot(1,2,1)
imshow(img),title('原图')
subplot(1,2,2)
imshow(img_rec),title(strcat('压缩后的图像(', '奇异值数量：10)'))

SVD for image compression

clear
img = imread('xinhuan.jpg');
img = im2double(img);
K = 10 ;
for i = 1:3
    [u,s,v] = svd(img(:,:,i)) ;
    s = diag(s) ;
    r = sum(abs(s(1:K)))/sum(abs(s)) ;
    s = diag(s(1:K)) ;
    W(:,:,i) = u(:,1:K)*s*v(:,1:K)' ;
end
figure(1)
subplot(2,2,1)
imshow(img),title('原图')
subplot(2,2,2)
imshow(P),title(strcat('压缩后的图像(', '奇异值数量：10)'))

奇异值分解图像压缩：

图5. SVD重构的图像与原图对比

PCA图像压缩：

图6. PCA重构的图像与原图对比

从图5和图6基本看不出来PCA和SVD压缩图像有多大区别，其特征值贡献比差别很大确不能说明什么问题，图像重建质量可使用峰值信噪比（PSNR）来评价，感兴趣的能够去试一试。

4. 讨论

PCA的主要优势有：

　　　　1）仅仅须要以方差衡量信息量，不受数据集之外的因素影响。　

　　　　2）各主成分之间正交，可消除原始数据成分间的相互影响的因素。

　　　　3）计算方法简单，主要运算是特征值分解，易于实现。

　　　　PCA算法的主要缺点有：

　　　　1）主成分各个特征维度的含义具备必定的模糊性，不如原始样本特征的解释性强。

　　　　2）方差小的非主成分也可能含有对样本差别的重要信息，因降维丢弃可能对后续数据处理有影响。

更多精彩内容请关注订阅号优化与算法和加入QQ讨论群1032493483获取更多资料