草庐IT

python - Scipy 稀疏矩阵 - 密集向量乘法性能 - block 与大型矩阵

coder 2023-08-24 原文

我有许多 scipy 稀疏矩阵(目前为 CSR 格式),我需要将它们与密集的 numpy 一维向量相乘。 该向量称为 G:

print G.shape, G.dtype
(2097152,) complex64

每个稀疏矩阵的形状都是 (16384,2097152) 并且非常稀疏。密度约为4.0e-6。 我有一个包含 100 个稀疏矩阵的列表,称为 spmats

我可以像这样轻松地将每个矩阵与 G 相乘:

res = [spmat.dot(G) for spmat in spmats]

这会按预期生成形状为 (16384,) 的密集向量列表。

我的应用程序对性能相当关键,所以我尝试了一种替代方法,即首先将所有稀疏矩阵连接成一个大的稀疏矩阵,然后只使用一次 dot() 调用像这样:

import scipy.sparse as sp
SPMAT = sp.vstack(spmats, format='csr')
RES = SPMAT.dot(G)

这会产生一个长向量 RES,其形状为 (1638400,),并且是 res 中所有结果向量的串联版本> 以上,正如预期的那样。我检查过结果是否相同。

也许我完全错了,但我预计第二种情况应该比第一种情况更快,因为 numpy 调用、内存分配、python 对象的创建、python 循环等要少得多。我不在乎关于连接稀疏矩阵所需的时间,只有计算结果的时间。然而,根据 %timeit:

%timeit res = [spmat.dot(G) for spmat in spmats]
10 loops, best of 3: 91.5 ms per loop
%timeit RES = SPMAT.dot(G)
1 loops, best of 3: 389 ms per loop

我已经检查过我在这两个操作中都没有用完内存,而且似乎没有任何可疑的事情发生。我是疯了,还是这真的很奇怪?这是否意味着所有稀疏矩阵向量乘积都应该分 block 进行,一次几行,以使其更快? 据我了解,稀疏矩阵与密集向量的乘法时间应该与非零元素的数量成线性关系,这在上述两种情况下都没有变化。是什么造成了如此大的不同?

我在一个单核 linux 机器上运行,4GB 内存,使用 EPD7.3

编辑:

这是一个为我重现问题的小例子:

import scipy.sparse as sp
import numpy as n

G = n.random.rand(128**3) + 1.0j*n.random.rand(128**3)

spmats = [sp.rand (128**2, 128**3, density = 4e-6, format = 'csr', dtype=float64) for i in range(100)]
SPMAT = sp.vstack(spmats, format='csr')

%timeit res = [spmat.dot(G) for spmat in spmats]
%timeit RES = SPMAT.dot(G)

我得到:

1 loops, best of 3: 704 ms per loop
1 loops, best of 3: 1.34 s per loop

这种情况下的性能差异不如我自己的具有某种结构的稀疏矩阵大(可能是因为缓存),但连接矩阵仍然更糟。

我已经尝试过 scipy 10.1 和 12.0。

最佳答案

我还没有找到问题中提到的奇怪行为的原因,但是我找到了一种可以显着加快计算速度的方法,这可能对其他人有用。

因为在我的特定情况下,我正在计算 float32 稀疏矩阵和 complex64 密集向量的乘积,所以我可以分别乘以实部和虚部。这为我提供了 4 倍的加速。

使用 SPMAT.shape == (16384000, 2097152) 需要 2.35 秒:

RES = SPMAT.dot(G)

虽然这只需要 541 毫秒:

RES = n.zeros((SPMAT.shape[0],),dtype=complex64)
RES.real = SPMAT.dot(G.real); RES.imag = SPMAT.dot(G.imag)

结果是一样的。 我想也许 n.zeros 预分配可能不是必需的,但我不知道该怎么做。

关于python - Scipy 稀疏矩阵 - 密集向量乘法性能 - block 与大型矩阵,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/16746974/

