几个特征值求解函数的性能比较

发表于2016-06-25由范雨

正如在scipy和numpy中与特征值问题相关的函数谈到的，Python的常用数值计算库numpy和scipy中一共有好几个功能和限制稍有不同的特征值问题求解函数，今天我们用一个8×8的“特殊”矩阵（Rosser矩阵）来看看这些函数的计算速度和结果精度如何。

Rosser矩阵

Rosser矩阵是这样的完全由整数构成的，对称的8×8矩阵：

\begin{bmatrix} 611 &196 & -192& 407& -8& -52& -49& 29&\\196& 899& 113& -192& -71& -43& -8& -44&\\-192& 113& 899& 196& 61& 49& 8& 52&\\407& -192& 196& 611& 8& 44& 59& -23&\\-8& -71& 61& 8& 411& -599& 208& 208&\\-52& -43& 49& 44& -599& 411& 208& 208&\\-49& -8& 8& 59& 208& 208& 99& -911&\\29& -44& 52& -23& 208& 208& -911& 99&\end{bmatrix}

可以通过符号运算得出该矩阵的精确特征值，见如下python代码：


import numpy as np
import sympy as sy
m = np.loadtxt('Rosser.txt')
print m
R = sy.Matrix(m)
re = R.eigenvals()
for key in re.keys():
 print '%30s,%5s,\t%f'%(key,re[key],key.evalf())

运行代码后，得：

                             0,    1,	0.000000
                          1000,    2,	1000.000000
                10*sqrt(10405),    1,	1020.049018
            100*sqrt(26) + 510,    1,	1019.901951
               -10*sqrt(10405),    1,	-1020.049018
                          1020,    1,	1020.000000
           -100*sqrt(26) + 510,    1,	0.098049

每一行的三个数据分别是精确的特征值表达式，特征值重数，和精确表达式的浮点数。可以发现该矩阵的特征值分布非常棘手：有1个零特征值，1个靠近零的特征值，1个值为1000的2重根，三个都很靠近1200的根（不是三重根），最大正特征值和最小负特征值的量级相同。这些特点使得它通常用来测试特征值算法的精确性。要用书本上最简单的方法来计算它的特征值确实是一个问题，但是对通用的数值计算库来说并不是问题，比如这里给出了MATLAB中的结果。

测试代码

import numpy
import scipy.linalg
import scipy.sparse.linalg
import time
import matplotlib.pyplot as plt

m = numpy.loadtxt('Rosser.txt')
print m

funcdic = {0:('numpy.linalg.eig',numpy.linalg.eig,{}),
           1:('numpy.linalg.eigh',numpy.linalg.eigh,{}),
           2:('scipy.linalg.eig',scipy.linalg.eig,{}),
           3:('scipy.linalg.eigh',scipy.linalg.eigh,{}),
           4:('scipy.sparse.linalg.eigs',scipy.sparse.linalg.eigs,{}),
           5:('scipy.sparse.linalg.eigsh',scipy.sparse.linalg.eigsh,{}),
           6:('scipy.sparse.linalg.eigs, sigma=0',scipy.sparse.linalg.eigs,{'sigma':0}),
           7:('scipy.sparse.linalg.eigsh, sigma=0',scipy.sparse.linalg.eigsh,{'sigma':0})}

def timing_val(func):
    def wrapper(*arg, **kw):
        t1 = time.time()
        for i in range(0,1000):
            res = func(*arg, **kw)
        t2 = time.time()
        print '%40s  %.3e sec' % (funcdic[arg[0]][0], t2-t1)
        return [res,t2-t1]
    return wrapper

@timing_val
def test_func(i):
    name,func,args = funcdic[i]
    eigval,eigvec = func(m,**args)
    return numpy.real(numpy.sort(eigval))

def compare():
    Resultslist = []
    for i in range(0,8):
        try:
            Resultslist.append(test_func(i))
        except Exception as e:
            print type(e)
            print e
        continue
    return Resultslist

if __name__ == '__main__':
    Resultslist = compare()

    exact = [-1020.049018,0,0.098049,1000,1000,1019.901951,1020.000000,1020.049018]

    for i in range(0,8):
        name,func,args = funcdic[i]
        print '='*20
        print name
        for j in range(0,8):
            if j < len(Resultslist[i][0]):
                print '%20f,%20f'%(exact[j],Resultslist[i][0][j])
            else:
                print '%20f,'%(exact[j],)
        print '='*20
        print '\n\n'

    info = 'method index - details\n'
    for i in range(0,8):
        info = info + '%d - %s\n'%(i,funcdic[i][0])

    plt.figure()
    tt = [Resultslist[i][1] for i in range(0,8)]
    plt.bar(numpy.arange(0,8)-0.4, tt, width=0.8,alpha=0.6)
    plt.xticks(range(0,8),[str(i) for i in range(0,8)])
    plt.xlabel('method index',fontsize = 16)
    plt.ylabel('time (s)',fontsize = 16)
    plt.text(0,0.5,info,fontsize = 14)
    plt.show()

