人工智能里的数学修炼 | 矩阵的花样分解：特征值分解（EVD）、相似对角化、QR分解、Schur分解、奇异值分解（SVD）的概念纠缠与详解

前言

在高等代数里，矩阵分解是一个十分基础与重要的内容，任何一个学校对于理工科的研究生教育都会开设相应的课程，如：矩阵分析、矩阵论、线性系统等。看了不少社区的问答、笔记和博客，在它们的基础上加入一些自己的理解，写下这篇概念详解,博客中借鉴了不少前人的观点，这里感谢他们的付出

一、特征值分解（EVD）

把特征值分解（EigenValue Decomposition，EVD）放在最前面说是因为这是一个相对来说更基础的概念，后面的一些分解方式会在这一基础上进行推导

1.1 特征值分解、特征值、特征向量

特征值分解，是指将矩阵分解为由其特征值和特征向量表示的矩阵之积的方法。一个NxN的方阵A可以有下式：

Av=λv

此时λ被称为方阵A的特征值，v被称为特征值 λ 对应的特征向量。上式可以化为

(λI−A)v=0

令 det|λI−A|=0 (此式被称为矩阵A的特征多项式)我们可以得到下式：

(λ−λ1)n1(λ−λ2)n2⋅⋅⋅(λ−λl)nl=0
∑li=1ni=N

这样我们便可以得到矩阵A的 l 种不同的特征值，每种特征值有 ni 个重复,共N个特征值。对于特征向量的求解，我们放在下一小节。从这里我们可以看出矩阵A可以有多个特征值和多个特征向量，每一个特征值对应一个或多个特征向量(此时该对应的特征值有重复)，不同特征值对应的特征向量之间线性不相关，重根对应的多个特征向量之间不一定不相关。特征值分解将矩阵分解成为了如下形式

A=QΛQ−1

对于没有重根的情况 Λ=diag(λ1,λ2,⋅⋅⋅,λn) ，对于有重根的情况 Λ 为Jordan标准型。这里的Q是以特征向量为列向量的NxN矩阵

1.2 特征向量的求解

求解特征向量有多种方法这里介绍最简单也是最常用的方法

• 计算A的特征多项式det|λI−A|=0，从而求得特征值 λi

• 对于单根特征值来说，求齐次方程 (λiI−A)vi=0，即可求得该特征值对应的特征向量

• 对于有重根的特征值来说，可以使用一下公式，依次迭代求解
  

  v1:Av1=λv1
  v2:Av2=v1+λv2
  v3:Av3=v2+λv3
  ⋅⋅⋅
  vN:AvN=vN−1+λvN

1.3 特征值与特征向量的意义解释

上面介绍了特征值分解、特征值与特征向量的数学定义与求解方式，但是看到这里，可能读者对于特征值分解的具体意义与作用还是模糊的，这也确实是一个比较难理解的地方。
我们知道，矩阵乘法其实是对应着一个线性变换，是把任意一个向量变成另一个方向或者长度的新向量。在这个变换中，原向量主要发生旋转、伸缩的变化。如果矩阵对某一个向量或某些向量只发生伸缩变换，而不对这些向量产生旋转效果，那么这些向量就称为这个矩阵的特征向量，伸缩的比例就是特征值
例如，对于矩阵 M=[3001] 它对应的线性变换为

M矩阵是对称的，所以这个变换是一个对x,y轴的拉伸变换，此处将原方块在x轴方向拉长了三倍。对于不是对称的情况，如 M=[1011] ，对应的线性变换为

看到这里，大家应该明白了，特征向量与特征值对于一个矩阵的意义，每一个特征向量都对应着这个矩阵在对另一个对象作线性变换时所产生变换的方向，而特征值则表示着这个变化的大小。也就是说，矩阵A的信息可以由其特征值和特征向量表示。对于矩阵为高维的情况下，这个变换就对应着很多方向的变化，我们通过特征值分解得到的对应特征值较大的前N个特征向量，便表示了这个矩阵最主要的N个变化方向。我们利用这前N个变化方向，就可以近似这个矩阵的变换，而著名的主成分分析（Principle Conponent Analysis,PCA）便是基于这一原理实现的。
总而言之，通过特征值分解我们可以得到特征值与特征向量，特征向量表示这个特征是什么，而特征值表示着这个特征有多重要。同时，我们要意识到特征值分解也有一定的局限性，它的对象必须是方阵，而实际中的矩阵往往不是方阵，后面要说的奇异值分解将其演化到了普通形式的矩阵

