再谈协方差矩阵之主成分分析

http://pinkyjie.com/2011/02/24/covariance-pca/
算法

自从上次谈了协方差矩阵以后,感受写这种科普性文章还不错,那我就再谈一把协方差矩阵吧。上次那篇文章在理论层次介绍了下协方差矩阵,没准不少人以为这东西用处不大,其实协方差矩阵在好多学科里都有很重要的做用,好比多维的正态分布,再好比今天咱们今天的主角——主成分分析(Principal Component Analysis,简称PCA)。结合PCA相信能对协方差矩阵有个更深刻的认识~数据结构

PCA的缘起ide

PCA大概是198x年提出来的吧,简单的说,它是一种通用的降维工具。在咱们处理高维数据的时候,为了能下降后续计算的复杂度,在“预处理”阶段一般要先对原始数据进行降维,而PCA就是干这个事的。本质上讲,PCA就是将高维的数据经过线性变换投影到低维空间上去,但这个投影可不是随便投投,要遵循一个指导思想,那就是:找出最可以表明原始数据的投影方法。这里怎么理解这个思想呢?“最能表明原始数据”但愿降维后的数据不能失真,也就是说,被PCA降掉的那些维度只能是那些噪声或是冗余的数据。这里的噪声和冗余我认为能够这样认识:wordpress

  • 噪声:咱们常说“噪音污染”,意思就是“噪声”干扰咱们想听到的真正声音。一样,假设样本中某个主要的维度A,它能表明原始数据,是“咱们真正想听到的东西”,它自己含有的“能量”(即该维度的方差,为啥?别急,后文该解释的时候就有啦~)原本应该是很大的,但因为它与其余维度有那么一些千丝万缕的相关性,受到这些个相关维度的干扰,它的能量被削弱了,咱们就但愿经过PCA处理后,使维度A与其余维度的相关性尽量减弱,进而恢复维度A应有的能量,让咱们“听的更清楚”!
  • 冗余:冗余也就是多余的意思,就是有它没它都同样,放着就是占地方。一样,假如样本中有些个维度,在全部的样本上变化不明显(极端状况:在全部的样本中该维度都等于同一个数),也就是说该维度上的方差接近于零,那么显然它对区分不一样的样本丝毫起不到任何做用,这个维度便是冗余的,有它没它一个样,因此PCA应该去掉这些维度。

这么一分析,那么PCA的最终目的就是“降噪”和消灭这些“冗余”的维度,以使下降维度的同时保存数据原有的特征不失真。后面咱们将结合例子继续讨论。函数

协方差矩阵——PCA实现的关键工具

前面咱们说了,PCA的目的就是“降噪”和“去冗余”。“降噪”的目的就是使保留下来的维度间的相关性尽量小,而“去冗余”的目的就是使保留下来的维度含有的“能量”即方差尽量大。那首先的首先,咱们得须要知道各维度间的相关性以及个维度上的方差啊!那有什么数据结构能同时表现不一样维度间的相关性以及各个维度上的方差呢?天然是非协方差矩阵莫属。回忆下《浅谈协方差矩阵》的内容,协方差矩阵度量的是维度与维度之间的关系,而非样本与样本之间。协方差矩阵的主对角线上的元素是各个维度上的方差(即能量),其余元素是两两维度间的协方差(即相关性)。咱们要的东西协方差矩阵都有了,先来看“降噪”,让保留下的不一样维度间的相关性尽量小,也就是说让协方差矩阵中非对角线元素都基本为零。达到这个目的的方式天然不用说,线代中讲的很明确——矩阵对角化。而对角化后获得的矩阵,其对角线上是协方差矩阵的特征值,它还有两个身份:首先,它仍是各个维度上的新方差;其次,它是各个维度自己应该拥有的能量(能量的概念伴随特征值而来)。这也就是咱们为什么在前面称“方差”为“能量”的缘由。也许第二点可能存在疑问,但咱们应该注意到这个事实,经过对角化后,剩余维度间的相关性已经减到最弱,已经不会再受“噪声”的影响了,故此时拥有的能量应该比先前大了。看完了“降噪”,咱们的“去冗余”还没完呢。对角化后的协方差矩阵,对角线上较小的新方差对应的就是那些该去掉的维度。因此咱们只取那些含有较大能量(特征值)的维度,其他的就舍掉便可。PCA的本质其实就是对角化协方差矩阵。spa

下面就让咱们跟着上面的感受来推推公式吧。假设咱们有一个样本集X,里面有N个样本,每一个样本的维度为d。即:code

X=\{X_1,\ldots,X_N\}  \quad X_i=(x_{i1},\ldots,x_{id})\in\mathcal{R}^d, i=1,\ldots,N

将这些样本组织成样本矩阵的形式,即每行为一个样本,每一列为一个维度,获得样本矩阵S:S\in\mathcal{R}^{N\times d}。咱们先将样本进行中心化,即保证每一个维度的均值为零,只需让矩阵的每一列除以减去对应的均值便可。不少算法都会先将样本中心化,以保证全部维度上的偏移都是以零为基点的。而后,对样本矩阵计算其协方差矩阵,按照《浅谈协方差矩阵》里末尾的update,咱们知道,协方差矩阵能够简单的按下式计算获得:ip

