一个Scipy稀疏矩阵的索引操作向量化

Question

尽管所有东西似乎都已经向量化了，但以下代码运行速度仍然太慢。

from numpy import *
from scipy.sparse import *

n = 100000;
i = xrange(n); j = xrange(n);
data = ones(n);

A=csr_matrix((data,(i,j)));

x = A[i,j]

问题似乎是索引操作作为 Python 函数实现的，调用 A[i,j] 导致以下分析输出。

         500033 function calls in 8.718 CPU seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   100000    7.933    0.000    8.156    0.000 csr.py:265(_get_single_element)
        1    0.271    0.271    8.705    8.705 csr.py:177(__getitem__)
(...)

换句话说，Python函数_get_single_element被调用了100000次，这非常低效。为什么不使用纯C实现呢？有没有人知道如何解决这个限制，并加速上述代码？ 我应该使用不同的稀疏矩阵类型吗？

- D R

Answer 1

您可以使用A.diagonal()更快地检索对角线（在我的机器上为0.0009秒，而不是3.8秒）。但是，如果您想进行任意索引，则这是一个更复杂的问题，因为您使用的不仅仅是切片，还有索引列表。_get_single_element函数被调用了100,000次，因为它只是遍历您传递给它的迭代器（i和j）。切片将是A [30：60,10]或类似的内容。

此外，我建议使用csr_matrix(eye(n,n))来创建与迭代器相同的矩阵，以简化操作。

更新：

好的，既然您的问题确实是关于如何快速访问大量随机元素的，我会尽力回答您的问题。

为什么没有用纯C实现？

答案很简单：还没有人去做。从我看到的情况来看，Scipy的稀疏矩阵模块仍有大量工作要做。其中一个部分是在不同稀疏矩阵格式之间进行转换。

是否有人知道解决此限制并加快上述代码的方法？

您可以尝试深入研究稀疏矩阵模块，并尝试加速它们。我这样做时，使用csr矩阵尝试了您上面的代码，以随机访问方式获得的时间少于原始时间的三分之一。我不得不直接访问_get_single_element并大大简化代码，包括删除边界检查。

然而，使用列表推导式进行csr_matrix的访问仍然比lil_matrix方法要慢得多。最终，lil矩阵在访问随机元素方面更快，因为它没有压缩。

在其原始形式下使用lil_matrix的运行时间约为上面列出的代码的五分之一。如果我去掉一些边界检查并直接调用lil_matrix的_get1()方法，则可以将时间进一步降低到原始时间的7％。为了清晰起见，这是从3.4-3.8秒加速到约0.261秒。

最后，我尝试编写自己的函数，直接访问lil矩阵的数据并避免重复的函数调用。这需要大约0.136秒的时间。这没有利用数据的排序，这是另一个潜在的优化（特别是如果您访问的许多元素在同一行）。

如果您想要更快的速度，那么您可能需要编写自己的C代码稀疏矩阵实现。

我应该使用不同的稀疏矩阵类型吗？

如果你的目的是访问大量元素，我建议使用lil matrix，但这取决于你需要做什么。例如，你是否需要乘以矩阵？请记住，在某些情况下，切换矩阵格式可能相当快速，因此不要排除改变矩阵格式以执行不同的操作。

如果你不需要对矩阵执行任何代数运算，那么可能应该使用defaultdict或类似的东西。defaultdict的危险在于，每当请求不存在于字典中的元素时，它会将该项设置为默认值并存储它，这可能会有问题。

- Justin Peel

---

我只是为了演示目的使用了一个对角矩阵。在我的情况下，指数i，j不沿对角线，并且矩阵不一定是对角线矩阵。 - D R

Answer 2

我认为只有在使用默认的'dtype'为'object'时才会调用_get_single_element。你是否尝试过提供'dtype'，例如csr_matrix((data, (i,j)), dtype=int32)

- Shane Holloway

---

只是让你知道，我尝试了一下，但似乎没有任何区别。不过还是值得一试的。 - Justin Peel