Abstract: 本文介绍SVD，奇异值分解，应该可以算是本章最后的高潮部分了，也是在机器学习中我们最常用的一种变换，我们经常需要求矩阵的特征值特征向量，比如联合贝叶斯，PCA等常规操作，本文还有两个线性代数的应用，在图像压缩上，以及互联网搜索上。 Keywords: Singular Value Decomposition,JPEG2000,Eigenvalues,Eigenvectors

SVD分解

今天的废话关于学习知识，最近看到一种说法，我觉的非常的形象，有个大神（是谁我忘了），他说已知的知识像一个圆圈，而自己能感受的未知就是紧邻圆圈，圆外部的区域，当你知道的知识越来越多，圆圈不断扩大，圆周也随之扩大，所以你会越来越发现自己无知，那么就会更努力的去学习，所以越有知识的人越谦逊，尤其是对待知识上，尊重知识，探索未知领域是人类文明存在的根本动力。 ## 奇异值分解 (Singular Value Decomposition) SVD，熟悉的名字，如果不学习线性代数，直接机器学习，可能最先接触的就是SVD，所以我记得在写上个系列的博客的时候（CSDN，图像处理算法）就说到过SVD，当时还调侃了下百度，每次搜SVD出来的都是一把枪（报告政府，这个枪是穿越火线里面的，没超过1.7J）

这张分解图是我无意中发现的，ak47的发明人说过，如果一把枪，零件已经精简到最少了，那么这个才是精品，类似的意思上篇博客也说过，矩阵变换到最简单的形式，能够体现出其最重要的性质。 SVD，奇异值分解，与QR，LU，\(S\Lambda S^{-1}\) 等变换类似，其经过变换后会得到一个结构特异性质非凡的矩阵，SVD分解的结果和形式与对角化都非常相似，只是在形式和思路上更复杂，或者说如果说Jordan 是矩阵的对角化的扩展，因为有些矩阵特征向量不完全，那么SVD也是对角化的扩展，因为有些矩阵并不是方的。 所以SVD也是对角化，并且拥有比 \(A=S\Lambda S^{-1}\) 更完美的性质，但却是也复杂了一些，\(A=S\Lambda S^{-1}\) 有以下几个问题，需要完善： 1. S中特征向量一般不是正交的，除非A是对称矩阵 2. A并不是总有足够的特征值，这个是Jordan解决的问题，多个特征值相等，其对应于一个特征向量的时候，Jordan可以写成一块一块的对角矩阵 3. A必须是方的方的方的

Singular Vectors作为eigenvectors 的替代品，可以完美解决上述问题，但是作为代价，我们的计算过程会变得复杂，并且Singular Vectors有两组，\(u\) 和 \(v\)

\(u\) 对应的是\(AA^T\) 的特征向量，因为 \(AA^T\) 对称，所以 \(u\) 们可以选择相互正交的一组。 同理 \(v\) 对应 \(A^TA\) 的特征向量，因为\(A^TA\) 对称，所以 \(v\) 们也可以选择相互正交的一组。 这里注意是选择，因为你也可以选择不正交的，但是不正交的可能就会很麻烦了。

铺垫的差不多 ，然后我们有下面的这条重要性质，为什么会成立后面有证明，现在就告诉你SVD究竟是个啥子鬼： \[ Av_1=\sigma_1u_1\\ Av_2=\sigma_2u_2\\ \vdots\\ Av_n=\sigma_nu_n\\ \]

\(v_1,\dots,v_n\) 是\(A^TA\) 的特征向量，所以 \(v\) 是矩阵A的Row Space \(u_1,\dots,u_n\) 是\(AA^T\) 的特征向量，所以 \(u\) 是矩阵A的Column Space \(\sigma_1,\dots,\sigma_n\) 全部为正数，称为矩阵A的奇异值。

然后下面我们把 \(u\) 和 \(v\) 组合成矩阵 \(U\) 和 \(V\) ,那么根据对称矩阵的性质，\(U^TU=I\) 同理 \(V^TV=I\) 那么接下来我们来组合一下：

\[ AV=U\Sigma \\ A \begin{bmatrix} &&\\ v_1&\dots&v_r\\ && \end{bmatrix}= \begin{bmatrix} &&\\ u_1&\dots&u_r\\ && \end{bmatrix} \begin{bmatrix} \sigma_1&&\\ &\ddots&\\ &&\sigma_r \end{bmatrix} \]

