主成分分析(PCA)以及python实现

前置知识

• 二次型
• 对角化
• 奇异值分解
• 施密特正交化

PCA

主成分分析是一种从原始数据中消除冗余信息，只需要一个或者几个属性合成的属性就可以提供大部分信息的方法。大概过程是找出一个特殊的线性组合，给要分析的元素赋予一个权值，综合得到一个新的数据。

假设每一个样本有p个属性，一共有n个样本，那么这些向量可以构成一个$p\times n$的矩阵，称为观测矩阵。

主成分分析实际上相当于变换坐标轴，并降低维度。
如一个二维数据集如下图：

经过PCA进行降维之后可以使用一个维度进行表示，如下图

均值和协方差

均值

令$[{\mathbf{X}}_1\cdots {\mathbf{X}}_N]$为一个$p\times N$观测矩阵，则其样本平均值M为：

$\mathbf{M}=\frac{1}{N}({\mathbf{X}}_1+\cdots +{\mathbf{X}}_N)$
令 ${\hat{\mathbf{X}}}_k={\mathbf{X}}_k-\mathbf{M},k=1,2,\cdots ,N$
则新矩阵 $B=[{\hat{\mathbf{X}}}_1\cdots {\hat{\mathbf{X}}}_N]$
具有零样本均值，这种矩阵 $B$被称为 平均偏差形式

协方差矩阵

协方差矩阵是一个$p\times p$的矩阵$S$,其定义为

$S=\frac{1}{N-1}B{B}^T$

由于任何具有$B{B}^T$形式的矩阵都是半正定的，所以$S$也是半正定的。

PS： 虽然协方差矩阵是半正定矩阵，但是由于计算误差，算出来的特征值可能是负数。

其中的对角线元素$s_{jj}$称为$x_j$的方差。方差用来度量$x_j$的分散性。

数据的总方差被称为矩阵的迹，即$S$的对角线元素之和，记做$tr(S)$

$S$中的元素$s_{ij},i\ne j$称为$x_i$和$x_j$的协方差。若协方差为0，则称$x_i$和$x_j$是无关的。当协方差矩阵是对角阵或几乎是对角阵时，${\mathbf{X}}_1\cdots {\mathbf{X}}_N$中多变量分析可以简化。

主成分分析

变量代换

设矩阵$[{\mathbf{X}}_1\cdots {\mathbf{X}}_N]$已经称为平均偏差形式，主成分分析的目标是找到一个$p\times p$的正交矩阵$P=[{\mathbf{u}}_1\cdots {\mathbf{u}}_p]$,确定一个变量代换$\mathbf{X}=P\mathbf{Y}$,或
$$\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_p\end{array}\right]=[{\mathbf{u}}_1\cdots {\mathbf{u}}_{\mathbf{p}}]\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_p\end{array}\right]$$
满足新的变量$y_1,\cdots ,y_p$两两无关，并保证整理后的方差具有递减顺序。
变量正交变换$X=PY$说明，每一个观测向量${\mathbf{X}}_k$得到一个新”名称“${\mathbf{Y}}_k$,即${\mathbf{Y}}_k$是${\mathbf{X}}_k$以$P$为基的坐标。
易得对任何正交矩阵$P$，${\mathbf{Y}}_1,\cdots ,{\mathbf{Y}}_{\mathbf{N}}$的协方差是$P^{T}SP$，于是可以构造符合要求的P。设$D$是对角矩阵且$S$的特征值在$D$的对角线上，重新整理使得${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_p\ge 0$即
$$D=\left[\begin{array}{cccc}{\lambda }_1& 0& \cdots & 0\\ 0& {\lambda }_2& \cdots & 0\\ ⋮& & \ddots & ⋮\\ 0& \cdots & \cdots & {\lambda }_p\end{array}\right]$$找到正交矩阵$P$，它的列对应单位特征向量${\mathbf{u}}_1,\cdots {\mathbf{u}}_p$,且有$S=PD{P}^T$。
协方差矩阵$S$的单位特征向量${\mathbf{u}}_1,\cdots {\mathbf{u}}_p$称为(观测矩阵中)数据的主成分。第一主成分是$S$中最大特征值对应的特征向量，第二主成分以此类推。
设$c1,c2,\cdots ,c_p$为${\mathbf{u}}_1$中的元素，第一主成分${\mathbf{u}}_1$可以通过公式$y_1={\mathbf{u}}_1^T\mathbf{X}=c_1{x}_1+c_2{x}_2+\cdots +c_p{x}_p$来确定新变量$y_1$，可以看出$y_1$是原变量$x_1,\cdots ,x_p$的线性组合。通过同样的方式，可以确定$y_2$等。
通过这种方式，可以将变量进行代换，但是优点并不明显。但是如果在对角矩阵$D$中，前$p^{\prime }$个方差和比其他方差大得多，那么我们就可以近似的把数据当成$p^{\prime }$维数据而不是$p$维的。

多维数据的降维

选取一个重构阈值$t$，选取使得下式成立的最小$p^{\prime }$值

$\frac{{\sum }_{i=1}^{p^{\prime }}{\lambda }_i}{{\sum }_{i=1}^p{\lambda }_i}\ge t$
PCA仅需保留 $P$和样本均值 $M$即可将新样本投影到低维空间。转换过程不可避免的造成数据的损失，但是舍弃这部分信息往往是必要的。一方面降维之后可以使样本的采样密度增大，这也是降维的主要动机，其次可以在一定程度上起到降噪的效果，因为最小的特征值对应的特征向量往往与噪声有关。

算法步骤

input

样本集$[{\mathbf{X}}_1\cdots {\mathbf{X}}_N]$，重构阈值$t$

process

1. 进行去中心化,得到矩阵$B$
2. 计算样本的协方差矩阵$S=B{B}^T$(常用对输入矩阵的奇异值分解代替特征值分解)
3. 对协方差矩阵进行特征值分解
4. 取出最大的$p^{\prime }$个特征值对应的特征向量 ${\mathbf{u}}_1,\cdots {\mathbf{u}}_{p^{\prime }}$

Output

投影矩阵$P$ = [${\mathbf{u}}_1\cdots {\mathbf{u}}_{p^{\prime }}$]

代码实现

特征值分解

from numpy import *

def PCA( data, t ):
    m, n = data.shape
    M = data.sum(1) / n
    M = tile(M,(1,n))
    B = data - M
    D, P = linalg.eig( B*B.T /(n - 1) )
    index = argsort(-D)
    D = D[index]; P = P[index];D[D<0] = 0
    denominator = sum(D); numerator = 0
    for i in range(len(D)):
        numerator += D[i]
        if numerator / denominator >= t:
            return diag(D[:i + 1]),P[:,:i + 1], M
    return diag(D), P, M

奇异值分解

from numpy import *

def PCA( data, k ):
    m, n = data.shape
    M = data.sum(1) / n
    M = tile(M,(1,n))
    B = data - M
    U, S, V = linalg.svd( B*B.T /(n - 1) )
    return U[:,:k+1], S[:,:k+1]