切片稀疏（scipy）矩阵

Question

我希望能得到一些帮助，来理解在使用scipy.sparse包中的lil_matrix进行切片时出现的行为。

其实，我想根据一个任意的行和列索引列表来提取一个子矩阵。

当我使用这两行代码时：

x1 = A[list 1,:]
x2 = x1[:,list 2]

一切都很好，我能够提取到正确的子矩阵。

但是，当我尝试用一行代码来做这件事时，它失败了（返回的矩阵是空的）。

x=A[list 1,list 2]

这是为什么呢？总的来说，我在matlab中使用过类似的命令，并且它是可以工作的。那么，为什么不继续使用第一种方式呢？因为它有效。看起来这样做很耗时间。由于我需要处理大量的数据，我希望能通过一个命令来加快速度。也许我使用了错误的稀疏矩阵类型……有什么想法吗？

Answer 1

同时索引像 B[arr1, arr2] 这样的用法确实有效，并且在我的机器上比listener的解决方案要快。你可以看看下面Jupyter示例中的 In [5]。如果想和提到的答案做比较，可以参考 In [6]。而且，我的解决方案不需要 .tocsc() 转换，这让我觉得它更易读。

请注意，要让 B[arr1, arr2] 正常工作，arr1 和 arr2 必须是可广播的 numpy 数组。

不过，有一个更快的解决方案是使用 B[list1][:, list2]，正如unutbu指出的。你可以看看下面的 In [7]。

In [1]: from scipy import sparse
      : import numpy as np
      : 
      : 

In [2]: B = sparse.rand(1000, 1000, .1, format='lil')
      : list1=[1,4,6,8]
      : list2=[2,4]
      : 
      : 

In [3]: arr1 = np.array(list1)[:, None]  # make arr1 a (n x 1)-array
      : arr1
      : 
      : 
Out[3]: 
array([[1],
       [4],
       [6],
       [8]])

In [4]: arr2 = np.array(list2)[None, :]  # make arr2 a (1 x m)-array
      : arr2
      : 
      : 
Out[4]: array([[2, 4]])

In [5]: %timeit A = B.tocsr()[arr1, arr2]
100 loops, best of 3: 13.1 ms per loop

In [6]: %timeit A = B.tocsr()[np.array(list1),:].tocsc()[:,np.array(list2)]
100 loops, best of 3: 14.6 ms per loop

In [7]: %timeit B[list1][:, list2]
1000 loops, best of 3: 205 µs per loop

Answer 2

对我来说，unutbu 提供的解决方案效果不错，但速度比较慢。

我找到了一种快速的替代方法，

A = B.tocsr()[np.array(list1),:].tocsc()[:,np.array(list2)]

你可以看到行和列是分别处理的，但每个都转换成了最快的稀疏格式，这样这次可以更快地获取索引。

在我的测试环境中，这段代码比之前的快了1000倍。

希望我没有说错什么或者犯什么错误。

Answer 3

你现在使用的方法，

A[list1, :][:, list2]

似乎是从稀疏矩阵中选择所需值的最快方式。下面有一个基准测试的结果。

不过，关于如何用一个索引从A中选择任意行和列的值，你需要使用所谓的“高级索引”：

A[np.array(list1)[:,np.newaxis], np.array(list2)]

使用高级索引时，如果arr1和arr2是多维数组，那么A[arr1, arr2]中的(i,j)部分等于

A[arr1[i,j], arr2[i,j]]

因此，你希望arr1[i,j]等于list1[i]，对于所有的j，而arr2[i,j]等于list2[j]，对于所有的i。

这可以通过广播（见下文）的帮助来实现，设置为arr1 = np.array(list1)[:,np.newaxis]，和arr2 = np.array(list2)。

arr1的形状是(len(list1), 1)，而arr2的形状是(len(list2), )，这会广播成(1, len(list2))，因为在需要时会自动在左侧添加新的轴。

每个数组可以进一步广播到形状(len(list1),len(list2))。这正是我们想要的，因为我们希望A[arr1[i,j],arr2[i,j]]能够有意义，因为我们希望(i,j)遍历所有可能的索引，以得到形状为(len(list1),len(list2))的结果数组。

---

这里有一个微基准测试的结果，显示A[list1, :][:, list2]是最快的选项：

In [32]: %timeit orig(A, list1, list2)
10 loops, best of 3: 110 ms per loop

In [34]: %timeit using_listener(A, list1, list2)
1 loop, best of 3: 1.29 s per loop

In [33]: %timeit using_advanced_indexing(A, list1, list2)
1 loop, best of 3: 1.8 s per loop

---

这是我用来进行基准测试的设置：

import numpy as np
import scipy.sparse as sparse
import random
random.seed(1)

def setup(N):
    A = sparse.rand(N, N, .1, format='lil')
    list1 = np.random.choice(N, size=N//10, replace=False).tolist()
    list2 = np.random.choice(N, size=N//20, replace=False).tolist()
    return A, list1, list2

def orig(A, list1, list2):
    return A[list1, :][:, list2]

def using_advanced_indexing(A, list1, list2):
    B = A.tocsc()  # or `.tocsr()`
    B = B[np.array(list1)[:, np.newaxis], np.array(list2)]
    return B

def using_listener(A, list1, list2):
    """https://stackoverflow.com/a/26592783/190597 (listener)"""
    B = A.tocsr()[list1, :].tocsc()[:, list2]
    return B

N = 10000
A, list1, list2 = setup(N)
B = orig(A, list1, list2)
C = using_advanced_indexing(A, list1, list2)
D = using_listener(A, list1, list2)
assert np.allclose(B.toarray(), C.toarray())
assert np.allclose(B.toarray(), D.toarray())