矩阵形式就是这样喽，没什么解释的，就是上面计算的组合形式，但是注意这里有个很重要的参数，\(r\) 没错，就是矩阵的rank，这里rank表示了矩阵A的Singular Values的数量，所以上面计算从规模上是： \[ (m\times n)(n\times r)=(m\times r)(r\times r)\\ m\times r=m\times r \] 从矩阵相乘的规模上也能看出等式没有问题，但是这个r有的问题，可以肯定的是，有效的Singular vector有r组，但是这样与原始矩阵形状差的有点多，那么就补一补，虽然补的都是没用的，但是也算是整齐划一了，首先 \(\Sigma\) 中缺少的只能补0 ，所以对应的V就只能补A的Nullspace了，因为这样 \(AV\) 的补充部分是0,同理，为了配合V，U添加的是left nullspace，并且这些添加的无用值也要选择orthonormal的，以保证\(U^TU=I\) 和\(V^TV=I\)。

其实这里隐藏了一个重要的知识点，就是四个空间的那课，矩阵的rowspace和nullspace正交column space与left nullspace正交，而V本来是A的行空间正交基，那么添加的一定是Nullspace中的正交基，以保证矩阵正交，所以完美结合，（如果忘了四个空间点击查看）

所以更一般化的表示： \[ AV=U\Sigma \\ A \begin{bmatrix} &&\\ v_1&\dots&v_n\\ && \end{bmatrix}= \begin{bmatrix} &&\\ u_1&\dots&u_m\\ && \end{bmatrix} \begin{bmatrix} \sigma_1&&&\\ &\ddots&&\\ &&\sigma_r&\\ &&& \end{bmatrix} \] 规模上是，注意 \(\Sigma\) 不是方阵： \[ (m\times n)(n\times n)=(m\times m)(m\times n)\\ m\times n=m\times n \] \(\Sigma\) 被填充成立 \(m\times n\) 通过在矩阵中加入0来实现，新的矩阵U和V依旧满足 \(V^TV=I\)以及 \(U^TU=I\)

那么我们的A就可以分解了 \[ AV=U\Sigma\\ for:\;VV^T=I\\ so:\\ A=U\Sigma V^T\\ SVD\,\,\, is:\\ A=u_1\sigma_1 v_1^T+\dots+u_r\sigma_r v_r^T \] 其中\(u\) 是\(m\times 1\) 的 \(v^T\) 是 \(1\times n\) 的，所以A是 \(m\times n\) 的没有问题，并且所有 \(u_i\sigma_r v_i^T\) d的rank都是1，这就是Sigular Values Decomposition了，这里反复的验证规模的原因是因为A不是方阵，所以，在做乘法的时候要非常小心矩阵规模。那个小的只有r个有用值的SVD我们叫他reduced SVD（其实我觉得这个更有实际意义，毕竟这里面才有最重要的信息，新增的那些最后奇异值都是0了，也就没有啥作用了）可以表示为: \[ A=U_r\Sigma_r V_r^T \] 写了这么多，我们到现在还不知道Singular是怎么计算出来的，那么我们先给出结论，后面继续证明： \[ \sigma_i^2=\lambda_i \] 其中\(\lambda_i\) 是\(A^TA\) 和\(AA^T\) 的特征值。 那么要问\(A^TA\) 和\(AA^T\) 拥有相同的特征值，为什么？ 这个我真没想明白怎么证明，所以这个地方算个坑，会了再回来填

然后我们得到Singular Values后，我们把他们按照从大到小的顺序排列，然后写成上面SVD的形式： \[ \sigma_1 \geq \sigma_2 \geq \sigma_3 \dots \geq \sigma_n \]

下面举个小🌰 ： 什么时候SVD和对角化相等？ 当A是半正定或者正定矩阵的时候，\(S=U\) 并且 \(S^T=V^T\) 此时 \(\Lambda=\Sigma\) ,因为正定矩阵特征值为正，而且是对称的，所以 \(U=V=Q\)

下面介绍一个应用，大应用，为什么我会把这个应用写出来呢？因为他和图像有关，所以我们可以简单实践一下，还是从感性上认识一下SVD，然后再理论上完整的证明一下。

图像压缩 (Image Compression)

在介绍SVD之前，我们先来分析一个问题：图像的存储空间，一个图像假如是512x512的大小，灰度图像，每个像素占一个字节的话，这张图片那么会占据硬盘262144个字节，也就260多k个字节，如果按照23帧每秒，十秒钟大概要60兆，一分钟大概3.6G，在想想这个尺寸这么小，我们看的一般都是720P的，这样算的话硬盘根本不够用，一个东京热以后就再也不能一本道了，所以必须要压缩一下子，怎么压缩呢，这里介绍下JPEG的一个大致思路，就是矩阵分解： \[ A=u_1\sigma_1 v_1^T+\dots+u_r\sigma_r v_r^T\\ A=\sigma_1 u_1 v_1^T+\dots+\sigma_r u_r v_r^T\\ A=\sigma_1 S_1+\dots+\sigma_r S_r\\ where:\,S_i=u_i v_i^T \] 这样就是按照singular的大小，给所有singular vector排序，奇异值越大证明这个奇异值向量对原始数据影响越大，所以这两奇异向量组成的这个基就越重要，当然要提前加进去，所以我们如果选择一部分奇异值和奇异向量，比如\(N=100\) 也就是200个奇异向量和100个奇异值，那么一共是\(200\times 512+100=102500\) 比原来的260k减少了一半，如果用的更少那么减少的更多，当然图像质量也就有所损失了，写了个python程序，可以观察一下： 使用多组S来还原原始数据，使用S越多还原度越高，但需要的存储空间就越大，S就是上面\(u_iv_i\)的结果，可以看做图像的一个切片，视频中第一幅为原图，第二幅为若干张S相加得到的结果，最后一张是他们之间的差，我们可以暂时理解为压缩误差。

