您的位置:首页 > 其它

scipy.sparse求稀疏矩阵前k个特征值

2012-12-13 11:46 302 查看
scipy.sparse求稀疏矩阵前k个特征值 - Waleking的专栏 - 博客频道 - CSDN.NET

scipy.sparse求稀疏矩阵前k个特征值

分类: python 机器学习 2012-05-30 13:03 448人阅读 评论(0) 收藏 举报 背景:

要在python中处理7000*7000的稀疏矩阵,计算前k小的特征值和相应的特征向量。不想在matlab中做这件事了,所有的数据预处理和展现工作都想在python中完成。然而一般的linalg提供的eig开销太大,要计算所有的特征值和特征向量,这个开销要达到 O(N^3),对于谱聚类来说,这个开销是不能忍受的。

所以要借助稀疏矩阵计算的工具包。

探索过程:

使用scipy.sparse

对于10*10大小的矩阵都没有问题:

[python] view plaincopyprint?import scipy as sp import scipy.sparse.linalg import bumpy as np vals, vecs = sp.sparse.linalg.eigs(np.identity(10), k=6)
import scipy as sp
import scipy.sparse.linalg
import bumpy as np
vals, vecs = sp.sparse.linalg.eigs(np.identity(10), k=6)


但是对于100*100的矩阵就出错了:

[python] view plaincopyprint?vals, vecs = sp.sparse.linalg.eigs(np.identity(100), k=6)
vals, vecs = sp.sparse.linalg.eigs(np.identity(100), k=6)


过程当中总是出现错误:

ArpackError: ARPACK error 3: No shifts could be applied during a cycle of the Implicitly restarted Arnoldi iteration. One possibility is to increase the size of NCV relative to NEV

查看了一下源代码,发现求特征值的方法引用了:

ARPACK USERS GUIDE: Solution of Large Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods

文章在这里 http://people.sc.fsu.edu/~jburkardt/pdf/arpack.pdf scipy.sparse.linalg引用它大意是说使用了Arnoldi方法求大规模矩阵的特征值问题。

查阅了一下手册,那个ArpackError错误的意思是要修改参数。

既然是这样,就查看一下scipy.sparse的手册吧 http://docs.scipy.org/doc/scipy/reference/sparse.html
里面介绍了几种类型的稀疏矩阵,sp.sparse.linalg.eigs的参数也应该是稀疏矩阵,于是修改了代码:

[python] view plaincopyprint?import scipy as sp import scipy.sparse.linalg import numpy as np import time A = sp.sparse.lil_matrix((10000, 10000)) A[0,:1000]=np.random.rand(1000) A.setdiag(np.random.rand(10000)) timeCheckin=time.clock() vals, vecs = sp.sparse.linalg.eigs(A, k=3) print("first 6 eigenvaluse cost %s" % (time.clock()-timeCheckin)) print vals
import scipy as sp
import scipy.sparse.linalg
import numpy as np
import time
A = sp.sparse.lil_matrix((10000, 10000))
A[0,:1000]=np.random.rand(1000)
A.setdiag(np.random.rand(10000))
timeCheckin=time.clock()
vals, vecs = sp.sparse.linalg.eigs(A, k=3)
print("first 6 eigenvaluse cost %s" % (time.clock()-timeCheckin))
print vals

构造了一个10000*10000的矩阵,第一行,第二行,第四行有 [0,1]上均匀分布的随机数。计算前3个特征值,耗时52.091 s. 又做了一次实验,因随机性,这次耗时314.156s,初步判断是和算法的收敛状况有关。

需要说明的是,如果对稀疏矩阵直接用矩阵的乘法运算,如dot(A,vecs[:,0]) 就会非常悲剧,貌似稀疏矩阵会完全展开成普通矩阵计算,在mac上运行了10个小时,内存最高峰时使用2.84GB,即便这样还是没有得到结果(不过证明了64位python的鲁棒性)。

另外使用MATLAB做对比,matlab代码是

[plain] view plaincopyprint?d=rand(10000,1) A=diag(d) A[1,1:1000]=rand(1,1000) eigs(A,3)
d=rand(10000,1)
A=diag(d)
A[1,1:1000]=rand(1,1000)
eigs(A,3)
查手册,是用Ritz方法,迭代了965次收敛,耗时30min。

但如果也让MATLAB对稀疏矩阵做操作,那么它也调用ARPACK工具包,也就是说和scipy.sparse是一样的,精度上没有区别,计算速度也和python的相当。

另外可以考虑pysparse:

这里是PySparse的官方网站
http://pysparse.sourceforge.net/spmatrix.html#vectorization
最终解决方案:

[python] view plaincopyprint?#coding=utf8
import scipy as sp
import scipy.sparse.linalg
import scipy.io as sio
import numpy as np
import time
import os
def generateTestMat(size):
'''''generate a test matrix for k-eigenvalues problem
' matlab and scipy.sparse.linalg.eigs seem use the same package ARPACK
' this matrix has two diagnol lines, one is on the diagnol and the other is off 1
'''
A = sp.sparse.lil_matrix((size, size))
A[0,:size]=np.random.rand(size)
A.setdiag(np.random.rand(size))
A.setdiag(np.random.rand(size-1),1)
#保存成matlab可以读取的格式,{"A":A}前面的A表示在matlab中的名字,后面的A表示在python中名字
sio.savemat("A.mat",{"A":A},oned_as='row')
print("generate a test matrix,with size %s by %s" %(size,size))

def calculateKEigs(tolerance=0):
if(os.path.exists(os.getcwd()+"/"+"A.mat")==False):
generateTestMat()
A=sio.loadmat("A.mat")["A"]#后面的A表示在matlab中名字
timeCheckin=time.clock()
#tol=0表示使用原先的计算精度
vals, vecs = sp.sparse.linalg.eigs(A, k=3,tol=tolerance,which="SM")
print("first 3 eigenvalues cost %s seconds" % (time.clock()-timeCheckin))
print("first 3 eigenvalues are %s" % vals)
k=len(vals)
nRow,nCol=A.get_shape()
for i in range(0,k):
print("error of lamda is %s " % (np.linalg.norm(A.dot(vecs[:,i])-vals[i]*vecs[:,i])))
for i in range(0,k):
v=np.random.rand(nCol)
#v must be normalization
v=v/np.linalg.norm(v)
print("random error is %s " % (np.linalg.norm(A.dot(v)-vals[i]*v)))

generateTestMat(10000)
calculateKEigs(tolerance=0.001)
#coding=utf8
import scipy as sp
import scipy.sparse.linalg
import scipy.io as sio
import numpy as np
import time
import os
def generateTestMat(size):
'''generate a test matrix for k-eigenvalues problem
'  matlab and scipy.sparse.linalg.eigs seem use the same package ARPACK
'  this matrix has two diagnol lines, one is on the diagnol and the other is off 1
'''
A = sp.sparse.lil_matrix((size, size))
A[0,:size]=np.random.rand(size)
A.setdiag(np.random.rand(size))
A.setdiag(np.random.rand(size-1),1)
#保存成matlab可以读取的格式,{"A":A}前面的A表示在matlab中的名字,后面的A表示在python中名字
sio.savemat("A.mat",{"A":A},oned_as='row')
print("generate a test matrix,with size %s by %s" %(size,size))

def calculateKEigs(tolerance=0):
if(os.path.exists(os.getcwd()+"/"+"A.mat")==False):
generateTestMat()
A=sio.loadmat("A.mat")["A"]#后面的A表示在matlab中名字
timeCheckin=time.clock()
#tol=0表示使用原先的计算精度
vals, vecs = sp.sparse.linalg.eigs(A, k=3,tol=tolerance,which="SM")
print("first 3 eigenvalues cost %s seconds" % (time.clock()-timeCheckin))
print("first 3 eigenvalues are %s" % vals)
k=len(vals)
nRow,nCol=A.get_shape()
for i in range(0,k):
print("error of lamda is %s " % (np.linalg.norm(A.dot(vecs[:,i])-vals[i]*vecs[:,i])))
for i in range(0,k):
v=np.random.rand(nCol)
#v must be normalization
v=v/np.linalg.norm(v)
print("random error is %s " % (np.linalg.norm(A.dot(v)-vals[i]*v)))

generateTestMat(10000)
calculateKEigs(tolerance=0.001)


结果如下,还是相当好的——10000*10000的矩阵耗时32.58秒,求出前3小的特征值,计算精度为0.001:

generate a test matrix,with size 10000 by 10000
first 3 eigenvalues cost 32.576715 seconds
first 3 eigenvalues are [-0.00547908+0.00513399j -0.00547908-0.00513399j -0.00547640+0.00724066j]
error of lamda is 3.44171050491e-06
error of lamda is 3.44171050491e-06
error of lamda is 7.36171280591e-06
random error is 43.5130937694
random error is 43.3207057797
random error is 43.3430705168

内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: