这几天在看有关ＰＣＡ的博客时，感觉文章中针对如何用ＳＶＤ解ＰＣＡ的过程的讲解不是很清晰，自己捋了捋思路，把自己对于ＳＶＤ解ＰＣＡ的步骤分享出来
关于ＰＣＡ的思想和原理，这里不阐述，推荐这篇文章，当初就是看这篇文章入坑ＰＣＡ的
/articles/pca-tutorial.html
1.直接解PCA
设有n组m维数据，组成 m*n 的矩阵 Xm∗n\Chi_{m*n}Xm∗n , 则 在去中心化后，X\ChiX 的协方差矩阵
C\mathcal{C}C = 1n\frac{1}{n}n1 X\ChiXXT\Chi^{T}XT
目标：将 C\mathcal{C}C 对角化
如何对角化 C\mathcal{C}C：
设P\mathcal{P}P是正交基按行组成的矩阵(最终希望得到的转换矩阵，P\mathcal{P}P的前k\mathcal{k}k行就是要寻找的基，可对X\ChiX降维)，则Y\mathcal{Y}Y=P\mathcal{P}PX\ChiX中，Y\mathcal{Y}Y是X\ChiX对P\mathcal{P}P做基变换得到的结果
则Y\mathcal{Y}Y 的协方差矩阵D\mathcal{D}D = 1n\frac{1}{n}n1Y\mathcal{Y}YYT\mathcal{Y}^{T}YT = P\mathcal{P}PC\mathcal{C}CPT\mathcal{P}^{T}PT 应该是一个对角阵。
已知知道一个m∗\mathcal{m}*m∗m\mathcal{m}m的实对称矩阵一定可以找到m\mathcal{m}m个正交特征向量，特征向量组成的矩阵为E\mathcal{E}E
则有 ET\mathcal{E}^{T}ETC\mathcal{C}CE\mathcal{E}E = Λ\LambdaΛ，对比上式，可以得到 P\mathcal{P}P = ET\mathcal{E}^{T}ET
即我们现在的目标是求出C\mathcal{C}C的特征向量矩阵，并将其转置
但是，高维数据的协方差矩阵非常庞大，比如单位样本是一副100*100的uint8图像，则一副图像就要表示为一个 10000 * 1 的列向量，当样本数量少于10000时，协方差矩阵的大小为10000 *10000(800m左右)，计算机求解这个矩阵会非常困难，我尝试过，用Matlab求解要用5分钟左右的时间，如果图像更大，比如时 70000 * 1 的情况，则协方差矩阵大小达到45G，求解成为了几乎不可能的事情
这时候我们就要借助一个工具 ＳＶＤ特征值分解，来曲线求解协方差矩阵的特征值矩阵
2. 用SVD特征值分解来解PCA
这里不阐述ＳＶＤ的原理，同样我推荐这篇文章
https://mp./s/Dv51K8JETakIKe5dPBAPVg
X\ChiX经过分解后，可以得到
Xm∗n\Chi_{m*n}Xm∗n = Um∗m\mathcal{U}_{m*m}Um∗mΣm∗n\Sigma_{m*n}Σm∗nVn∗nT\mathcal{V}^{T}_{n*n}Vn∗nT
其中Um∗m\mathcal{U}_{m*m}Um∗m是X\ChiXXT\Chi^{T}XT的特征向量按列组成的矩阵
Vn∗n\mathcal{V}_{n*n}Vn∗n是XT\Chi^{T}XTX\ChiX的特征向量按列组成的矩阵
Σm∗n\Sigma_{m*n}Σm∗n的对角线上是奇异值σ\sigmaσ，奇异值与特征值的关系为σ2\sigma^{2}σ2 = λ\lambdaλ
按以上的性质，可以得到 P\mathcal{P}P = UT\mathcal{U}^{T}UT
所以我们现在的目标是求出U\mathcal{U}U
SVD的分解思想是
Xm∗n\Chi_{m*n}Xm∗n≈\approx≈Um∗r\mathcal{U}_{m*r}Um∗rΣr∗r\Sigma_{r*r}Σr∗rVr∗nT\mathcal{V}^{T}_{r*n}Vr∗nT
我们现在将Xm∗n\Chi_{m*n}Xm∗n近似为
Xm∗n\Chi_{m*n}Xm∗n≈\approx≈Um∗n\mathcal{U}_{m*n}Um∗nΣn∗n\Sigma_{n*n}Σn∗nVn∗nT\mathcal{V}^{T}_{n*n}Vn∗nT （对于n <= m）
Vn∗n\mathcal{V}_{n*n}Vn∗n的特征值矩阵为Σn∗n2\Sigma_{n*n}^{2}Σn∗n2
所以我们可以求得
Um∗n\mathcal{U}_{m*n}Um∗n = Xm∗n\Chi_{m*n}Xm∗n(Σn∗n(\Sigma_{n*n}(Σn∗nVn∗n)−1\mathcal{V}_{n*n})^{-1}Vn∗n)−1
从而就可以依据特征值大小将Um∗n\mathcal{U}_{m*n}Um∗n的特征值排序(U\mathcal{U}U和V\mathcal{V}V拥有相同的特征值)，就可以降维至Um∗k\mathcal{U}_{m*k}Um∗k，再将其转置就可以得到变换矩阵P\mathcal{P}P
这时可能会有疑问
Um∗m\mathcal{U}_{m*m}Um∗m和Vn∗n\mathcal{V}_{n*n}Vn∗n维数不一样，为什么特征值个数是min{m, n}，其实，若直接对前者进行特征值求解，排序后，会发现，前者只有前n个特征值是实数，后面都是虚数
可以看到，在用ＳＶＤ求解ＰＣＡ时，我们只需要计算Vn∗n\mathcal{V}_{n*n}Vn∗n的特征值和特征向量，通常情况下，数据样本 n 是远小于数据维数 m 的，因此计算 n*n 矩阵的特征值和特征向量，将可以大大减小计算量，ＳＶＤ“曲线救国”的目的也就达到了
3.Python实现PCA
在实际进行PCA的求解中，需要根据数据的维度与样本数来进行直接求解PCA与用SVD求解PCA的权衡
import numpy as npclass PCA():"""PCA降维"""def __init__(self):self.X = Noneself.k = Nonedef 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 = Xself.k = km = np.mean(self.X, 1)m = np.array([m])self.X = self.X - m.Tif 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.Telse:# 直接解PCAC = 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.95total = 0for i in range(0, len(sortV[0])):total += sortV[0][i]if total > threshold:breakself.k = ielse:passif __name__ == "__main__":pass