计算结果对比

运行测试代码，在控制台输出：

====================
numpy.linalg.eig
        -1020.049018,        -1020.049018
            0.000000,           -0.000000
            0.098049,            0.098049
         1000.000000,         1000.000000
         1000.000000,         1000.000000
         1019.901951,         1019.901951
         1020.000000,         1020.000000
         1020.049018,         1020.049018
====================

====================
numpy.linalg.eigh
        -1020.049018,        -1020.049018
            0.000000,           -0.000000
            0.098049,            0.098049
         1000.000000,         1000.000000
         1000.000000,         1000.000000
         1019.901951,         1019.901951
         1020.000000,         1020.000000
         1020.049018,         1020.049018
====================

====================
scipy.linalg.eig
        -1020.049018,        -1020.049018
            0.000000,           -0.000000
            0.098049,            0.098049
         1000.000000,         1000.000000
         1000.000000,         1000.000000
         1019.901951,         1019.901951
         1020.000000,         1020.000000
         1020.049018,         1020.049018
====================

====================
scipy.linalg.eigh
        -1020.049018,        -1020.049018
            0.000000,            0.000000
            0.098049,            0.098049
         1000.000000,         1000.000000
         1000.000000,         1000.000000
         1019.901951,         1019.901951
         1020.000000,         1020.000000
         1020.049018,         1020.049018
====================

====================
scipy.sparse.linalg.eigs
        -1020.049018,        -1020.049018
            0.000000,         1000.000000
            0.098049,         1000.000000
         1000.000000,         1019.901951
         1000.000000,         1020.000000
         1019.901951,         1020.049018
         1020.000000,
         1020.049018,
====================

====================
scipy.sparse.linalg.eigsh
        -1020.049018,        -1020.049018
            0.000000,         1000.000000
            0.098049,         1000.000000
         1000.000000,         1019.901951
         1000.000000,         1020.000000
         1019.901951,         1020.049018
         1020.000000,
         1020.049018,
====================

====================
scipy.sparse.linalg.eigs, sigma=0
        -1020.049018,        -1000.694183
            0.000000,            0.000000
            0.098049,            0.098048
         1000.000000,         1000.000000
         1000.000000,         1020.029187
         1019.901951,         1020.133333
         1020.000000,
         1020.049018,
====================

====================
scipy.sparse.linalg.eigsh, sigma=0
        -1020.049018,            0.000000
            0.000000,            0.098023
            0.098049,           57.341793
         1000.000000,         1000.070411
         1000.000000,         1007.980761
         1019.901951,         1019.958703
         1020.000000,
         1020.049018,
====================

其中第一列是精确的特征值，第二行是对应方法得出的结果。可以发现numpy.linalg 和scipy.linalg中的函数（最开始的4组结果）都完美地计算出了所有特征值.对于scipy.sparse.linalg中的函数，由于采用子空间迭代法，本质上无法求解出所有特征值，所有这次采用默认设置，求解出6个特征值。可以发现采用默认参数的scipy.sparse.linalg.eigs和scipy.sparse.linalg.eigsh虽然没能计算出所有特征值，但是给出的6个特征值都是精确的。不足是未能求出靠近零的那两个特征值，为了求出它们需要通过sigma参数设置偏移量（最后两组结果）。这样，虽然算出了靠近零的那两个特征值，但是其他特征值的结果出现了误差。对于scipy.sparse.linalg.eigs 函数，算出了不存在的-1000.69，而scipy.sparse.linalg.eigsh和居然算出了57.3这个偏差很大的特征值！

计算时间对比

figure_1

讨论

在运算时间上（比较计算1000次的累积时间），scipy.sparse库中的函数明显要慢于scipy.linalg和numpy.linalg的函数。各个库中名字带’h’的函数都比不带的快，这不难理解，因为有’h’的函数考虑了矩阵的对称性，加快了计算。一个有意思的结果是numpy.linalg中的函数居然要略快于scipy.linalg，要知道在scipy官方网站的解释中可都是推荐scipy.linalg的。这里的具体的原因不明，可能是一些函数调用过程中的时间损耗。可以用随机数生成一个更大的矩阵来做运算，这两者的差异很小。

根据各个函数在Rosser矩阵上的表现，说明如果矩阵规模较小，而又需要几乎所有的特征值的话，还是不要用稀疏矩阵库的函数了。他们无论从求解效率和精度来说，都不占优势。当然，当矩阵规模很大而需要的仅仅是前几阶特征值的时候，还是应当用稀疏矩阵库的函数的。

未知 的头像

About 范雨

Talk is cheap. Show me the code.
此条目发表在特征值问题, Dynamics, Python, 数值计算分类目录。将固定链接加入收藏夹。

留下评论