python - 大型稀疏矩阵的所有特征向量都为零

关于numpy.linalg.eig()并且scipy.linalg.eig()可能是与矩阵大小相关的内存错误的问题：50000x50000 双精度密集矩阵占用 18Go 内存。

numpy.linalg.eig()并且scipy.linalg.eig()依赖DGEEV()LAPACK的套路。LAPACKDGEEV()并DGEEVX()利用QR 算法计算密集矩阵的所有特征值和特征向量。首先，使用将矩阵简化为上 Hessenberg 形式dgehrd()，然后在中执行 QR 迭代dhseqr()。在DGEEVX()中，首先对矩阵进行平衡和缩放。

另一方面，scipy.sparse.linalg.eigs()依靠scipy.sparse.linalg.eigsh()ARPACK 函数在稀疏矩阵上实现隐式重启 Arnoldi 方法和隐式重启 Lanczos 方法。两者都是幂方法的改进，并且在计算最大特征值/特征向量时非常有效，并且精度更高。如果可以快速求解 Ax=b，则这些迭代方法在找到给定值附近的最小特征值/特征向量或特征值/特征向量方面也非常有效。

Lloyd N. Trefethen 和 David Bau, NUMERICAL LINEAR ALGEBRA, Lecture 33. The Arnoldi Iteration解释了这些方法之间的区别。

...而Arnoldi迭代基于矩阵的QR分解（33.7），其列是b，A b，...，A^{n-1} b，同时迭代和QR算法基于列为 A^n e_1 ... , A^n e_n 的矩阵的 QR 分解 (28.16)。

从实际的角度来看，Arnoldi 迭代总是用于检索有限数量的特征向量，该特征向量明显低于矩阵的大小。然而，或的 论点能够在附近找到内部特征值和特征向量。因此，可以使用不同的. sigmascipy.sparse.linalg.eigs()scipy.sparse.linalg.eigsh()sigmascipy.sparse.linalg.eigsh()sigma如果特征值都具有有限的多重性，则可以恢复所有特征向量。通过分离特征值并跟踪它们的多重性来避免潜在的重复。

使用sigma写入的基本调用：

sparsevalue, sparsevector = sp.sparse.linalg.eigsh(M.copy(), k = kkeig,sigma=sig, which='LM', mode='normal')

如果矩阵是 Hermitian semi 这是一个示例代码，基于您之前删除的帖子，它计算正半正定稀疏矩阵的所有特征值。由于特征值都是实数和正数，sigma逐渐增加以找到所有特征值：

import numpy as np
import networkx as nx
import scipy as sp




def tolist(sparsevalue, keeplast=False):
    listofeigval=[]
    listofmultiplicity=[]

    curreig=sparsevalue[0]-1.1*np.abs(sparsevalue[0])
    for i in range(len(sparsevalue)):
           #print curreig, sparsevalue[i], len(listofeigval)
           if np.abs(curreig-sparsevalue[i])<1e-11*(np.abs(curreig)+np.abs(sparsevalue[i])):
                listofmultiplicity[-1]=listofmultiplicity[-1]+1
           else :
                curreig=sparsevalue[i]
                listofeigval.append(curreig)
                listofmultiplicity.append(1)

    if keeplast==False:
        #the last one is not sure regarding multiplicity:
        listofeigval.pop()
        listofmultiplicity.pop()

    return listofeigval,listofmultiplicity

