大型稀疏矩阵的所有特征向量都为零
All eigenvectors of large sparse matrices are zero
我有一个 50,000 x 50,000 或更大的密集矩阵。如果我使用 numpy 或 scipy- 包，我所有特征向量的条目都是 0。如果我使用 scipy.sparse 来计算仅 1000-8000 个特征向量，我得到正确的特征向量。但是我都需要。
这会不会是内存问题？或者出现这样的问题是什么原因？
我可以使用 LAPACK 或 ARPACK 来计算正确的特征向量吗？
请注意我的矩阵是 networkx 图的表示，因此是稀疏矩阵。我将它们转换为密集矩阵以使用 numpy.linalg，否则我使用 scipy.sparse.linalg.
关于numpy.linalg.eig() and scipy.linalg.eig()的问题可能是与矩阵大小有关的内存错误：50000x50000双精度密集矩阵占用18Go内存。
numpy.linalg.eig() and scipy.linalg.eig()依赖LAPACK的套路DGEEV()。 LAPACK DGEEV() 和 DGEEVX() 使用 QR 算法 来计算 [=58= 的所有特征值和特征向量]密集矩阵。首先，使用 dgehrd() 将矩阵简化为上 Hessenberg 形式，然后在 dhseqr() 中执行 QR 迭代。在DGEEVX()中，首先对矩阵进行平衡和缩放。
另一方面，scipy.sparse.linalg.eigs() and scipy.sparse.linalg.eigsh() rely on ARPACK functions implementing the Implicitly Restarted Arnoldi Method and Implicitly Restarted Lanczos Method 在稀疏矩阵上。两者都是幂法的改进，在计算最大的 eigenvalues/eigenvectors 时非常有效，精度更高。如果 Ax=b 可以快速求解，这些迭代方法在找到给定值附近的最小 eigenvalues/eigenvectors 或 eigenvalues/eigenvectors 时也非常有效。
这些方法之间的区别在 Lloyd N. Trefethen and David Bau, NUMERICAL LINEAR ALGEBRA, Lecture 33. The Arnoldi Iteration 中有解释。
... whereas the Arnoldi iteration is based upon the QR factorization (33.7) of the matrix ,whose columns are b, A b, ... ,A^{n-1} b, simultaneous iteration and the QR algorithm are based upon the QR factorization (28.16) of the matrix whose columns are A^n e_1 ... , A^n e_n .
从实际的角度来看，Arnoldi 迭代总是用于检索显着小于矩阵大小的有限数量的特征向量。然而，scipy.sparse.linalg.eigs() or scipy.sparse.linalg.eigsh() enables finding interior eigenvalues and eigenvectors near sigma. Hence, it is possible to call scipy.sparse.linalg.eigsh()的参数sigma多次使用不同的sigma。如果特征值都具有有限的多重性，则可以恢复所有特征向量。通过分离特征值并跟踪它们的多重性来避免可能的重复。
使用 sigma 的基本调用写道：
sparsevalue, sparsevector = sp.sparse.linalg.eigsh(M.copy(), k = kkeig,sigma=sig, which='LM', mode='normal')
如果矩阵是 Hermitian semi
这是一个示例代码，基于您之前删除的 post，它计算半正定稀疏矩阵的所有特征值。由于特征值都是实数且为正，因此 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()
对于具有实数正负特征值的厄尔米特矩阵，sigma 的正值和负值都有待探索。如果矩阵不是 Hermitian 矩阵，则特征值可能不是实数，并且要选择复平面上 sigma 的值。首先搜索 A 的最大特征值的大小，将区域限制为圆盘。
建议的方法非常慢，可能并不总是有效。它对 20000x20000 矩阵工作一次，使用 1Go 内存。
我有一个 50,000 x 50,000 或更大的密集矩阵。如果我使用 numpy 或 scipy- 包，我所有特征向量的条目都是 0。如果我使用 scipy.sparse 来计算仅 1000-8000 个特征向量，我得到正确的特征向量。但是我都需要。
这会不会是内存问题？或者出现这样的问题是什么原因？ 我可以使用 LAPACK 或 ARPACK 来计算正确的特征向量吗？
请注意我的矩阵是 networkx 图的表示，因此是稀疏矩阵。我将它们转换为密集矩阵以使用 numpy.linalg，否则我使用 scipy.sparse.linalg.
关于numpy.linalg.eig() and scipy.linalg.eig()的问题可能是与矩阵大小有关的内存错误：50000x50000双精度密集矩阵占用18Go内存。
numpy.linalg.eig() and scipy.linalg.eig()依赖LAPACK的套路DGEEV()。 LAPACK DGEEV() 和 DGEEVX() 使用 QR 算法 来计算 [=58= 的所有特征值和特征向量]密集矩阵。首先，使用 dgehrd() 将矩阵简化为上 Hessenberg 形式，然后在 dhseqr() 中执行 QR 迭代。在DGEEVX()中，首先对矩阵进行平衡和缩放。
另一方面，scipy.sparse.linalg.eigs() and scipy.sparse.linalg.eigsh() rely on ARPACK functions implementing the Implicitly Restarted Arnoldi Method and Implicitly Restarted Lanczos Method 在稀疏矩阵上。两者都是幂法的改进，在计算最大的 eigenvalues/eigenvectors 时非常有效，精度更高。如果 Ax=b 可以快速求解，这些迭代方法在找到给定值附近的最小 eigenvalues/eigenvectors 或 eigenvalues/eigenvectors 时也非常有效。
这些方法之间的区别在 Lloyd N. Trefethen and David Bau, NUMERICAL LINEAR ALGEBRA, Lecture 33. The Arnoldi Iteration 中有解释。
... whereas the Arnoldi iteration is based upon the QR factorization (33.7) of the matrix ,whose columns are b, A b, ... ,A^{n-1} b, simultaneous iteration and the QR algorithm are based upon the QR factorization (28.16) of the matrix whose columns are A^n e_1 ... , A^n e_n .
从实际的角度来看，Arnoldi 迭代总是用于检索显着小于矩阵大小的有限数量的特征向量。然而，scipy.sparse.linalg.eigs() or scipy.sparse.linalg.eigsh() enables finding interior eigenvalues and eigenvectors near sigma. Hence, it is possible to call scipy.sparse.linalg.eigsh()的参数sigma多次使用不同的sigma。如果特征值都具有有限的多重性，则可以恢复所有特征向量。通过分离特征值并跟踪它们的多重性来避免可能的重复。
使用 sigma 的基本调用写道：
sparsevalue, sparsevector = sp.sparse.linalg.eigsh(M.copy(), k = kkeig,sigma=sig, which='LM', mode='normal')
如果矩阵是 Hermitian semi
这是一个示例代码，基于您之前删除的 post，它计算半正定稀疏矩阵的所有特征值。由于特征值都是实数且为正，因此 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()
对于具有实数正负特征值的厄尔米特矩阵，sigma 的正值和负值都有待探索。如果矩阵不是 Hermitian 矩阵，则特征值可能不是实数，并且要选择复平面上 sigma 的值。首先搜索 A 的最大特征值的大小，将区域限制为圆盘。
建议的方法非常慢，可能并不总是有效。它对 20000x20000 矩阵工作一次，使用 1Go 内存。