关于numpy.linalg.eig()和scipy.linalg.eig()的问题可能是与矩阵大小有关的内存错误：50000x50000双精度密集矩阵占用18Go内存。
numpy.linalg.eig() 和 scipy.linalg.eig() 依赖于 LAPACK 的例程 DGEEV()。 LAPACK DGEEV() 和 DGEEVX() 使用 QR algorithm 计算 密集矩阵 的所有特征值和特征向量。首先，使用dgehrd() 将矩阵简化为上 Hessenberg 形式，然后在dhseqr() 中执行 QR 迭代。在DGEEVX() 中，首先对矩阵进行平衡和缩放。
另一方面，scipy.sparse.linalg.eigs() 和 scipy.sparse.linalg.eigsh() 依赖于在稀疏矩阵上实现 Implicitly Restarted Arnoldi Method 和 Implicitly Restarted Lanczos Method 的 ARPACK 函数。两者都是幂方法的改进，并且在计算最大特征值/特征向量时非常有效，并且精度更高。如果 Ax=b 可以快速求解，这些迭代方法在找到最小特征值/特征向量，或给定值附近的特征值/特征向量方面也非常有效。
Lloyd N. Trefethen and David Bau, NUMERICAL LINEAR ALGEBRA, Lecture 33. The Arnoldi Iteration 中解释了这些方法之间的区别。
...而Arnoldi迭代基于矩阵的QR分解（33.7），其列是b，A b，...，A^{n-1} b，同时迭代和QR算法基于矩阵的 QR 分解 (28.16)，其列为 A^n e_1 ... , A^n e_n 。
从实际的角度来看，Arnoldi 迭代总是用于检索有限数量的特征向量，该特征向量明显低于矩阵的大小。然而，scipy.sparse.linalg.eigs() 或scipy.sparse.linalg.eigsh() 的参数sigma 可以在sigma 附近找到内部特征值和特征向量。因此，可以使用不同的sigma 多次调用scipy.sparse.linalg.eigsh()。 如果特征值都具有有限的多重性，则可以恢复所有特征向量。通过分离特征值并跟踪它们的多重性来避免潜在的重复。
使用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 内存。