scipy.sparse求稀疏矩阵前k个特征值
背景:
要在python中处理7000*7000的稀疏矩阵,计算前k小的特征值和相应的特征向量。不想在matlab中做这件事了,所有的数据预处理和展现工作都想在python中完成。然而一般的linalg提供的eig开销太大,要计算所有的特征值和特征向量,这个开销要达到 O(N^3),对于谱聚类来说,这个开销是不能忍受的。
所以要借助稀疏矩阵计算的工具包。
探索过程:
使用scipy.sparse
对于10*10大小的矩阵都没有问题:
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的矩阵就出错了:
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的参数也应该是稀疏矩阵,于是修改了代码:
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,初步判断是和算法的收敛状况有关。
另外使用MATLAB做对比,matlab代码是
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
最终解决方案:
#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 normalizationv=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
scipy.sparse求稀疏矩阵前k个特征值相关推荐
- 关于稀疏矩阵转化为稠密矩阵问题 (scipy.sparse格式和tensor稀疏张量格式)
scipy.sparse: todense() pytorch中的稀疏张量tensor: to_dense()
- scipy.sparse、pandas.sparse、sklearn稀疏矩阵的使用
单机环境下,如果特征较为稀疏且矩阵较大,那么就会出现内存问题,如果不上分布式 + 不用Mars/Dask/CuPy等工具,那么稀疏矩阵就是一条比较容易实现的路. 文章目录 1 scipy.sparse ...
- scipy.sparse使用简例
CDIMC-Net[1] 中有个对整个数据集求 kNN 图的函数 get_kNNgraph2[2],是用 dense 的 numpy.ndarray 存的,空间复杂度 O ( n 2 ) O(n^2) ...
- 【scipy.sparse中csr.matrix的用法】
scipy.sparse中csr.matrix的用法 作用:用于压缩稀疏行矩阵 1.csr_matrix(D) with a dense matrix or rank-2 ndarray D 2.cs ...
- matlab cell2mat 一维,MATLAB cell2mat(...)與具有一束稀疏矩陣的單元陣列功能,內存溢出意外...
我想什麼做的是: cell_array_outer = cell(1,N) parfor k = 1:N cell_array_inner = cell(1,M); for i = 1:M A = d ...
- scipy.sparse.coo_matrix、csr_matrix、lil_matrix、dia_matrix
文章目录 coo_matrix csr_matrix lil_matrix dia_matrix coo_matrix 1.coo啥意思?COOrdinate(坐标) 2.那么coo_matrix又是 ...
- matlab计算原点矩,关于用matlab求样本均值方差以及k阶原点矩的matlab程序
关于用matlab求样本均值方差以及k阶原点矩的matlab 程序 关于用matlab求样本均值和方差以及matlab程 序 1n1. 样本均值,公式xX,(其中X为样本).程序如下: ,i,1in ...
- scipy.sparse.csr_matrix函数和coo_matrix函数
Scipy高级科学计算库:和Numpy联系很密切,Scipy一般都是操控Numpy数组来进行科学计算.统计分析,所以可以说是基于Numpy之上了.Scipy有很多子模块可以应对不同的应用,例如插值运算 ...
- python 的csr_Python scipy.sparse.csr_matrix()[csc_matrix()]
csr_matrix是Compressed Sparse Row matrix的缩写组合,下面介绍其两种初始化方法 csr_matrix((data, (row_ind, col_ind)), [sh ...
最新文章
- vue 如何判断两个数组相同_如何判断车头与障碍物的距离,教你两个办法,轻松靠墙10公分...
- Openstack Paste.ini 文件详解
- Android开发之Mediaplayer
- STM32应用实例六:与MS5837压力传感器的I2C通讯
- ASP.NET发送电子邮件
- 华为重磅反击,鸿蒙来了!
- 金融项目app服务器配置,云在金融的应用
- 【Elasticsearch】es 7 Failed to parse value [analyzed] as only [true] or [false] are allowed
- PHP 实现定时任务的几种方法
- loadrunner录制脚本参数化之间的关联设置
- html让时间只展示年月日_JS 如何动态显示当前年月日时分秒-百度经验
- C#毕业设计——基于C#+asp.net+cs的即时通信系统设计与实现(毕业论文+程序源码)——即时通信系统
- 数据安全“考题”怎么破解?11月2日厦门站算力私享会开启
- 服务器查看光模块信息的命令,通过命令行界面(CLI)查看在交换机的光模块状态...
- Sublime Text 2 - 性感无比的代码编辑器!程序员必备神器!跨平台支持Win/Mac/Linux,支持32与64位,支持各种流行编程语言的语法高亮、代码补全等...
- c语言判断字符是否对称,2020-07-23(C语言)数据结构-试设计算法判断该链表的全部n个字符是否中心对称。...
- 支付宝 收款通知 mysql_基于支付宝微信通知的一种个人收款回调方案(转)
- 垃圾填埋场渗滤液厌氧处理过程中沼气的综合利用
- 2019南京帆软春招
- 腾讯开发平台php,腾讯AI开放平台 Tencent AI open platform
热门文章
- java图像处理-(指定区域内)灰度化、透明化(alpha通道)处理
- C(判断一个字母是否为英语字母)
- 四个程序员编辑器,编程必备!!!
- python数据分析与可视化从入门到精通_零基础学Python爬虫、数据分析与可视化从入门到精通...
- 全网最全ebay大数据面经合集
- CSS图片保持原比例
- zxing扫描条形码 ios
- 零基础学Python———求一个字符串的每个字符重新组合排列python排列组合的数学运算(递归法)
- 华为nova5里面有用到鸿蒙吗,华为nova7se是不是鸿蒙系统?
- 类似ftp文件服务器有哪些,FTP的替代品有哪些,你知道吗?