MIT 18.065—机器学习中的矩阵方法12 计算特征值和奇异值

数据分析、信号处理和机器学习中的矩阵方法

第12讲 计算特征值和奇异值

这并不是数值线性代数课程，但我们需要探讨如何计算特征值和奇异值。你可以调用eig或svd或Python以及Julia中的等效程序，但仍需要了解相关知识。

特征值计算

首先谈QR分解，矩阵 \[{{\boldsymbol A}_{0}}={{\boldsymbol Q}_{0}}{{\boldsymbol R}_{0}}\] 。矩阵R在非正交矩阵A和正交矩阵Q之间建立了联系。

取 \[{{\boldsymbol A}_{1}}={{\boldsymbol R}_{0}}{{\boldsymbol Q}_{0}}\] ，则 \[{{\boldsymbol A}_{1}}={{\boldsymbol R}_{0}}{{\boldsymbol Q}_{0}}={{\boldsymbol R}_{0}}{{\boldsymbol A}_{0}}\boldsymbol R_{0}^{-1}\] ，两矩阵相似，因此矩阵A1和矩阵A0具有相同的特征值。继续在这个操作构造矩阵，矩阵An对角线以下的元素会变得越来越小，而对角线元素会越来越接近于特征值。

例： \[{{\boldsymbol A}_{0}}=\left[ \begin{matrix} \cos \theta & \sin \theta \\ \sin \theta & 0 \\ \end{matrix} \right]\] ，则有 \[{{\boldsymbol A}_{1}}=\left[ \begin{matrix} \cos \theta (1\text{+}{{\sin }^{2}}\theta ) & -{{\sin }^{3}}\theta \\ -{{\sin }^{3}}\theta & -{{\sin }^{2}}\theta \cos \theta \\ \end{matrix} \right]\] ，对角线下元素取了立方，再变换一步会再取立方变为9次方……

这个方法迅速淘汰了所有其他用于计算特征值的方法。但做数值计算的人希望进一步改进。

引入一个平移矩阵，得到 \[{{\boldsymbol A}_{0}}-s\boldsymbol I\] ，它具有和原矩阵相同的特征向量，而特征值改变了s， \[({{\boldsymbol A}_{0}}-s\boldsymbol I)\bold v={{\boldsymbol A}_{0}}\bold v-s\bold v=(\lambda -s)\bold v\] 。

对该矩阵进行QR分解，得到 \[{{\boldsymbol A}_{0}}-s\boldsymbol I={{\boldsymbol Q}_{0}}{{\boldsymbol R}_{0}}\] ，将分解矩阵反向相乘并反向平移变换构造矩阵 \[{{\boldsymbol A}_{1}}={{\boldsymbol R}_{0}}{{\boldsymbol Q}_{0}}+s\boldsymbol I\] 。则矩阵A1和矩阵A0为相似矩阵，具有相同特征值。

\[{{\boldsymbol A}_{1}}={{\boldsymbol R}_{0}}{{\boldsymbol Q}_{0}}+s\boldsymbol I={{\boldsymbol R}_{0}}({{\boldsymbol A}_{0}}-s\boldsymbol I)\boldsymbol R_{0}^{-1}+s\boldsymbol I={{\boldsymbol R}_{0}}{{\boldsymbol A}_{0}}\boldsymbol R_{0}^{-1}\]

这样处理的好处是使得特征值收敛得更快。

假设我们的矩阵已经有一些零，比如说通过一些操作得到了有一条下对角线的矩阵。

实际上是很难把矩阵通过相似变换直接得到上三角矩阵的，因为那意味着我能够很容易的得到特征值，而特征值计算对应一个一元高次方程，而高于5次的方程是没有代数解法的。从这个角度说求特征值比用高斯消元解Ax=b（得到上对角阵U）、求逆矩阵等更困难。我们通过QR方法，是尽可能逼近特征值。

因此求取特征值的第一步是通过相似变换得到如上图的上Hessenburg矩阵，它由一个上三角阵和一条下对角线组成。

第二步是通过带有平移操作的QR分解得到特征值。

这也是MATLAB中eig(A)的操作内容，你可以从线代数据包LAPACK中获得代码，了解数值计算过程。

