矩阵特征值分解与奇异值分解含义解析及应用

原文在此，仅仅将原文的Matlab代码改为Python3版本。

特征值与特征向量的几何意义

矩阵的乘法是什么，别只告诉我只是“前一个矩阵的行乘以后一个矩阵的列”，还会一点的可能还会说“前一个矩阵的列数等于后一个矩阵的行数才能相乘”，然而，这里却会和你说——那都是表象。
矩阵乘法真正的含义是变换，我们学《线性代数》一开始就学行变换列变换，那才是线代的核心——别会了点猫腻就忘了本——对，矩阵乘法：

从上图可知，y方向进行了2倍的拉伸，x方向进行了3倍的拉伸，这就是B=[3 0; 0 2]的功劳,3和2就是伸缩比例。请注意，这时B除了对角线元素为各个维度的倍数外，非正对角线元素都为0，因为下面将要看到，对角线元素非0则将会发生切变及旋转的效果。
2. 切变

A = np.array([[0, 1, 1, 0, 0],[1, 1, 0, 0, 1]]) #原空间

plt.clf()
B1=np.array([[1,0],[1,1]])
B2=np.array([[1,0],[-1,1]])
B3=np.array([[1,1],[0,1]])
B4=np.array([[1,-1],[0,1]])

Y1 = np.dot(B1,A)
Y2 = np.dot(B2,A)
Y3 = np.dot(B3,A)
Y4 = np.dot(B4,A)

plt.subplot(221)
plt.plot(A[0],A[1],'-*',lw=2)
plt.plot(Y1[0],Y1[1],'-r*',lw=2)
plt.axis([-1,3,-1,3]);plt.grid(True)
plt.subplot(222)
plt.plot(A[0],A[1],'-*',lw=2)
plt.plot(Y2[0],Y2[1],'-r*',lw=2)
plt.axis([-1,3,-1,3]);plt.grid(True)
plt.subplot(223)
plt.plot(A[0],A[1],'-*',lw=2)
plt.plot(Y3[0],Y3[1],'-r*',lw=2)
plt.axis([-1,3,-1,3]);plt.grid(True)
plt.subplot(224)
plt.plot(A[0],A[1],'-*',lw=2)
plt.plot(Y4[0],Y4[1],'-r*',lw=2)
plt.axis([-1,3,-1,3]);plt.grid(True)
plt.show()

3. 旋转
所有的变换其实都可以通过上面的伸缩和切变变换的到，如果合理地对变换矩阵B取值，能得到图形旋转的效果，如下：

A = np.array([[0, 1, 1, 0, 0],[1, 1, 0, 0, 1]]) #原空间

plt.clf()
theta = np.pi/6
Bt=np.array([[np.cos(theta),np.sin(theta)],[-np.sin(theta),np.cos(theta)]])
Yt = np.dot(Bt,A)
plt.plot(A[0],A[1],'-*',lw=2)
plt.plot(Yt[0],Yt[1],'-r*',lw=2)
plt.axis([-1,3,-1,3]);plt.grid(True)
plt.show()

好，关于矩阵乘就这些了。那么，我们接着就进入主题了，对特定的向量，经过一种方阵变换，经过该变换后，向量的方向不变（或只是反向），而只是进行伸缩变化（伸缩值可以是负值，相当于向量的方向反向）？这个时候我们不妨将书上对特征向量的定义对照一遍：

数学教材定义： 设A是n阶方阵，如果存在 λ 和n维非零向量X，使
SVD分解实例
SVD之所以很有效，是因为：在很多情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上了。在这里，我们用图像简单的实践一下SVD的大妙处，下面是python对图像进行SVD分解的例子：

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
#设置中文字体
myfont = FontProperties(fname='C:/Windows/Fonts/msyh.ttc')

im = np.array(Image.open('lena.jpg').convert('L'),dtype=np.float32)
u,s,v = np.linalg.svd(im)

