使用SVD奇异值分解求解PCA+Python实现

这几天在看有关ＰＣＡ的博客时，感觉文章中针对如何用ＳＶＤ解ＰＣＡ的过程的讲解不是很清晰，自己捋了捋思路，把自己对于ＳＶＤ解ＰＣＡ的步骤分享出来

关于ＰＣＡ的思想和原理，这里不阐述，推荐这篇文章，当初就是看这篇文章入坑ＰＣＡ的
http://blog.codinglabs.org/articles/pca-tutorial.html

1.直接解PCA

设有n组m维数据，组成 m*n 的矩阵 ${\mathrm{X}}_{m\ast n}$ , 则 在去中心化后，$\mathrm{X}$ 的协方差矩阵
　　　　　　　　　　　 $\mathcal{C}$ = $\frac{1}{n}$ $\mathrm{X}$${\mathrm{X}}^T$
目标：将 $\mathcal{C}$ 对角化
如何对角化 $\mathcal{C}$：
设$\mathcal{P}$是正交基按行组成的矩阵(最终希望得到的转换矩阵，$\mathcal{P}$的前$k$行就是要寻找的基，可对$\mathrm{X}$降维)，则$\mathcal{Y}$=$\mathcal{P}$$\mathrm{X}$中，$\mathcal{Y}$是$\mathrm{X}$对$\mathcal{P}$做基变换得到的结果
则$\mathcal{Y}$ 的协方差矩阵$\mathcal{D}$ = $\frac{1}{n}$$\mathcal{Y}$${\mathcal{Y}}^T$ = $\mathcal{P}$$\mathcal{C}$${\mathcal{P}}^T$ 应该是一个对角阵。

已知知道一个$m\ast$$m$的实对称矩阵一定可以找到$m$个正交特征向量，特征向量组成的矩阵为$\mathcal{E}$

则有 　${\mathcal{E}}^T$$\mathcal{C}$$\mathcal{E}$ = $\Lambda$，对比上式，可以得到 　$\mathcal{P}$ = ${\mathcal{E}}^T$

即我们现在的目标是求出$\mathcal{C}$的特征向量矩阵，并将其转置

但是，高维数据的协方差矩阵非常庞大，比如单位样本是一副100*100的uint8图像，则一副图像就要表示为一个 10000 * 1 的列向量，当样本数量少于10000时，协方差矩阵的大小为10000 *10000(800m左右)，计算机求解这个矩阵会非常困难，我尝试过，用Matlab求解要用5分钟左右的时间，如果图像更大，比如时 70000 * 1 的情况，则协方差矩阵大小达到45G，求解成为了几乎不可能的事情

这时候我们就要借助一个工具 ＳＶＤ特征值分解，来曲线求解协方差矩阵的特征值矩阵

2. 用SVD特征值分解来解PCA

这里不阐述ＳＶＤ的原理，同样我推荐这篇文章
https://mp.weixin.qq.com/s/Dv51K8JETakIKe5dPBAPVg

$\mathrm{X}$经过分解后，可以得到　
　　　　　　　${\mathrm{X}}_{m\ast n}$ = ${\mathcal{U}}_{m\ast m}$${\Sigma }_{m\ast n}$${\mathcal{V}}_{n\ast n}^T$

其中${\mathcal{U}}_{m\ast m}$是$\mathrm{X}$${\mathrm{X}}^T$的特征向量按列组成的矩阵
${\mathcal{V}}_{n\ast n}$是${\mathrm{X}}^T$$\mathrm{X}$的特征向量按列组成的矩阵
${\Sigma }_{m\ast n}$的对角线上是奇异值$\sigma$，奇异值与特征值的关系为${\sigma }^2$ = $\lambda$

按以上的性质，可以得到 $\mathcal{P}$ = ${\mathcal{U}}^T$
所以我们现在的目标是求出$\mathcal{U}$

SVD的分解思想是
　　　　　　　${\mathrm{X}}_{m\ast n}$$\approx$${\mathcal{U}}_{m\ast r}$${\Sigma }_{r\ast r}$${\mathcal{V}}_{r\ast n}^T$

我们现在将${\mathrm{X}}_{m\ast n}$近似为
${\mathrm{X}}_{m\ast n}$$\approx$${\mathcal{U}}_{m\ast n}$${\Sigma }_{n\ast n}$${\mathcal{V}}_{n\ast n}^T$ （对于n <= m）
${\mathcal{V}}_{n\ast n}$的特征值矩阵为${\Sigma }_{n\ast n}^2$

