Sep 082013
 

矩阵最大特征值的求法

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的。

 

 

 Leave a Reply

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

(required)

(required)

*