recv1 = np.dot(np.dot(u[:,:20],np.diag(s[:50])[:20,:50]),v[:50,:])  #svd取最高的20个特征值进行恢复
recv2 = np.dot(np.dot(u[:,:50],np.diag(s[:100])[:50,:100]),v[:100,:])  #svd取最高的50个特征值进行恢复
recv3 = np.dot(np.dot(u[:,:200],np.diag(s[:200])),v[:200,:])  #svd取最高的200个特征值进行恢复
plt.clf()
plt.figure(1)
plt.subplot(221)
plt.imshow(im,cmap=plt.get_cmap('gray'))
plt.title(u'原图',fontsize ='large',fontproperties=myfont)

plt.subplot(222)
plt.imshow(recv1,cmap=plt.get_cmap('gray'))
plt.title(u'恢复：左奇异20、右奇异50',fontsize ='large',fontproperties=myfont)
plt.subplot(223)
plt.imshow(recv2,cmap=plt.get_cmap('gray'))
plt.title(u'恢复：左奇异50、右奇异100',fontsize ='large',fontproperties=myfont)
plt.subplot(224)
plt.imshow(recv3,cmap=plt.get_cmap('gray'))
plt.title(u'恢复：左奇异200、右奇异200',fontsize ='large',fontproperties=myfont)
plt.show()

图注：SVD二维图像压缩恢复

如果按左下角的方式压缩原图，则存储量变为：50x50+100x100+50=12500，而存储原图像的存储量为512x512=262144，则压缩比为262144/12500=20.97，这里没有考虑存储数据类型的差异。

SVD分解相对于特征值分解的优势就是：

• 分解的矩阵可以是任意矩阵
• 在恢复信号的时候左右奇异值可以选择不同的维度

另外值得注意的一点：不论是奇异值分解还是特征值分解，分解出来的特征向量都是正交的。

奇异值分解与PCA

（节选自小节标题链接的原文，感谢原作者）

主成分分析在上一节里面也讲了一些，这里主要谈谈如何用SVD去解PCA的问题。PCA的问题其实是一个基的变换，使得变换后的数据有着最大的方差。方差的大小描述的是一个变量的信息量，我们在讲一个东西的稳定性的时候，往往说要减小方差，如果一个模型的方差很大，那就说明模型不稳定了。但是对于我们用于机器学习的数据（主要是训练数据），方差大才有意义，不然输入的数据都是同一个点，那方差就为0了，这样输入的多个数据就等同于一个数据了。以下面这张图为例子：

这个假设是一个摄像机采集一个物体运动得到的图片，上面的点表示物体运动的位置，假如我们想要用一条直线去拟合这些点，那我们会选择什么方向的线呢？当然是图上标有signal的那条线。如果我们把这些点单纯的投影到x轴或者y轴上，最后在x轴与y轴上得到的方差是相似的（因为这些点的趋势是在45度左右的方向，所以投影到x轴或者y轴上都是类似的），如果我们使用原来的xy坐标系去看这些点，容易看不出来这些点真正的方向是什么。但是如果我们进行坐标系的变化，横轴变成了signal的方向，纵轴变成了noise的方向，则就很容易发现什么方向的方差大，什么方向的方差小了。
一般来说，方差大的方向是信号的方向，方差小的方向是噪声的方向，我们在数据挖掘中或者数字信号处理中，往往要提高信号与噪声的比例，也就是信噪比。对上图来说，如果我们只保留signal方向的数据，也可以对原数据进行不错的近似了。
PCA的全部工作简单点说，就是对原始的空间中顺序地找一组相互正交的坐标轴，第一个轴是使得方差最大的，第二个轴是在与第一个轴正交的平面中使得方差最大的，第三个轴是在与第1、2个轴正交的平面中方差最大的，这样假设在N维空间中，我们可以找到N个这样的坐标轴，我们取前r个去近似这个空间，这样就从一个N维的空间压缩到r维的空间了，但是我们选择的r个坐标轴能够使得空间的压缩使得数据的损失最小。
还是假设我们矩阵每一行表示一个样本，每一列表示一个feature，用矩阵的语言来表示，将一个m * n的矩阵A的进行坐标轴的变化，P就是一个变换的矩阵从一个N维的空间变换到另一个N维的空间，在空间中就会进行一些类似于旋转、拉伸的变化。