二、相似对角化

2.1 相似矩阵的定义

设A、B是两个n阶方阵，如果存在可逆矩阵T，使得 T−1AT=B 则称A与B相似，记作A～B，从A到B的这种变化称为相似变换，T称为相似变换矩阵。
矩阵的相似关系是一种等价关系（并不是相等），相似矩阵满足以下特性

• 自反性：A～A
• 对称性：若A～B，则B～A
• 传递性：若A～B，B～A，则A～C

2.2 相似对角化的条件与推论

N阶方阵A与对角阵相似的前提是：A有N个线性无关的特征向量。可以对角化即意味着存在某组基，使得这个矩阵所代表的线性变换在这组基的每一个方向上都是伸缩变换（复向量上的复“伸缩变换“近似于在某种意义上非刚性但依然线性的伸缩旋转），不能对角化即意味着找不到这样一组基

注：对于相似变换来说，有一些东西是变换前后没有改变的

• 若 A∼Λ=diag(λ1,λ2,⋅⋅⋅,λN) ,则A与 Λ 的特征值相同，Λ 的主对角线元素 λ1,λ2,⋅⋅⋅,λN 为A的全部特征值。相似变换的变换矩阵为 P=(p1,p2,⋅⋅⋅,pN)，则列向量p1,p2,⋅⋅⋅,pN依次是λ1,λ2,⋅⋅⋅,λN对应的特征向量

• 相似变换矩阵不唯一，因为特征向量的排列顺序可以发生变化

• A∼Λ，若不计Λi的排列顺序，则Λ唯一，称为A的相似标准型

基于相似转换的定义以及以上特性，我们可以得到一些重要的推论

2.2.1 推论一

若N阶方阵A的n个特征值互不相同，则A与对角阵相似，这是显然的，因为对应于不同特征值的特征向量线性无关，由此N个特征向量可以产生N个线性无关的基向量

2.2.2 推论二

设λ1,λ2,⋅⋅⋅,λN是A的l个互不相同的特征值，重数依次为r1,r2,⋅⋅⋅,rN，且有r1+r2+⋅⋅⋅+rl=N,则A可以相似对角化的充分必要条件为：A的ri重特征值λi恰有ri个线性无关的特征向量(i=1,2,⋅⋅⋅,l)

注：此处的ri又被称为代数重数，而实际的线性无关的特征向量的个数称为几何重数Ri。我们有Ri<=ri总是成立的，推论二是指只有当Ri=ri时方阵A才可以相似对角化

推论二的证明较为繁琐，感兴趣的可以点击这里

2.2.3 推论三

如果N阶方阵A可以对角化，则A的秩等于A的非零特征值的个数。这也是很容易理解的，若A可以对角化，设与其相似的对角阵为Λ，即存在可逆矩阵P，使得P−1AP=Λ。因此A与Λ等价，由此 rank(A) = rank(Λ) ，所以Λ对角线上的非零元素个数为rank（A）。又因为A与Λ相似，所以A的特征值与Λ的特征值相同，所以A的秩矩阵的秩等于A的非零特征值的个数

那么，对于那些不能对角化的矩阵我们该如何理解呢？这里我借用知乎上一位匿名答主的回答向大家解释:

可对角化的矩阵举例如下，左图为原图，右图为经过可以对角化的矩阵线性变换后的结果，箭头表示着伸缩的方向，长度表示变换的大小

不能对角化的两个变换如下，注意这种时候发生了切变，下图的变换均不能表示为各个方向独立的伸缩，也不能表示为带伸缩的旋转。图中不从原点出发的箭头表示切变的大小和方向

同时我们也应该注意到以上的四幅图中，第二幅图可以对角化的矩阵，它的jordan标准型是对角化的，而三、四幅图，它的Jordan标准型不是对角化的。实际上第二幅图的Jordan标准型就是变换矩阵对角相似化后的对角矩阵（在这里我们也称它为对角标准型），对角标准型是Jordan标准型的特例
相似对角化是矩阵分析当中最基础也是最重要的内容，在高等代数与工程问题中被广泛运用，可以极大的简化矩阵的运算，如计算方阵的高次幂、求矩阵行列式等

2.3 实对称矩阵与相似对角化