C=\frac{S^T S}{N-1} \quad C\in\mathcal{R}^{d\times d}

下面,根据咱们上文的推理,将协方差矩阵C对角化。注意到,这里的矩阵C是是对称矩阵,对称矩阵对角化就是找到一个正交矩阵P,知足:P^TCP=\Lambda。具体操做是:先对C进行特征值分解,获得特征值矩阵(对角阵)即为\Lambda,获得特征向量矩阵并正交化即为P。显然,P,\Lambda\in\mathcal{R}^{d\times d}。假如咱们取最大的前p(p<d)个特征值对应的维度,那么这个p个特征值组成了新的对角阵\Lambda_1\in\mathcal{R}^{p\times p},对应的p个特征向量组成了新的特征向量矩阵P_1\in\mathcal{R}^{d\times p}ci

实际上,这个新的特征向量矩阵P_1就是投影矩阵,为何这么说呢?假设PCA降维后的样本矩阵为S_1,显然,根据PCA的目的,S_1中的各个维度间的协方差基本为零,也就是说,S_1的协方差矩阵应该为\Lambda_1。即知足:

\frac{S_1^TS_1}{N-1}=\Lambda_1

而咱们又有公式:

P^TCP=\Lambda \Rightarrow P_1^TCP_1=\Lambda_1

代入可得:

\frac{S_1^TS_1}{N-1}=\Lambda_1=P_1^TCP_1=P_1^T\left(\frac{S^TS}{N-1}\right)P_1=\frac{(SP_1)^T(SP_1)}{N-1}

\Rightarrow S_1=SP_1 \quad S_1\in\mathcal{R}^{N\times p}

因为样本矩阵S_{N\times d}的每一行是一个样本,特征向量矩阵P_{1(d\times p)}的每一列是一个特征向量。右乘P_1至关于每一个样本以P_1的特征向量为基进行线性变换,获得的新样本矩阵S_1\in\mathcal{R}^{N\times p}中每一个样本的维数变为了p,完成了降维操做。实际上,P_1中的特征向量就是低维空间新的坐标系,称之为“主成分”。这就是“主成分分析”的名称由来。同时,S_1的协方差矩阵\Lambda_1为近对角阵,说明不一样维度间已经基本独立,噪声和冗余的数据已经不见了。至此,整个PCA的过程已经结束,小小总结一下:

  1. 造成样本矩阵,样本中心化
  2. 计算样本矩阵的协方差矩阵
  3. 对协方差矩阵进行特征值分解,选取最大的p个特征值对应的特征向量组成投影矩阵
  4. 对原始样本矩阵进行投影,获得降维后的新样本矩阵

Matlab中PCA实战

首先,随机产生一个10*3维的整数矩阵做为样本集,10为样本的个数,3为样本的维数。

1
S = fix(rand(10,3)*50);

计算协方差矩阵:

1
2
3
4
S = S - repmat(mean(S),10,1);
C = (S'*S)./(size(S,1)-1);
or
C = cov(S);

对协方差矩阵进行特征值分解:

1
[P,Lambda] = eig(C);

这里因为三个方差没有明显特别小的,因此咱们都保留下来,虽然维度没有降,但观察Lambda(即PCA后的样本协方差矩阵)和C(即原始的协方差矩阵),能够发现,3个维度上的方差都有增大,也就是能量都比原来增大了,3个维度上的方差有所变化,但对角线之和没有变,能量从新获得了分配,这就是“降噪”的功劳。最后咱们获得降维后的样本矩阵:

1
S1 = S*P;

为了验证,咱们调用matlab自带的主成分分析函数princomp:

1
[COEFF,SCORE] = princomp(S) % COEFF表示投影矩阵,SCORE表示投影后新样本矩阵

对比,能够发现,SCORE和S1在不考虑维度顺序和正负的状况下是彻底吻合的,之因此咱们计算的S1的维度顺序不一样,是由于一般都是将投影矩阵P按能量(特征值)的降序排列的,而刚才咱们用eig函数获得的结果是升序。另外,在一般的应用中,咱们通常是不使用matlab的princomp函数的,由于它不能真正的降维(不提供相关参数,仍是我没发现?)。通常状况下,咱们都是按照协方差矩阵分解后特征值所包含的能量来算的,好比取90%的能量,那就从最大的特征值开始加,一直到部分和占特征值总和的90%为止,此时部分和含有的特征值个数即为p。

通过了一番推公式加敲代码的过程,相信你们对主成分分析应该不陌生了吧,同时对协方差矩阵也有了更深层次的认识了吧,它可不仅是花花枪啊。我我的以为PCA在数学上的理论仍是很完备的,想必这也是它能在多种应用中博得鳌头的缘由吧。