视频演示地址http://player.youku.com/embed/XMzE5NTIyNjQ4MA==

代码

import cv2
import numpy as np
from numpy import linalg
import copy
N=500
image=cv2.imread('lena.jpg',0)
cv2.imshow('src',image)
U,Sigma,V=linalg.svd(image)
sigma_size=len(Sigma)
image_s=[]
waittime=0
for i in range(N):

    if i>sigma_size:
        break
    image_s.append(U[:,i].reshape(len(U[:,i]),1)*V[i]*Sigma[i])
    if i>0:
        image_s[i]+=image_s[i-1]
    print 'total: ',i,' Singular Value'
    show_image=copy.deepcopy(np.uint8(image_s[i]))
    font = cv2.FONT_HERSHEY_TRIPLEX
    cv2.putText(show_image, `i`, (10, 500), font, 4, (255, 255, 0), 1, False)
    cv2.imshow('Compression',show_image)
    cv2.imshow('Different',np.uint8(image-image_s[i]))
    cv2.waitKey(waittime)
    waittime=100

上面公式中每个S都是一个rank=1的矩阵，如果两个这样的矩阵相加，rank就变成了2，n个相加就是rank=n，所以SVD分解后再组合有点像切片，每一片都是有规律的，同样的道理，SVD换成小波，那就有了JPEG2000，小波的基矩阵也是满足一些特殊性质的，后面可能会将，但是不确定，所以数据压缩可以理解为最关键的一步就是两个向量相乘，能够得到一个矩阵，这个矩阵组成了基础的片，然后多个片加权求和就得到了还原数据。 SVD的应用应该非常多这里就写了一个，比较重要直观的，下面还是要继续研究原理。

基和SVD(The Bases and the SVD)

书中并没有给出严格的证明和推到过程，老爷子还是走的启发式套路，我们来从基开始看，假设一个矩阵A是个2x2的矩阵，这样比较好计算，并且A是非奇异矩阵，也就是A可逆，rank=2，span成整个 \(\Re^2\) 空间，那么我们应该可以找到两个向量正交的单位向量 \(v_1,v_2\) 满足 \(u_1=\frac{Av_1}{|Av_1|},u_2=\frac{Av_2}{|Av_2|}\) 使得\(u_2\)和\(u_1\) 正交，并且是单位向量，那么就有： \[ A \begin{bmatrix}v_1&v_2\end{bmatrix}= \begin{bmatrix}Av_1&Av_2\end{bmatrix}= \begin{bmatrix}|Av_1|\cdot u_1&|Av_2|\cdot u_2\end{bmatrix}= \begin{bmatrix}u_1&u_2\end{bmatrix} \begin{bmatrix}|Av_1|&\\&|Av_2|\end{bmatrix}\\ AV=U\Sigma \] 总结下，这个构造基的过程是瞄准了目标去的，也就是说目标就是类似于构造 \(Ax=\alpha y\) 形状的一种形式，但是x和y要满足不同的关系，如果相等那么就是特征值，如果不相等就可以构造出多组相互正交形成奇异值，并且通过上面的构造过程，我们可以得知奇异值等于缩放u到单位向量的缩放比例 \(|Av|\) 其实可以用一种更直观的方法解答SVD的存在： 假设对任意矩阵A存在分解 \(A=U\Sigma V^T\) 并且其中\(U^TU=I\) \(V^TV=I\) 那么我们要证明U和V的存在即可： \[ A^TA=(U\Sigma V^T)^T(U\Sigma V^T)\\ A^TA=V\Sigma (U^TU) \Sigma V^T\\ A^TA=V\Sigma^2 V^T \] 虽然A不是对称的，但因为\(A^TA\) 是对称矩阵，存在正交矩阵Q使得 \(A^TA=Q\Lambda Q^T\)，那么这时候的V就是 \(A^TA\) 的Q这个是没问题的，至于 \(\lambda_i=\sigma_i^2 \geq 0\) 成立的原因是 \(A^TA\) 是个正定矩阵（A中各列线性无关），或者半正定矩阵（A中各列线性相关），所以其特征值 \(\lambda\) 必然非负数，所以根号后能得到奇异值，根据\(Av_i=\sigma u_i\) 可以求出剩下的 \(u_i\) ,当然这是理论上的方法，实际上的数值计算过程中可以避免 \(A^TA\) 这种大规模矩阵乘法。 回忆一下正定矩阵关于椭圆的那个例子 一个2x2正定矩阵对应二维空间一个椭圆（或者圆），其正交特征矩阵Q矩阵是对椭圆轴的旋转，特征值矩阵 \(\Lambda\) 是对轴的拉伸，那么我们的SVD有同样的功效，而且有过之无不及，思考： 作为\(A^TA\) 的正交特征矩阵\(V\)也是一个旋转矩阵，旋转的是圆的轴，\(V^T\) 当然就是反方向旋转，\(\Sigma\) 是对图形的拉伸，圆的拉成长的，接着 \(AA^T\) 的正交特征矩阵也是旋转，整个过程如下图：