有关python - Scipy 稀疏矩阵 - 密集向量乘法性能 - block 与大型矩阵的更多相关文章

  1. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声? - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题?更新问题,以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点?我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等),但我确实想创建一个输出文件。

  2. ruby - RSpec - 使用测试替身作为 block 参数 - 2

    我有一些Ruby代码,如下所示:Something.createdo|x|x.foo=barend我想编写一个测试,它使用double代替block参数x,这样我就可以调用:x_double.should_receive(:foo).with("whatever").这可能吗? 最佳答案 specify'something'dox=doublex.should_receive(:foo=).with("whatever")Something.should_receive(:create).and_yield(x)#callthere

  3. ruby-on-rails - Enumerator.new 如何处理已通过的 block ? - 2

    我在理解Enumerator.new方法的工作原理时遇到了一些困难。假设文档中的示例:fib=Enumerator.newdo|y|a=b=1loopdoy[1,1,2,3,5,8,13,21,34,55]循环中断条件在哪里,它如何知道循环应该迭代多少次(因为它没有任何明确的中断条件并且看起来像无限循环)? 最佳答案 Enumerator使用Fibers在内部。您的示例等效于:require'fiber'fiber=Fiber.newdoa=b=1loopdoFiber.yieldaa,b=b,a+bendend10.times.m

  4. ruby - 在匿名 block 中产生 - 2

    我没有理解以下行为(另请参阅inthisSOthread):defdef_testputs'def_test.in'yieldifblock_given?puts'def_test.out'enddef_testdoputs'def_testok'endblock_test=procdo|&block|puts'block_test.in'block.callifblockputs'block_test.out'endblock_test.calldoputs'block_test'endproc_test=procdoputs'proc_test.in'yieldifblock_gi

  5. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意,但是谷歌搜索||=并不是很有帮助;)Python中是否有与Ruby和Perl中的||=语句等效的语句?例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外,类似这样的东西的通用术语是什么?条件分配是我的第一个猜测,但Wikipediapage跟我想的不太一样。

  6. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi? - 2

    什么是ruby​​的rack或python的Java的wsgi?还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

  7. 旋转矩阵的几何意义 - 2

    点向量坐标矩阵的几何意义介绍旋转矩阵的几何含义之前,先介绍一下点向量坐标矩阵的几何含义点:在一维空间下就是一个标量,如同一条直线上,以任意某一个位置为0点,以一定的尺度间隔为1,2,3...,相反方向为-1,-2,-3...;如此就形成了一维坐标系,这时候任何一个点都可以用一个数值表示,如点p1=5,即即从原点出发沿着x轴正方向移动5个尺度;点p2=-3,负方向移动3个尺度;     在一维坐标系上过原点做垂直于一维坐标系的直线,则形成了二维坐标系,此时描述一个点需要两个数值来表示点p3=(3,2),即从原点出发沿着x轴正方向移动3个尺度,在此基础上沿着y轴正方向移动两个尺度的位置就是点p3。

  8. 华为OD机试用Python实现 -【明明的随机数】 2023Q1A - 2

    华为OD机试题本篇题目:明明的随机数题目输入描述输出描述:示例1输入输出说明代码编写思路最近更新的博客华为od2023|什么是华为od,od薪资待遇,od机试题清单华为OD机试真题大全,用Python解华为机试题|机试宝典【华为OD机试】全流程解析+经验分享,题型分享,防作弊指南华为o

  9. python - 如何读取 MIDI 文件、更改其乐器并将其写回? - 2

    我想解析一个已经存在的.mid文件,改变它的乐器,例如从“acousticgrandpiano”到“violin”,然后将它保存回去或作为另一个.mid文件。根据我在文档中看到的内容,该乐器通过program_change或patch_change指令进行了更改,但我找不到任何在已经存在的MIDI文件中执行此操作的库.他们似乎都只支持从头开始创建的MIDI文件。 最佳答案 MIDIpackage会为您完成此操作,但具体方法取决于midi文件的原始内容。一个MIDI文件由一个或多个音轨组成,每个音轨是十六个channel中任何一个上的

  10. 「Python|Selenium|场景案例」如何定位iframe中的元素? - 2

    本文主要介绍在使用Selenium进行自动化测试或者任务时,对于使用了iframe的页面,如何定位iframe中的元素文章目录场景描述解决方案具体代码场景描述当我们在使用Selenium进行自动化测试的时候,可能会遇到一些界面或者窗体是使用HTML的iframe标签进行承载的。对于iframe中的标签,如果直接查找是无法找到的,会抛出没有找到元素的异常。比如近在咫尺的例子就是,CSDN的登录窗体就是使用的iframe,大家可以尝试通过F12开发者模式查看到的tag_name,class_name,id或者xpath来定位中的页面元素,会抛出NoSuchElementException异常。解决

随机推荐