看了第二节之后，我们知道，对于一般的方阵我们常常无法进行相似对角化来简化矩阵，同时对于高维矩阵来说，对角化的条件难于判断。在这一小节，要介绍一类一定可以实现对角化的矩阵——实对称矩阵

2.3.1 实对称矩阵的特征值与特征向量

实对称矩阵的特征值为实数，对应的特征向量为实向量。设λ是A的特征值，v是对应λ的特征向量，若A为实对称矩阵，则有一下性质

AT=A , A¯=A
Av=λv,v≠0
可以据此推导
v¯TAv=v¯TATv=(Av¯)Tv=(Av⎯⎯⎯⎯)Tv
λv¯Tv=λ¯v¯Tv
所以
λ=λ¯
λ为实数，因此det|A−λE|x=0必有实的基础解系，从而对应的特征向量可以取实向量

2.3.2 实对称矩阵正交相似于对角矩阵

若A是N阶实对称矩阵，λ1,λ2,⋅⋅⋅,λN是A的特征值，则有正交矩阵Q，使得

Q−1AQ=QTAQ=Λ=diag(λ1,λ2,⋅⋅⋅,λN)
此时称A正交相似于Λ。这里有一个非常重要的结论：实对称矩阵的所有特征向量之间是线性无关的，我们之前提到的方阵A的特征向量之间线性无关，仅仅限于不同特征值之间的特征向量，对于同一特征值的不同特征向量之间，普通的方阵并不能保证他们线性无关

那么，这个正交变换矩阵该如何求出来呢？可以按照以下步骤：

• 求出A的全部特征值
• 求出每个特征值所对应的全部特征向量
• 以特征向量为列向量写出变换矩阵
• 使用Schmidt 正交化将变换矩阵正交化，单位化，得到正交矩阵Q

2.4 相似对角化与特征值分解的区别

相似对角化与特征值分解乍看上去是极为相似的，它们都需要用到特征值与特征向量的概念，但其实有较大差别

• 目的：特征值分解的目的在于矩阵分解，求得矩阵对应的特征值与特征向量；而相似对角化的目的在于通过变换矩阵的线性变换将原方阵转换为对角矩阵

• 条件：所有的方阵都可以进行特征值分解得到对应的特征值与特征向量；只有当方阵的几何重数与代数重数相等（方阵的最小多项式无重根）时，方阵才可以实现对角化

• 结果：通过特征值分解得到的特征向量与特征值可以构成对角标准型与jordan标准型（前者是后者的特例），其中Jordan标准型不是对角矩阵；而相似对角化得到的矩阵一定是对角矩阵

三、QR分解

QR分解是目前求取一般矩阵全部特征值的最有效并且广泛应用的办法,它是将矩阵分解成为一个正交矩阵Q和一个上三角矩阵R，所以称为QR分解。这一分解方法除了用于特征值计算外，在参数估计和通信领域也有着广泛的应用

3.1 QR分解的定义与推导

若AϵRn×m，且 n>=m，则存在列正交矩阵 QϵRn×m和上三角矩阵 RϵRm×m，使得 A=QR。
当 m=n 时 ，Q为正交矩阵。如果A是非奇异的n×n 矩阵，则R的所有对角线元素均为正，并且在这种情况下，Q和R是唯一的。若A是复矩阵，则Q和R取复值

关于QR分解的证明，这里根据Schmidt 正交化的方法给出当A为列满秩情况下的证明：

• 将 A 表示为 A=[x1,x2,⋅⋅⋅,xm]
• 由于A满秩，所以 xi 之间线性独立，通过Schmidt 正交化我们可以得到一组正交向量和一个上三角矩阵如下

[u1u2⋅⋅⋅um]=[x1x2⋅⋅⋅xm]⎡⎣⎢⎢⎢t1100⋅⋅⋅⋱0t1m⋮tmm⎤⎦⎥⎥⎥
U=AT

• 这里的T矩阵是Schmidt 正交化的变换矩阵，由于
  tii=∥∥xi−∑i−1j=1<uj,xi>uj∥∥−1
  矩阵T是非奇异的，同时T−1也同样为上三角矩阵，令Q=U，R=T−1,我们便可以得到A=QR

对于矩阵的QR分解其实还有很多种方法除了上面提到的Schmidt 正交化，还有利用矩阵的初等变换、利用Givens变换求QR分解等方法

