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

关于PCA的思想和原理,这里不阐述,推荐这篇文章,当初就是看这篇文章入坑PCA的
http://blog.codinglabs.org/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}n1​Y\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,求解成为了几乎不可能的事情

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

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

这里不阐述SVD的原理,同样我推荐这篇文章
https://mp.weixin.qq.com/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∗n​Vn∗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∗r​Vr∗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∗n​Vn∗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∗n​Vn∗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个特征值是实数,后面都是虚数

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

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

使用SVD奇异值分解求解PCA+Python实现相关推荐

  1. SVD奇异值分解(PCA,LSI)

    版权声明: 本文由LeftNotEasy发布于http://leftnoteasy.cnblogs.com, 本文可以被全部的转载或者部分使用,但请注明出处,如果有问题,请联系wheeleast@gm ...

  2. 【数学和算法】最小二乘法,SVD奇异值分解、LU分解的应用场景

    多项式拟合曲线时使用最小二乘法,把问题化简为A*x=B的线性方程组: 然后使用LU矩阵分解算法求解线性方程组A*x=B,具体做法是: 先对矩阵A进行行初等变换得到上三角矩阵U: 再求出下三角矩阵L,就 ...

  3. svd奇异值分解_传统推荐算法(一)SVD推荐(1)解读奇异值分解

    文章目录 写在前面 1. 从几何变换到奇异值分解 2. 代数角度理解奇异值与奇异向量 2.1 从正交基映射推导SVD 2.2 特征值分解求解奇异值和奇异向量 2.2.1 求解过程 2.2.2 推论 2 ...

  4. SVD奇异值分解(标题重复率过高)

    本篇为<深度学习>系列博客的第三篇,该系列博客主要记录深度学习相关知识的学习过程和自己的理解,方便以后查阅. 之前在看slam的时候就遇到过SVD但是没有做系统的学习,现在接触到深度学习, ...

  5. svd奇异值分解_SVD(奇异值分解)到底在干什么

    奇异值分解就是在低维空间中寻找最接近原矩阵 的低维矩阵 ,说白了就是数据降维. 奇异值分解是一种十分重要但又难以理解的矩阵处理技术,据人工智能的大牛吴恩达老师所说,在机器学习中是最重要的分解没有之一的 ...

  6. svd奇异值分解_Lecture 28 | 奇异值分解

    01 Singular Value Decomposition 奇异值分解 奇异值分解指任一mxn的矩阵A都可以分解为一个mxm酉矩阵U乘一个mxn对角阵Σ再乘一个nxn酉矩阵V共轭转置的形式. 下面 ...

  7. SVD的原理及python实现——正本清源

    SVD的原理及python实现 希望这篇文章能达到正本清源的目的.

  8. 三硬币问题建模及Gibbs采样求解(Python实现)

    大纲 三硬币问题建模及Gibbs采样求解(Python实现) 三硬币问题 二项分布与Beta分布 基础概念 二项分布 Beta分布 二项分布-Beta分布 三硬币问题建模过程 Gibbs 采样求解 p ...

  9. python 字典处理_python numpy求解积分python中的字典操作及字典函数

    字典 dict_fruit = {'apple':'苹果','banana':'香蕉','cherry':'樱桃','avocado':'牛油果','watermelon':'西瓜'} 字典的操作 W ...

  10. 大M法的python编程求解和python包求解

    大M法的python编程求解和python包求解 一.大M算法的求解步骤讲解 二.python编程求解 三.利用python包scipy的优化包optimize 四.用excel求解 五.分析结果 一 ...

最新文章

  1. Entity Framework简介
  2. pfx证书密码怎么查询_2019成考成绩查询时间通知!忘记账号和密码怎么办?
  3. 2.IDA-数据显示窗口(反汇编窗口、函数窗口、十六进制窗口)
  4. aws emr 大数据分析_DataOps —使用AWS Lambda和Amazon EMR的全自动,低成本数据管道
  5. hashtable和hashmap的区别?
  6. Python 内置模块之 asyncio(异步iO)
  7. (王道408考研操作系统)第三章内存管理-第二节3:页面置换算法2
  8. 【Linux入门学习之】grep命令详解
  9. MySQL绿色版mysql-5.7.17-winx64简洁安装教程
  10. 成人教育考试报名照片的尺寸是多少?大一寸照片怎么做?
  11. Maze环境以及DQN的实现
  12. access里面的表达式运用_Access 如何使用表达式生成器
  13. 2018-01-05-医药行业的IT革命探讨
  14. LibVLC —— 常用对象解析
  15. 运行阶段及面向对象技巧
  16. 服务器内存不足导致程序(tomcat)崩溃
  17. 活久见!这么好的图文电子书制作工具我竟然才发现
  18. 备忘3:爬取博主热门信息以及所有热门微博的点赞的用户信息
  19. C# 科学计数法的转换
  20. implode() 函数

热门文章

  1. python成语接龙_python——成语接龙小游戏
  2. 电脑打开应用程序提示配置系统未能初始化--解决方案
  3. C语言各个符号优先级(全)
  4. Verilog 实现四选一选择器
  5. 2021微信公开课PRO:微信视频号首次公布运营规则,“点赞”表情成视频号年度表情
  6. 医院病历html模板,三甲医院电子病历模板参考
  7. Adobe reader 的书签问题
  8. Python chardet模块
  9. 计算机驱动程序检测,检测到计算机制造商图形驱动程序对于显卡驱动程序
  10. matlab 滤波器设计 coe_巴特沃斯滤波器