切片稀疏(scipy)矩阵

Question

感谢任何帮助,以便在从scipy.sparse包中切割lil_matrix(A)时理解以下行为.

实际上,我想基于行和列的任意索引列表提取子矩阵.

当我使用这两行代码时:

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

一切都很好,我可以提取正确的子矩阵.

当我尝试在一行中执行此操作时,它失败了(返回的矩阵为空)

x=A[list 1,list 2]

为什么会这样？总的来说,我在matlab中使用了类似的命令,并在那里工作.那么,为什么不使用第一个,因为它有效？这似乎非常耗时.由于我必须经历大量的条目,我想使用单个命令加速它.也许我使用错误的稀疏矩阵类型......任何想法？

Answer 1

你已经使用的方法,

A[list1, :][:, list2]

似乎是从备件矩阵中选择所需值的最快方法.请参阅下面的基准.

但是,要回答有关如何A 使用单个索引从任意行和列中选择值的问题,您需要使用所谓的"高级索引":

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

使用高级索引,如果是arr1和arr2NDarrays,则等于(i,j)组件A[arr1, arr2]

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())

Answer 2

切片使用以下语法进行：

a[1:4]

对于 a = array([1,2,3,4,5,6,7,8,9])，结果是

array([2, 3, 4])

元组的第一个参数表示第一个要保留的值，第二个参数表示第一个不保留的值。

如果两侧都使用列表，则意味着数组的维度与列表长度一样多。

因此，根据您的语法，您可能需要这样的东西：

x = A[list1,:,list2]

取决于A的形状。

希望它确实对你有帮助。