而将一个m * n的矩阵A变换成一个m * r的矩阵，这样就会使得本来有n个feature的，变成了有r个feature了（r < n)，这r个其实就是对n个feature的一种提炼，我们就把这个称为feature的压缩。用数学语言表示就是：
但是这个怎么和SVD扯上关系呢？之前谈到，SVD得出的奇异向量也是从奇异值由大到小排列的，按PCA的观点来看，就是方差最大的坐标轴就是第一个奇异向量，方差次大的坐标轴就是第二个奇异向量…我们回忆一下之前得到的SVD式子：
在矩阵的两边同时乘上一个矩阵V，由于V是一个正交的矩阵，所以V转置乘以V得到单位阵I，所以可以化成后面的式子
将后面的式子与A * P那个m * n的矩阵变换为m * r的矩阵的式子对照看看，在这里，其实V就是P，也就是一个变化的向量。这里是将一个m * n 的矩阵压缩到一个m * r的矩阵，也就是对列进行压缩，如果我们想对行进行压缩（在PCA的观点下，对行进行压缩可以理解为，将一些相似的sample合并在一起，或者将一些没有太大价值的sample去掉）怎么办呢？同样我们写出一个通用的行压缩例子：
这样就从一个m行的矩阵压缩到一个r行的矩阵了，对SVD来说也是一样的，我们对SVD分解的式子两边乘以U的转置U'
这样我们就得到了对行进行压缩的式子。可以看出，其实PCA几乎可以说是对SVD的一个包装，如果我们实现了SVD，那也就实现了PCA了，而且更好的地方是，有了SVD，我们就可以得到两个方向的PCA，如果我们对A’A进行特征值的分解，只能得到一个方向的PCA。
PCA就是一种用于对数据进行降维的方法（降维肯定会丢失数据，只不过能在减少大量存储量的同时损失尽可能少），参见之前matlab对图像进行SVD分解的例子，更容易理解：实现了SVD就实现了PCA，PCA仅是SVD的包装。
PCA的应用很广，主要用在机器学习中对特征进行降维，还能用于去噪，下面两图是PCA降维和PCA去噪的例子（图片来自邹博PPT：北京9月秋季班·机器学习初步）
图注：PCA降维

降维说白了就是将信息通过投影到更低得多维度，这样必然会带来信息的丢失，但就如上图，这种信息的丢失却有时对分类没有影响，反而能降低识别算法的维度，提高速度，缓解所谓的维度灾难。

图注：PCA去噪

PCA去噪的前提是噪声的特征值会比信号的特征值小，即信噪比高的情况，否则PCA去噪会产生逆效果——把信号去掉了而噪声没去掉。

SVD其它

SVD还有其它很多方面的应用，通过查找资料，这里先做个简单的罗列，有机会再一个个研究：

• 求伪逆。我们知道，矩阵求逆要求矩阵必须是方阵，SVD可以用来求任意矩阵的逆运算，求任意矩阵的逆矩阵称为求伪逆
• 最小二乘法求解。凭着对《矩阵论》的零星的记忆，SVD算法就是因为能求伪逆所以用来求解最小二乘法。
• 基于SVD的文本分类。首先接触是从吴军老师的《数学之美》一书上看到的，大致是：通过TF/IDF(term frequency/inverse document frequency)构建“一百万篇文章和五十万词的关联性”的矩阵 A1000000x500000
  ，然后对A矩阵使用SVD分解之后，存储量减小，而且左奇异值矩阵和右奇异值矩阵以及奇异值矩阵的物理含义将非常明晰，此处博文 有简单介绍，更多参见吴军老师的《数学之美》

另外，开源视觉库OpenCV中也提供SVD分解的算法。

参考

1. 强大的矩阵奇异值分解(SVD)及其应用
2. http://blog.sciencenet.cn/blog-696950-699432.html ，关于SVD部分的讲解即出自此处。
3. http://www.cise.ufl.edu/~srajaman/Rajamanickam_S.pdf
4. July培训班PPT：“北京9月秋季班：贝叶斯、PCA、SVD”.
5. 戴华编著《矩阵论》 教材，科学出版社. 有时候明白一点之后再回头参考数学教材要比网上大把的抓好得多，豁然开朗。