我将不胜感激任何帮助,以了解从 scipy.sparse 包中切片 lil_matrix (A) 时的以下行为。
实际上,我想根据行和列的任意索引列表提取子矩阵。
当我使用这两行代码时:
x1 = A[list 1,:]
x2 = x1[:,list 2]
一切都很好,我可以提取正确的子矩阵。
当我试图在一行中执行此操作时,它失败了(返回矩阵为空)
x=A[list 1,list 2]
为什么会这样?总的来说,我在 matlab 中使用了类似的命令并且它可以工作。 那么,为什么不使用第一个,因为它有效?这似乎很费时间。由于我必须浏览大量条目,因此我想使用单个命令加快速度。也许我使用了错误的稀疏矩阵类型...有什么想法吗?
最佳答案
你已经在使用的方法,
A[list1, :][:, list2]
似乎是从备用矩阵中选择所需值的最快方法。请参阅下面的基准。
但是,要回答有关如何使用单个索引从 A
的任意行和列中选择值的问题,
你需要使用所谓的 "advanced indexing" :
A[np.array(list1)[:,np.newaxis], np.array(list2)]
使用高级索引,如果 arr1
和 arr2
是 NDarray,A[arr1, arr2]
等于
A[arr1[i,j], arr2[i,j]]
因此,对于所有 j
,您希望 arr1[i,j]
等于 list1[i]
,并且
arr2[i,j]
等于 list2[j]
对于所有 i
。
这可以在 broadcasting 的帮助下安排(见下文)通过设置
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)
遍历 a 的所有可能索引形状 (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())
关于python - 切片稀疏(scipy)矩阵,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/7609108/