PCA主成分分析

PCA主成分分析：

主成分分析是一种矩阵的压缩算法，在减少矩阵维数的同时尽可能的保留原矩阵的信息，简单来说就是将 n×m的矩阵转换成n×k的矩阵（通常k是小于m的），仅保留矩阵中所存在的主要特性，从而可以大大节省空间和数据量。 然后假期学习到图像主成分分析的书上直接给出了程序实现的过程，而其中原理并不是特别清楚。去网上找别人写的博客去看的时候，感觉都写的过于片面。然后我就想在这些基础之上总结一下。

PCA数学原理分析：

PCA数学原理分析

分析这篇文章写的很好：可以参考这篇文章。

我们生活所面对的数据都可以完全被抽象为一组向量。

内积和 投影：

两个维数相同的向量的内积被定义为：
$$(a_1,a_2...a_n)\ast (b_1,b_2...b_n)=a_1{b}_1+a_2{b}_2+...+a_n{b}_n$$

设A=(x1,y1),B=(x2,y2).好，现在我们从A点向B所在直线引一条垂线。我们知道垂线与B的交点叫做A在B上的投影，再设A与B的夹角是a，则投影的矢量长度为
$$|A|cosa$$
接下来内积还有另外一种表示方法：
$$A\ast B=|A|\ast |B|\ast cosa$$
也就是说，设向量B的模为1，则A与B的内积值等于A向B所在直线投影的矢量长度！

基：

这个红色线段的向量我们通常都直接认为是(3,2).但实际上这样说只是相对于当前的这个坐标系来说。所以来说并不是特别准确。

实际要准确表示这个向量的话，应该按照下面的去表示：
$$3(1,0{)}^T+2(0,1{)}^T$$

此处(1,0)和(0,1)叫做二维空间中的一组基。所以，要准确描述向量，首先要确定一组基，然后给出在基所在的各个直线上的投影值，就可以了。只不过我们经常省略第一步，而默认以(1,0)和(0,1)为基。我们之所以默认选择(1,0)和(0,1)为基，当然是比较方便，因为它们分别是x和y轴正方向上的单位向量，因此就使得二维平面上点坐标和向量一一对应，非常方便。但实际上任何两个线性无关的二维向量都可以成为一组基，所谓线性无关在二维平面内可以直观认为是两个不在一条直线上的向量。

，(1,1)和(-1,1)也可以成为一组基。一般来说，我们希望基的模是1，因为从内积的意义可以看到，如果基的模是1，那么就可以方便的用向量点乘基而直接获得其在新基上的坐标了！实际上，对应任何一个向量我们总可以找到其同方向上模为1的向量，只要让两个分量分别除以模就好了。这一步通常在线性代数也叫做单位化。所以上面的基可以变为：
$$(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})和(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}）$$
现在我们要想获得(3,2)在新的基上的坐标，我们只需要分别计算(3,2)与新的基的内积即可，

所以可以得到：
$$(\frac{5}{\sqrt{2}},\frac{-1}{\sqrt{2}})$$

我们所选的大多数情况下选的基都是正交的。

基变换的矩阵表示

因为进行基变换的方法坐标分别与基做内积运算。所以我们可以用矩阵去表示这种变换，我们写的时候通常将原来的坐标写为列向量的形式。
$$\left(\begin{array}{cc}1/\sqrt{2}& 1/\sqrt{2}\\ -1/\sqrt{2}& 1/\sqrt{2}\end{array}\right)\left(\begin{array}{c}3\\ 2\end{array}\right)=\left(\begin{array}{c}5/\sqrt{2}\\ -1/\sqrt{2}\end{array}\right)$$
如果我们有m个二维向量，只要将二维向量按列排成一个两行m列矩阵，然后用“基矩阵”乘以这个矩阵，就得到了所有这些向量在新基下的值。

一般的，如果我们有M个N维向量，想将其变换为由R个N维向量表示的新空间中，那么首先将R个基按行组成矩阵A，然后将向量按列组成矩阵B，那么两矩阵的乘积AB就是变换结果，其中AB的第m列为A中第m列变换后的结果。

矩阵的相乘在看完那篇博客后好像明白他们的意义。其实矩阵变换无非就是两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。

协方差矩阵及其表示：

如果进行矩阵变换的时候基的数量少于向量本身的维数的时候我们就达到了降维的目的。但是如何选取基又是一个问题。看了别人的博客写的后感觉这种变换的思想真的巧妙。

选取一个合适的基的意思是变换后的坐标尽量分散。而分散的程度我们在数学上通常使用方差来表示，一个字段的方差可以看做是每个元素与字段均值的差的平方和的均值，
$$Var(a)=\frac{1}{m}\sum _{i=1}^m(a_i-\mu {)}^2$$
在这里为了简化运算我们通常提前对数据进行处理，每个数据减去这个字段的均值。

处理过程很简单。举个例子：
$$1,2,3,4,5这五个数据的均值为15/5=3所以每个数据都减去3.所以处理后的数据为：-2，-1,0，1,2.这五个数据均值一定为0.$$
这样处理后，
$$这样\mu 可以为0了。$$
所以寻找一个基，使得所有数据变换为这个基上的坐标表示后，方差值最大。

对于上面二维降成一维的问题来说，找到那个使得方差最大的方向就可以了。不过对于更高维，还有一个问题需要解决。考虑三维降到二维问题。与之前相同，首先我们希望找到一个方向使得投影后方差最大，这样就完成了第一个方向的选择，继而我们选择第二个投影方向。

如果我们还是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一起”，显然这样的维度是没有用的，因此，应该有其他约束条件。从直观上说，让两个字段尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个字段不是完全独立，必然存在重复表示的信息。

数学上可以用两个字段的协方差表示其相关性，由于已经让每个字段均值为0，协方差原始公式为：

则协方差的公式如下：
$$Cov(a,b)=\frac{1}{m}\sum _{i=1}^m{a}_i{b}_i$$
可以看到，在字段均值为0的情况下，两个字段的协方差简洁的表示为其内积除以元素数m。当协方差为0时，表示两个字段完全独立。为了让协方差为0，我们选择第二个基时只能在与第一个基正交的方向上选择。因此最终选择的两个方向一定是正交的。为什么正交的话？因为线性代数里面也说正交的话内积为0，所以协方差就为0.

将一组N维向量降为K维（K大于0，小于N），其目标是选择K个单位（模为1）正交基，使得原始数据变换到这组基上后，各字段两两间协方差为0，而字段的方差则尽可能大（在正交的约束下，取最大的K个方差）

上面就是我们要得到的目标，但是就是上面的这个结论并不能告诉我们要具体怎么去做。所以说数学的思想真的巧妙。

我们设有a和b俩个字段。
$$X=\left(\begin{array}{cccc}{a}_1& a_2& \cdots & a_m\\ b_1& b_2& \cdots & b_m\end{array}\right)$$
然后我们用X乘以X的转置，并乘上系数1/m：
$$\frac{1}{m}X{X}^{\mathsf{T}}=\left(\begin{array}{cc}\frac{1}{m}{\sum }_{i=1}^m{a}_i^2& \frac{1}{m}{\sum }_{i=1}^m{a}_i{b}_i\\ \frac{1}{m}{\sum }_{i=1}^m{a}_i{b}_i& \frac{1}{m}{\sum }_{i=1}^m{b}_i^2\end{array}\right)$$
然后神奇的事情发生了，主对角线上是方差，斜对角线是协方差。

设我们有m个n维数据记录，将其按列排成n乘m的矩阵X，设
$$C=\frac{1}{m}X{X}^{\mathsf{T}}$$
，则C是一个对称矩阵，其对角线分别个各个字段的方差，而第i行j列和j行i列元素相同，表示i和j两个字段的协方差

协方差矩阵对角化：

为什么要将协方差矩阵对角化呢。首先我们了解一下相似矩阵。设A和B为n阶矩阵如果存在n阶可逆矩阵P，使得：
$$P^{-1}AP=B$$
我们就称A相似于B，或者说A和B相似。矩阵相似具有非常多好用的性质：

A和B有相等的行列式。A和B有相等的迹。

相似矩阵具有许多共同的性质，对于N阶矩阵。我们希望与A相似的矩阵中寻找一个较简单的矩阵。所以我们考虑n阶矩阵是否与一个对角矩阵相似的问题的时候。这时候就用到了矩阵对角化。

设原始数据矩阵X对应的协方差矩阵为C，而P是一组基按行组成的矩阵，设Y=PX，则Y为X对P做基变换后的数据。设Y的协方差矩阵为D，我们推导一下D与C的关系：
$$\begin{array}{lll}D& =& \frac{1}{m}Y{Y}^{\mathsf{T}}\\ & =& \frac{1}{m}(PX)(PX{)}^{\mathsf{T}}\\ & =& \frac{1}{m}PX{X}^{\mathsf{T}}{P}^{\mathsf{T}}\\ & =& P(\frac{1}{m}X{X}^{\mathsf{T}})P^{\mathsf{T}}\\ & =& PC{P}^{\mathsf{T}}\end{array}$$
现在事情很明白了！我们要找的P不是别的，而是能让原始协方差矩阵对角化的P。换句话说，优化目标变成了寻找一个矩阵P，满足
$$PC{P}^{\mathsf{T}}$$

---

是一个对角矩阵，并且对角元素按从大到小依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件。

，在线性代数上，实对称矩阵有一系列非常好的性质：

1）实对称矩阵不同特征值对应的特征向量必然正交。
$$
2）设特征向量\lambda重数为r，则必然存在r个线性无关的特征向量对应于\lambda\

因此可以将这r个特征向量单位正交化。
$$

一个n行n列的实对称矩阵一定可以找到n个单位正交特征向量，设这n个特征向量为:
$$e_1,e_2,\cdots ,e_n$$
则对协方差矩阵C有如下结论：
$$E^{\mathsf{T}}CE=\Lambda =\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_n\end{array}\right)$$

$$P=E^{\mathsf{T}}$$

矩阵P就已经得到了。

P是协方差矩阵的特征向量单位化后按行排列出的矩阵，其中每一行都是C的一个特征向量。如果设P按照特征值的从大到小，将特征向量从上到下排列，则用P的前K行组成的矩阵乘以原始数据矩阵X，就得到了我们需要的降维后的数据矩阵Y。

PCA算法

总结一下PCA的算法步骤：

设有m条n维数据。

1）将原始数据按列组成n行m列矩阵X

2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值

3）求出协方差矩阵
$$C=\frac{1}{m}X{X}^{\mathsf{T}}$$
4）求出协方差矩阵的特征值及对应的特征向量

5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P

6）Y=PX即为降维到k维后的数据

代码如下：

from PIL import Image
from numpy import *
from pylab import *
from imtools import get_imlist
import pickle
def pca(X):
    """主成分分析
    输入：矩阵X其中该矩阵中存储训练数据，每一行为一条训练数据
    返回：投影矩阵（按照维度的重要性排序），方差和均值
    """
    num_data,dim = X.shape#获取维数
    #数据中心化
    mean_x = X.mean(axis=0)#axis = 0：压缩行，对各列求均值，返回 1* n 矩阵axis =1 ：压缩列，对各行求均值，返回 m *1 矩阵
    X = X - mean_x
    if dim > num_data:
        M = dot(X,X.T) #PCA-使用紧致技巧得到协方差矩阵
        e,EV = linalg.eigh(M)#计算特征值和特征向量
        tmp = dot(X.T,EV).T#这就是紧致技巧
        V  = tmp[::-1]#由于最后的特征向量就是我们需要的，所以将其逆转
        S = sqrt(e)[::-1]#由于特征值是按照递增顺序排列的，所以需要将其逆转。
        for i in range(V.shape[1]):
            V[:,i] /= S
    else:
        U,S,V = linalg.svd(X)
        V = V[:num_data]#仅仅返回num_data维的数据才合理
    return V,S,mean_x#返回投影矩阵，方差和均值
imlist = get_imlist(r'..\image\1.jpg','.jpg')
im = array(Image.open(imlist[0]))#打开一幅图像，获得其大小
m,n = im.shape[0:2]
imnum = len(imlist)
#创建矩阵，保存所有压平后的图像数据
immatrix = array([array(Image.open(im)).flatten() for im in imlist],'f')
#执行PCA操作
V,S,immean = pca(immatrix)
#显示一些图像（均值图像和前10个模式）
figure()
gray()
subplot(2,5,1)#m和n代表在一个图像窗口中显示m行n列个图像
imshow(immean.reshape(m,n))
for i in range(9):
    subplot(2,5,i+2)
    imshow(V[i].reshape(m,n))