矩阵最大特征值的求法
Largest Eigenvalue of Matrix
矩阵的特征值分解是在科学计算常见的线性代数方法。在一些情况下，我们不需要全部特征值，只需要一部分特征值。或者因为矩阵的一些性质，我们只对一部分特征值感兴趣。比如对于矩阵 B = A’ * A， 如果A是m by n的实矩阵，m很大，n很小。根据线性代数理论，A, B矩阵的Rank最多是n，并且B矩阵是半正定矩阵。对矩阵B做特征值分解时，其特征值都是非负实数，并且最多有n个正数。因此，如果求解B矩阵的特征值，只需要知道其最大的n个特征值，因为其余特征值都是零。
求解一般矩阵的特征值主要有两种：（1）特征多项式（2）Power Method（包括 Arnoldi方法和QR分解)。 对于实对称正定矩阵，可以使用Choleskey分解。这些方法可以求出所有的特征值。而对于我们特别关心的实对称(Hermitain/Symmetric)正定(Positive Definite)矩阵，比如协方差矩阵，这里我们关心的是其最大（或者最小）的特征值。这里常用的方法的方法叫做Linczos方法，可以看做是Arnoldi方法的特例。Linczos方法计算A的特征值是（1）先随机生成一个向量v，（2）再利用一系列向量： v, A*v, A*A*v, …， A^(n-1) * v 来得到一个特殊矩阵T，这个矩阵的大小是m by m，这里m往往远大于n。（3）得到T之后，可以以O(m^2)的复杂度来的到特征值。
实现Lanczos方法时会使用Multiple restart，最有名的实现是ARPACK。ARPACK是用Fortran实现的，对Fortran不熟悉就不容易直接调用ARPACK的函数。不过很多其他软件都对ARPACK的功能进行了封装。在R里面使用ARPACK，可以用igraph包，R-blogger有一个这方面的例子（链接）。其他软件中一般也都有用ARPACK来实现求一部分特征值的功能：Matlab的eigs，Scipy的eigs和eighs（分别对应实对称矩阵和Hermitian矩阵， 链接），Julia的eigs（链接）。
后话：现在开源的科学计算软件很多，当我们想寻找合适高效的数值计算工具时，可以从源代码出发。比如我们知道Matlab有eigs函数，就可以用“julia eigs”来看julia语言里是怎么实现eigs的。