def test():
    N_1 = 2000
    R_1 = 0.1
    k = 0
    iterations = 1
    while k < iterations:
      G = nx.random_geometric_graph(N_1, R_1)
      if nx.is_connected(G) == True:
          print 'got connected network'
          k = k+1
          M=nx.laplacian_matrix(G)      #M is here a sparse matrix
          M = M.astype(float)
          #M[0,0]=M[0,0]+1.  # this makes the laplacian_matrix positive definite.
          #sparsevalue, sparsevector = sp.sparse.linalg.eigsh(M, k = N_1-2)
          kkeig=400
          sparsevalue, sparsevector = sp.sparse.linalg.eigsh(M.copy(), k = kkeig, which='SM')
          print sparsevalue

          listofeigval=[]
          listofmultiplicity=[]
          listofeigval,listofmultiplicity=tolist(sparsevalue)
          print len(listofeigval), len(listofmultiplicity)

          nbeigfound=0
          for mul in listofmultiplicity:
              nbeigfound=nbeigfound+mul
          keepgoing=True
          while( nbeigfound<N_1):
               print '(',nbeigfound,'/',N_1,')  is ', listofeigval[-1]
               meanspacing=0.
               meanspacingnb=0.

               for i in range(10):
                     meanspacingnb=meanspacingnb+listofmultiplicity[len(listofeigval)-i-1]
               meanspacing=(listofeigval[-1]-listofeigval[len(listofeigval)-10])/meanspacingnb
               sig=listofeigval[-1]+0.1*kkeig*meanspacing

               keeplast=False
               if nbeigfound<N_1-0.5*kkeig:
                   sparsevalue, sparsevector = sp.sparse.linalg.eigsh(M.copy(), k = kkeig,sigma=sig, which='LM', mode='normal')
               else:
                   keeplast=True
                   sparsevalue, sparsevector = sp.sparse.linalg.eigsh(M.copy(), k = kkeig, which='LM')
               listofneweigval,listofnewmultiplicity=tolist(sparsevalue,keeplast)
               #need to retreive the starting point
               index=len(listofeigval)-2*kkeig/10
               if listofneweigval[1]>listofeigval[index]:
                   while(np.abs(listofneweigval[1]-listofeigval[index])>1e-10*(np.abs(listofneweigval[1])+np.abs(listofeigval[index]))):
                       index=index+1
               else:
                   while(np.abs(listofneweigval[1]-listofeigval[index])>1e-10*(np.abs(listofneweigval[1])+np.abs(listofeigval[index]))):
                       index=index-1
               #The first matching eigenvalue is found.
               #zipping the list and checking if it works.
               i=1
               while index<len(listofeigval) and i<len(listofneweigval) :
                    if (np.abs(listofneweigval[i]-listofeigval[index])>1e-10*(np.abs(listofneweigval[i])+np.abs(listofeigval[index]))):
                         print 'failed after ', index, ' different eigenvalues', ': wrong eigenvalue'
                         return listofeigval[0:index], listofmultiplicity[0:index], 1
                    if listofmultiplicity[index] != listofnewmultiplicity[i] :
                         print 'failed after ', index, ' different eigenvalues', ': wrong multiplicity'
                         return listofeigval[0:index], listofmultiplicity[0:index], 2
                    index=index+1
                    i=i+1
               #adding the remaining eigenvalues.
               while i<len(listofneweigval) :
                    listofeigval.append(listofneweigval[i])
                    listofmultiplicity.append(listofnewmultiplicity[i])
                    nbeigfound=nbeigfound+listofnewmultiplicity[i]
                    i=i+1
          print 'number of eigenvalues: ', nbeigfound
          nbl=0
          for i in range(len(listofeigval)):
               print 'eigenvalue ',i,' (',nbl,'/',N_1,')  is %14.8f'% listofeigval[i], ' with multiplicity', listofmultiplicity[i]
               nbl=nbl+listofmultiplicity[i]

          return listofeigval, listofmultiplicity, 0


          #sig=39.1
          #sparsevalue, sparsevector = sp.sparse.linalg.eigsh(M.copy(), k = 100,sigma=sig, which='LM', mode='normal')
          #print sparsevalue

test()

对于具有真实正负特征值的 Hermitian 矩阵，需要探索 sigma 的正负值。如果矩阵不是 Hermitian，则特征值可能不是实数，sigma需要选择复平面上的值。首先搜索最大特征值的大小，A将区域限制为一个圆盘。

建议的方法非常慢，可能并不总是有效。它为 20000x20000 矩阵工作了一次，使用 1Go 内存。

python - 大型稀疏矩阵的所有特征向量都为零

问题描述

解决方案

推荐阅读