所以一个2x2的可逆矩阵，对一个圆的操作就是先拉伸，然后旋转。图中来自Cliff Long and Tom Herm

上面讲解了如何求V，同样的道理也可以求U， \[ AA^T=(U\Sigma V^T)(U\Sigma V^T)^T\\ AA^T=(U\Sigma V^T)(V\Sigma U)\\ AA^T=U\Sigma^2 U^T \] U是\(AA^T\) 的正交特征矩阵，也就是说同一个 \(\Sigma^2\) 既是 \(AA^T\) 又是 \(A^TA\) 的特征值，所以上面那个疑问也得到了证明。 我们重新梳理一下这个证明过程，我们首先假设结论成立，来找到使结论成立的条件，也就是V和U矩阵，结果很理想的，我们找到了，所以原来结论成立，成立的条件就是我们刚找到的这两个矩阵（如果算上奇异值，可以说是三个矩阵），思考过程可以通过2x2构造那里来推出简单情况下，来验证我们的结论，然后再推广到任意矩阵。 至此SVD的基本来路我们已经算是摸着门了，所以可以深入开发下V和U矩阵，这两个矩阵是方阵，那么组成这些矩阵的向量门也是很有来历的，总结个表格，前面在第一部分有说过，就是reduced SVD那部分，这里在啰嗦一边：

前r列对应的是 \(A^TA\) 和 \(AA^T\) 的特征向量，因为其间的正交关系，所以SVD是没有什么问题的。 严格的证明： \[ A^TAv_i=\sigma_i^2v_i \] 这个是可以作为条件使用的，其成立的必然原因是 \(A^TA\) 是正定或半正定矩阵，所以必然存在n个大于等于0的实数特征值，这样可以得到\(v_i,\sigma_i\) ，在这个条件下我们目标是证明：存在\(u_i\) 使得\(Av_i=\sigma_i u_i\) 成立

对条件两边同时乘上 \(v_i^T\) 后： \[ v_i^TA^TAv_i=\sigma_i^2v_i^Tv_i\\ ||Av_i||^2=\sigma_i^2 so:\\ ||Av_i||=\sigma_i\\ \] 对条件两边同时乘A得到： \[ AA^TAv_i=\sigma_i^2Av_i\\ set\,unit\,vector:\;u_i=\frac{Av_i}{\sigma}\\ AA^Tu_i=\sigma_i^2 u_i\\ \] 上面两个过程及其精妙，尤其是下面这个，用\(u_i=\frac{Av_i}{\sigma}\) 进行置换后，得到\(u_i\) 是\(AA^T\) 的特征向量，然后得出结论，存在 \(u_i=\frac{Av_i}{\sigma}\) 使得命题成立，而且u就是\(AA^T\)的特征向量，并且相互正交： 下面证明u相互正交: \[ u_i^Tu_j=(Av_i)^T(Av_j)=v_i^T(A^TAv_j)=v_i^T(\sigma^2v_j)=0 \] QED

然后老爷子在书上发话了，这是这本书最高潮的部分了，也是基础理论的最后一步了，因为他用到了所有的前面的理论： 1. 四个子空间的维度（我们上面那个表） 2. 正交 3. 正交基来对角化矩阵A 4. 最后得到SVD \(A=U\Sigma V^T\)

于是我写了一天这篇博客，比之前看书收获了更多，也可能有写纰漏，难免的因为水平有限，而且这个是比较精髓的部分，自然难度也比较大

搜索网络(Searching the Web)

这里讲了个应用，但是我实在不想写了，已经精疲力尽了，所以我打算后面选几个应用写，这个作为候选，这里略 ## Conclusion 线性代数的高潮算是来了，但是后面还有一个比较有意思的主题，也是很常用的，叫做线性变换，也可以作为切入线性代数的一个切入点，我们这个系列的博客是从线性方程组开始的，当然也可以通过线性变换引出矩阵，我们下一篇来讲解这个，待续。。