3.2 QR分解的应用

QR分解在实际工程问题中得到了广泛的应用，其核心还是围绕着利用QR分解求解矩阵的特征值进行的，这里列举出一些常见的例子

• 利用QR分解求取矩阵的全部特征值

• 使用QR分解解决最小二乘法

• QR分解在参数估计中的应用

• QR分解在通信系统中的应用

• 简化PCA算法用于人脸识别

四、Schur分解

基于QR分解我们可以推导出Schur分解，同时，Schur分解也是很多重要理论推导的基石，是很多重要定理证明的出发点。在介绍Schur分解之前，我们需要先了解一下什么是酉矩阵(unitary matrix)

4.1 什么是酉矩阵？

4.1.1 “等距”(isometry)

对于一个矩阵 UϵFn×m ,如果 UHU=I (H为共轭转置)我们便称 U 为一个等距(isometry)，它的所有列向量是正交的。等距作为一个线性变换时是保内积运算，同时也是保模长运算，即

<Ux,Uy>=<x,y>
∥∥Ux∥∥=<Ux,Uy>1/2=<x,y>1/2=∥∥x∥∥

4.1.2 “协等距”(co-isometry)

对于一个矩阵 UϵFn×m ,如果 UUH=I 我们便称 U 为一个协等距(isometry),它的所有的行向量是正交的。

4.1.3 酉矩阵(unitary matrix)

一个矩阵如果既满足等距，又满足协等距，我们便就称它为酉矩阵，它的最大特点在于U−1=UH。酉矩阵其实是正交矩阵在复数域上的推广，它满足

UUH=UHU=I

4.2 Schur分解的定义与推导

方阵 AϵFn×n 具有特征值 λ1,λ2,⋅⋅⋅,λn ，则存在一个酉矩阵 UϵCn×n 使得

UHAU=T
T为一个上三角矩阵，它的对角线元素 tii=λi。现在来证明Schur分解的存在

• 记 xi 为对应于特征值 λi 的特征向量，令 X1=[x1,x2,⋅⋅⋅,xn]

• 对 X1 进行QR分解，可以得到X1=Q1R1，Q1这里是酉矩阵，R1是上三角矩阵。要注意的是 Q1 的第一列仍然是 A 对应于特征值 λi 的特征向量，因此有
  

  QH1AQ1=[λ10∗A1]
  这里 A1ϵC(n−1)×(n−1)，它的特征值为 λ2,⋅⋅⋅,λn

• 使用同样的步骤，我们又可以得到一个酉矩阵 Q2ϵC(n−1)×(n−1)，得到
  

  QH2A1Q2=[λ20∗A2]
  再令
  U2=[100Q2]
  于是有
  U_{2}^{H}Q_{1}^{H}AQ_{1}U_{2} =
  ⎡⎣⎢⎢λ100∗λ20∗∗A2⎤⎦⎥⎥

  -重复上述步骤，得到酉矩阵 QiϵC(n−i+1)×(n−i+1) 可以使
  

  QHiAi−1Qi=[λi0∗Ai]
  以及 UiϵCn×n
  Ui=[I00Qi]

• 最后矩阵U=Q1U2⋅⋅⋅Un−1即为所求的酉矩阵

4.3 Schur分解的缺陷

Schur分解将原方阵转化为了一个对角线为特征值的上三角矩阵，在这一章节的开头已经说过Schur分解是很多重要定理推导的基石与出发点。但是矩阵的Schur分解，在更多意义上是一种理论上的存在，在实际中通常不方便通过有限次运算得到，真正要计算时，一般也是通过迭代的方法进行逼近

五、奇异值分解（SVD）

我们在第一章节中就介绍了特征值分解，这是一个很好的提取矩阵特征的方法，但是，特征值分解只能对方阵使用，在实际问题中，我们解决问题所遇到的矩阵往往都不是方阵，那么，我们怎么样来描述一个普通矩阵的重要特征呢？奇异值分解为我们提供了一种可以用于任意矩阵分解的方法，这是这篇文章中最重要也是最核心的部分

5.1奇异值分解的定义与推导

对于一个矩阵 AϵFm×n，可将其写成如下形式

A=UΛVT
其中 UϵFm×m 的酉矩阵，被称为左奇异向量；ΛϵFm×n的半正定对角矩阵； VHϵFn×n 是V的共轭转置，被称为右奇异向量。这样的分解就叫做奇异值分解，Λ对角线上的元素λi即为原矩阵A的奇异值,奇异值一般按照从大到小排列，即
λ1>=λ2>=⋅⋅⋅>=λmin(n,m)

奇异值分解的推导可以从特征值分解开始

• 首先，我们对n阶对称方阵 ATA 作特征值分解，得到
  

  ATA=VΛVT

• 通过特征值分解我们得到一组正交基 V=(v1,v2,⋅⋅⋅,vn)，满足如下性质
  

  (ATA)vi=λivi
  由于ATA为对称矩阵，vi之间两两相互正交，所以有
  <Avi,Avj>=vTi(ATA)vj=vTiλjvj=λjvTivj=0

• 因为rank(ATA)=rank(A)=r，我们可以得到另一组正交基 Av1,Av1,⋅⋅⋅,Avr将其标准化有
  

  ui=Avi|Avi|=1λ√Avi
  Avi=λi‾‾√ui=δiui
  注：
  |Avi|2=<Avi,Avi>=λivTivi=λi

• 将向量组 (u1,u2,⋅⋅⋅,ur)扩充为Fm中的标准正交基(u1,u2,⋅⋅⋅,ur,⋅⋅⋅,um)则：
  

  AV=A(v1,v2,⋅⋅⋅,vn)=(Av1,Av2,⋅⋅⋅,Avr,0,⋅⋅⋅,0)=(δ1u1,δ2u2,⋅⋅⋅,δrur,0,⋅⋅⋅,0)=UΛ
  由此，可以得到奇异值分解的形式
  A=UΛVT

5.2 奇异值分解的求解

我们现在已经知道了奇异值分解的具体形式，那么奇异值和奇异向量到底怎样求解呢？

5.2.1奇异值的计算

• 对于较小维度的矩阵，我们可以从奇异值分解的推导中看出，奇异值δi=λi‾‾√。于是可以通过求解原矩阵的转置与其自身相乘得到的矩阵的特征值，再对该特征值求平方根的方法求得矩阵的奇异值

• 高纬度的矩阵的奇异值的计算是一个难题，是一个O(N^3)的算法，随着规模的增长，计算的复杂度会呈现出3次方的扩大，感兴趣的朋友可以看这里

5.2.1奇异向量的计算

在奇异值分解中，有一个十分重要的推论，那就是在式A=UΛVT里，U的列向量为AAT的特征向量，V的列向量为ATA的特征向量。知道这一推论，我们在计算出特征值之后就可以较为方便的求出矩阵的特征向量

5.3 奇异值分解的意义

奇异值分解的目的在于，找到一组正交基，使得矩阵在变换过后是正交的，这是奇异值分解的精髓所在。

5.3.1 数据降维压缩

奇异值往往对应着矩阵中隐含的重要信息，且重要性和奇异值的大小正相关，对于这一解释的最大应用就在与图像的压缩。可以将矩阵表示为若干个秩一矩阵之和

A=δ1u1vT1+δ2u2vT2+⋅⋅⋅+δrurvTr
我们知道，矩阵的奇异值一般按照降序排列即
λ1>=λ2>=⋅⋅⋅>=λmin(n,m)>0
一般来说，前10%甚至1%的奇异值之和就可以占到全部奇异值的99%以上，也就是说，我们可以使用奇异值较大的一些特征来表示图像，省略去较小的奇异值(绝大多数奇异值)，来实现图像的降维压缩，这里以知乎上的一名匿名网友的回答为例

左上：原图 ； 右上：保留前五项 ； 左下：保留前二十项 ； 右下：保留前五十项

原图的维度远远超过10000维，而通过奇异值分解，从上图可以看出，我们只需要保留前50项，就可以很好的复原图像，即实现图像的压缩。除了实现对图像的压缩外，奇异值分解在好友推荐算法，工业过程故障诊断等领域均有广泛应用。

5.3.2 几何的线性变换

奇异值分解的几何意义与特征值分解也极为相似，即奇异向量代表着线性变换的方向，而奇异值表示着在这个方向上变化的大小。这里举一个较为有名的椭圆变换为例

假设矩阵A的奇异值分解为

