我有一个 50,000 x 50,000 或更大的密集矩阵。如果我使用 numpy 或 scipy- 包，我所有的特征向量的条目都是 0。如果我使用 scipy.sparse 只计算 1000-8000 个特征向量，我会得到正确的特征向量。但我需要所有这些。
这可能是内存问题吗？或者这个问题的原因是什么？ 我可以使用 LAPACK 或 ARPACK 来计算正确的特征向量吗？
请注意，我的矩阵是 networkx 图的表示，因此是稀疏矩阵。我将它们转换为密集矩阵以供使用numpy.linalg，否则我使用scipy.sparse.linalg.
关于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 内存。