高效的 numpy 子矩阵视图

Question

我希望将 Hungarian algorithm 应用于由列表 row_ind、col_ind 的叉积索引的 numpy 矩阵 C 的许多子集。目前，我看到以下选项：

双切片：

linear_sum_assignment(C[row_ind,:][:,col_ind])

问题： 每个子集操作两个副本。

高级切片通过 np.ix_:

linear_sum_assignment(C[np.ix_(row_ind, col_ind)])

问题：每个子集一个副本，np.ix_效率低下（分配n x n矩阵）。

更新：正如@hpaulj 所指出的，np.ix_ 它实际上不是分配 n x n 矩阵，但它仍然比 1 慢。

Masked array.

问题： 不适用于 linear_sum_assignment。

所以，没有一个选项是令人满意的。

理想情况下，我们希望能够使用矩阵 C 和几个分别用于行和列的一维掩码来指定子矩阵视图，因此可以将这样的视图传递给 linear_sum_assignment .对于另一个 linear_sum_assignment 调用，我会快速调整遮罩但从不修改或 copy/subset 完整矩阵。

numpy 中是否已经有类似的东西？

处理同一个大矩阵的多个子矩阵的最有效方法是什么（尽可能少 copies/memory 分配）？

Answer 1

索引数组的不同方式与lists/arrays时间差不多。他们都生产副本，而不是视图。

例如

In [99]: arr = np.ones((1000,1000),int)
In [100]: id1=np.arange(0,1000,10)
In [101]: id2=np.arange(0,1000,20)

In [105]: timeit arr[id1,:][:,id2].shape
52.5 µs ± 243 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [106]: timeit arr[np.ix_(id1,id2)].shape
66.5 µs ± 47.4 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

相比之下，如果我使用切片（在本例中选择相同的元素），我会得到 view，速度要快得多：

In [107]: timeit arr[::10,::20].shape
661 ns ± 18.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

ix_ 不会创建 (m,n) 数组；它 returns 调整后的一维数组的元组。相当于

In [108]: timeit arr[id1[:,None], id2].shape
54.5 µs ± 1.6 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

时间差异主要是由于多了一层函数调用。

你的 scipy link 有一个 [source] link:

https://github.com/scipy/scipy/blob/v0.19.1/scipy/optimize/_hungarian.py#L13-L107

此 optimize.linear_sum_assignment 函数使用 cost_matrix 创建一个 _Hungary 对象。这会制作一个副本，并通过搜索和操作其值来解决问题。

使用文档示例：

In [110]: optimize.linear_sum_assignment(cost)
Out[110]: (array([0, 1, 2], dtype=int32), array([1, 0, 2], dtype=int32))

它所做的是创建一个状态对象：

In [111]: H=optimize._hungarian._Hungary(cost)
In [112]: vars(H)
Out[112]: 
{'C': array([[4, 1, 3],
        [2, 0, 5],
        [3, 2, 2]]),
 'Z0_c': 0,
 'Z0_r': 0,
 'col_uncovered': array([ True,  True,  True], dtype=bool),
 'marked': array([[0, 0, 0],
        [0, 0, 0],
        [0, 0, 0]]),
 'path': array([[0, 0],
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0]]),
 'row_uncovered': array([ True,  True,  True], dtype=bool)}

迭代，

In [113]: step=optimize._hungarian._step1
In [114]: while step is not None:
     ...:     step = step(H)
     ...:

结果状态是：

In [115]: vars(H)
Out[115]: 
{'C': array([[1, 0, 1],
        [0, 0, 4],
        [0, 1, 0]]),
 'Z0_c': 0,
 'Z0_r': 1,
 'col_uncovered': array([False, False, False], dtype=bool),
 'marked': array([[0, 1, 0],
        [1, 0, 0],
        [0, 0, 1]]),
 'path': array([[1, 0],
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0]]),
 'row_uncovered': array([ True,  True,  True], dtype=bool)}

解决方案是从 marked 数组中提取的

In [116]: np.where(H.marked)
Out[116]: (array([0, 1, 2], dtype=int32), array([1, 0, 2], dtype=int32))

总成本是这些值的总和：

In [122]: cost[np.where(H.marked)]
Out[122]: array([1, 2, 2])

但是最终状态下C数组的代价是0:

In [124]: H.C[np.where(H.marked)]
Out[124]: array([0, 0, 0])

所以即使你给optimize.linear_sum_assignment的子矩阵是一个视图，搜索仍然涉及一个副本。搜索 space 和时间随着成本矩阵的大小而显着增加。

高效的 numpy 子矩阵视图

Efficient numpy submatrix view

numpy

scipy

hungarian-algorithm