A=[u1u2][3001][vT1vT2]
其中 u1,u2,v1,v2 是二维平面的向量，根据奇异值分解的性质，u1,u2 线性无关，v1,v2线性无关。那么对二维平面上任意的向量x，都可以表示为：x=a1v1+a2v2，当A作用在x上时
y=Ax=A[v1v2][aT1aT2]=[u1u2][3001][vT1vT2][v1v2][aT1aT2]=3a1u1+a2u2

令η1=3a1, η2=a2，我们可以得出结论：如果x是在单位圆ai21+ai22=1上，那么y正好在椭圆η21/32+η22/12=1上。这表明：矩阵A将二维平面中单位圆变换成椭圆，而两个奇异值正好是椭圆的两个半轴长，长轴所在的直线是span{u1}，短轴所在的直线是span{u2}

推广到一般情形：一般矩阵A将单位球∥x∥2=1变换为超椭球面Em={y∈Fm: y=Ax, x∈Fn, ∥x∥2=1}，那么矩阵A的每个奇异值恰好就是超椭球的每条半轴长度。

---

参考文献：

• 特征值分解部分
  [1] http://www.cnblogs.com/LeftNotEasy/archive/2011/01/19/svd-and-applications.html
  [2] http://blog.csdn.net/jinshengtao/article/details/18448355
  [3] https://wenku.baidu.com/view/3ec0a4ddaeaad1f346933f42.html
  [4] http://www.doc88.com/p-9009713157157.html
  [5] https://wenku.baidu.com/view/f14c18215901020207409c97.html
  [6] https://www.zhihu.com/question/22548386
  [7] https://github.com/LiangjunFeng/Machine-Learning/blob/master/8.PCA.py
  [8] Bhushan Datta K. Linear system theory and design, by Chi‐Tsong Chen, Oxford University Press, New York, 1999, 334 pages, ISBN 0‐19‐511777‐8[J]. International Journal of Robust & Nonlinear Control, 2015, 10(15):1360-1362.
• 相似对角化部分
  [1] https://wenku.baidu.com/view/c22c4a708e9951e79b892760.html
  [2] https://www.zhihu.com/question/36187051
  [3] https://wenku.baidu.com/view/347b97466edb6f1aff001f98.html
  [4] https://wenku.baidu.com/view/21cd4a9f32d4b14e852458fb770bf78a65293ac2.html
  [5] https://en.wikipedia.org/wiki/Shear_stress
  [6] http://www.doc88.com/p-7178359679199.html
  [7] https://wenku.baidu.com/view/f83d600084254b35effd3401.html
  [8] https://wenku.baidu.com/view/41a43f0316fc700abb68fca5.html
  [9] https://en.wikipedia.org/wiki/Orthogonal_matrix
• QR分解部分
  [1] https://wenku.baidu.com/view/bf00c82cf8c75fbfc77db2da.html
  [2] http://blog.sina.com.cn/s/blog_3f41287a0101ke2s.html
  [3] http://blog.csdn.net/zhaogang1993/article/details/42562009
  [4] http://johnhany.net/2016/05/from-qr-decomposition-to-pca-to-face-recognition/
• Schur分解部分
  [1] https://baike.baidu.com/item/%E9%85%89%E7%9F%A9%E9%98%B5/2967660
  [2] https://wenku.baidu.com/view/65aff9174431b90d6c85c7b9.html
  [3] https://wenku.baidu.com/view/257d4dc10722192e4436f654.html
  [4] http://www.doc88.com/p-6791891160524.html
  [5] http://blog.csdn.net/makeway123/article/details/17803991
  [6] https://www.zhihu.com/question/20903131
• 奇异值分解部分
  [1] http://blog.csdn.net/zhongkejingwang/article/details/43053513
  [2] https://www.zhihu.com/question/22237507
  [3] http://www.cnblogs.com/LeftNotEasy/archive/2011/01/19/svd-and-applications.html
  [4] http://blog.csdn.net/jinshengtao/article/details/18448355
  [5] https://wenku.baidu.com/view/38693ef2e109581b6bd97f19227916888486b916.html
  [6] https://wenku.baidu.com/view/3ec0a4ddaeaad1f346933f42.html
  [7] http://www.cnblogs.com/liangzh/archive/2013/03/05/2841025.html
  [8] http://www.ams.org/samplings/feature-column/fcarc-svd
  [9] http://charleshm.github.io/2016/03/Singularly-Valuable-Decomposition/
  [10] http://blog.csdn.net/zhuiqiuk/article/details/69390357

---