切片稀疏（scipy）矩阵

Question

切片稀疏（scipy）矩阵

13

我需要帮助理解对scipy.sparse包中的lil_matrix(A)进行切片时出现的以下行为。

实际上，我想要基于任意行和列的索引列表提取一个子矩阵。

当我使用以下两行代码时：

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

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

当我尝试在一行中完成此操作时失败了（返回的矩阵为空）。

x=A[list 1,list 2]

为什么会这样？总体而言，我在matlab中使用了类似的命令，那里它可以工作。那么为什么不使用第一个呢，因为这样似乎非常耗时。由于我需要查看大量的数据，我希望能够使用单个命令加速。也许我使用了错误的稀疏矩阵类型...有什么想法吗？

- user972858

list1和list2的内容是什么？A[list1:list2]的含义是什么？ - Louis

list1和list2是包含整数的Python列表对象，例如[1,4,6,8]。A[list1:list2]是空的（<1x3稀疏矩阵，类型为'<type 'numpy.int32'>'，在链表格式中没有存储元素>）。 - user972858

4个回答

4

对我来说，unutbu提供的解决方案很好，但速度较慢。

我发现一个快速替代方法，

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

你可以看到行和列被单独切割，但是每个都转换成最快的稀疏格式，这次获取索引。

在我的测试环境中，此代码比另一个代码快1000倍。

我希望我没有说错或犯错。

- listener

2

同时索引，如B[arr1，arr2]是可行的，并且在我的机器上比listener's solution更快。请参见下面Jupyter示例中的In [5]。要将其与提到的答案进行比较，请参考In [6]。此外，我的解决方案不需要.tocsc()转换，这使得它更易读（在我看来）。

请注意，对于B[arr1，arr2]起作用，arr1和arr2必须是broadcastable numpy数组。

然而，一个更快的解决方案是使用B[list1][:，list2]，如pointed out by 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

- yogabonito

-1

切片使用这种语法进行操作：

a[1:4]

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

array([2, 3, 4])

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

如果两边都使用列表，则意味着您的数组具有与列表长度相同的维数。

因此，使用您的语法，您可能需要类似于以下内容的东西：

x = A[list1,:,list2]

根据A的形状而定。

希望对您有所帮助。

- Louis

1

问题不是关于array()的。 - Will

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- unutbu · Accepted Answer

你已经在使用的方法，

A[list1, :][:, list2]

似乎是从备用矩阵中选择所需值的最快方式。请参见下面的基准测试。

但是，如果要使用单个索引从任意行和列中选择A的值，您需要使用所谓的"高级索引":

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

使用高级索引，如果arr1和arr2是NDarrays，那么A[arr1, arr2]的(i,j)分量等于

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

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

通过使用广播（见下文），可以通过设置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://dev59.com/S2sz5IYBdhLWcg3w3bxb#26592783 (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())