如果矩阵A0是对称的，会发现A1也对称，最后会得到“对称的Hessenburg矩阵”，即三对角矩阵。矩阵不是n的平方个数组成的，只有2n个数，因为上对角线和下对角线是相同的。

这就是eig方法求特征值的过程，这算法的真正核心是QR，不要再试着去求解行列式 \[\det (\boldsymbol A-\lambda \boldsymbol I)=0\] ，它把n平方个信息压缩到n个系数里，所以你会丢失很多信息。

奇异值计算

下面讨论奇异值，不要一上来就用 \[{{\boldsymbol A}^{T}}\boldsymbol A\] 进行计算。前面已经讨论过，我们可以将对称矩阵S通过相似矩阵 \[\boldsymbol Q\boldsymbol S{{\boldsymbol Q}^{T}}=\boldsymbol Q\boldsymbol S{{\boldsymbol Q}^{-1}}\] 变成一个三对角阵，而不改变其特征值，而奇异值有没有对应的处理方法？

矩阵的SVD分解为 \[\boldsymbol A=\boldsymbol U\boldsymbol \varSigma {{\boldsymbol V}^{T}}\] ，会发现\[\boldsymbol Q\boldsymbol A=\boldsymbol Q\boldsymbol U\boldsymbol \varSigma {{\boldsymbol V}^{T}}\]不会改变矩阵的奇异值，因为左侧两个正交矩阵的乘积QU依旧是正交矩阵。

\[{{(\boldsymbol Q\boldsymbol U)}^{-1}}={{\boldsymbol U}^{-1}}{{\boldsymbol Q}^{-1}}={{\boldsymbol U}^{T}}{{\boldsymbol Q}^{T}}={{(\boldsymbol Q\boldsymbol U)}^{T}}\] 可验证QU是正交矩阵。

而SVD分解的右侧可以乘以任意正交矩阵，而并非Q矩阵的转置，因为奇异值矩阵没有发生变化，而两侧仍是正交矩阵，即 \[{{\boldsymbol Q}_{1}}\boldsymbol A{{\boldsymbol Q}_{2}}={{\boldsymbol Q}_{1}}\boldsymbol U\boldsymbol \varSigma {{\boldsymbol V}^{T}}{{\boldsymbol Q}_{2}}\] 。这比之前求特征值的方法有了更大的灵活性，最终矩阵A可以被处理为双对角线矩阵。

而对于\[{{\boldsymbol A}^{T}}\boldsymbol A\]，如果直接处理得到的就是三对角矩阵，如果将矩阵A变换得到的二对角矩阵与其转置相乘也得到三对角矩阵。

两套处理方法都是先通过相似变换得到具有很多零的稀疏矩阵，然后再通过带平移操作的QR分解得到特征值或者奇异值。对于1000阶以下的矩阵，这些就是eig和SVD命令的基本操作。

Krylov方法

当n太大时（例如矩阵A的大小为100万），需要用Krylov方法。从向量b开始，将其乘以A得到Ab……然后一直到 \[{{\boldsymbol A}^{99}}\bold {b}\] ，得到100维的Krylov空间。它将矩阵限制在这个100维的空间中，而这个子空间就可以捕获特征向量。向量v表示为线性组合：

\[\bold v={{c}_{1}}\bold b+{{c}_{2}}\boldsymbol A\bold b+{{c}_{3}}{\boldsymbol {A}^{2}}\bold b+\cdots +{{c}_{100}}{\boldsymbol {A}^{99}}\bold b+error\]

这个误差是可以忽略的。当向量乘以A时，所得向量有存在于Krylov空间k100中的部分，也有在k100之外的部分，而我们忽略了后一部分。

我只是花几分钟来解释Krylov的想法可以做些什么。通过构造Krylov向量得到这种特定类型的子空间，通过Gram-schmidt快速获得它的基向量，考察一下矩阵在该空间中的作用，可以寻找到限于该空间的特征值。所以转化为100 x 100的问题（矩阵A通过100个Krylov空间中的正交向量变换），然后找到该矩阵的特征值，它们是原特征值的一个很好的近似。但是我可能不确定这个想法会给前100个特征值，它们不是完美的特征值，但有一定精度。

下一步将讨论随机抽样。如果您的矩阵太大了，需要进行随机抽样，这是一个非常新的想法，是数值线性代数中的一个非常不同的想法。