所以我们可以求得
　　　　　　　${\mathcal{U}}_{m\ast n}$ = ${\mathrm{X}}_{m\ast n}$$({\Sigma }_{n\ast n}$${\mathcal{V}}_{n\ast n}{)}^{-1}$
从而就可以依据特征值大小将${\mathcal{U}}_{m\ast n}$的特征值排序($\mathcal{U}$和$\mathcal{V}$拥有相同的特征值)，就可以降维至${\mathcal{U}}_{m\ast k}$，再将其转置就可以得到变换矩阵$\mathcal{P}$

这时可能会有疑问
${\mathcal{U}}_{m\ast m}$和${\mathcal{V}}_{n\ast n}$维数不一样，为什么特征值个数是min{m, n}，其实，若直接对前者进行特征值求解，排序后，会发现，前者只有前n个特征值是实数，后面都是虚数

可以看到，在用ＳＶＤ求解ＰＣＡ时，我们只需要计算${\mathcal{V}}_{n\ast n}$的特征值和特征向量，通常情况下，数据样本 n 是远小于数据维数 m 的，因此计算 n*n 矩阵的特征值和特征向量，将可以大大减小计算量，ＳＶＤ“曲线救国”的目的也就达到了

3.Python实现PCA

在实际进行PCA的求解中，需要根据数据的维度与样本数来进行直接求解PCA与用SVD求解PCA的权衡

import numpy as np

class PCA():
    """
    PCA降维
    """
    def __init__(self):
        self.X = None
        self.k = None

    def getTransMat(self, X, k, svd):
        """
        训练得到特征空间中的主成分转换矩阵
        args:
        X -- 训练数据 shape=(features, samples)
        k -- 要降维的维度，若小于0，则自动选择, 选择阈值为特征值总和95%
        svd -- 是否使用svd求解PCA, true or false 如果训练样本数多于特征数，选择直接解更省内存(faster)，反之选择svd奇异值分解更好
        Returns:
        transMat -- 转换矩阵 shape=(降维后维度, 降维前维度)
        """
        # 将X的每一维进行零均值化
        self.X = X
        self.k = k
        m = np.mean(self.X, 1)
        m = np.array([m])
        self.X = self.X - m.T

        if svd:
            # 用SVD特征值分解解出变换矩阵
            # 求出 X' * X 的特征值矩阵V，并求出奇异值对角阵Σ（奇异值从大到小排序）
            C = np.dot(self.X.T, self.X)

            # 求出 X' * X 的特征值矩阵及对应的特征向量和奇异值矩阵
            [eigValue, V] = np.linalg.eig(C)
            eigValue = abs(np.array([eigValue]))
            sinValue = np.sqrt(eigValue)
            temp = np.zeros([np.size(sinValue, 1), np.size(sinValue, 1)])
            for i in range(0, np.size(sinValue, 1)):
                temp[i, i] = sinValue[0, i]
            sinValue = temp

            # 将特征向量归一化
            V = V / np.linalg.norm(V, axis=0)

            # 求出变换矩阵
            U = np.dot(self.X, np.linalg.inv(np.dot(sinValue, V.T)))
            U = U.T
        
        else:
            # 直接解PCA
            C = np.dot(self.X, self.X.T)
            [eigValue, V] = np.linalg.eig(C)
            eigValue = np.array([eigValue])

            # 将特征向量归一化
            V = V / np.linalg.norm(V, axis=0)

            # 求出变换矩阵
            U = V.T

        # 按特征值大小给变换矩阵行向量重新排序
        index = np.argsort(-eigValue)
        eigValue.tolist()[0].reverse()
        temp = U.copy()
        for i in range(0, np.size(U, 0)):
            U[i, :] = temp[index[0][i], :]

        # 自动选择降维维度
        self.autoChooseDim(eigValue)

        # 获得降维变换矩阵
        transMat = U[0:self.k, :]
        return np.real(transMat)
        
    def autoChooseDim(self, sortV):
        """
        自适应选择降维维度
        args:
        sortV -- 已经排好序的特征值序列
        returns:
        None
        """
        if self.k <= 0:
            threshold = sum(sum(sortV)) * 0.95
            total = 0
            for i in range(0, len(sortV[0])):
                total += sortV[0][i]
                if total > threshold:
                    break
            self.k = i
        else:
            pass

if __name__ == "__main__":
    pass