使用布尔矩阵替换 for 循环以执行高级索引答案

【问题标题】：Replace for loop by using a boolean matrix to perform advanced indexing使用布尔矩阵替换 for 循环以执行高级索引
【发布时间】：2018-02-20 06:21:33
【问题描述】：

当处理维度为 (A, B, C) 的 3 维矩阵“M”时，可以使用 2 个向量 X 索引 M，其中元素在 [0, A) 中，Y 中元素在 [0, B) 中同维D。

更具体地说，我在写作时明白这一点

M[X,Y,:]

对于 D 中的每个“i”，我们正在取

M[X[i], Y[i], :],

从而最终生成一个 DxC 矩阵。

现在假设

X is a numpy array of dim U, same concept as before
this time Y is a matrix UxL, where each row correspond to a Boolean numpy array 
(a mask)

看看下面的代码

for u in U:
    my_matrix[Y[u], X[u], :] += 1  # Y[u] is the mask that selects specific elements of the first dimension

我想编写没有 for 循环的相同代码。类似这样的东西

np.add.at(my_matrix, (Y, X), 1) # i use numpy.ufunc.at since same elements could occur multiple times in X or Y.

不幸的是返回以下错误

IndexError: 布尔索引与第 0 维的索引数组不匹配；维度是 L 但对应的布尔维度是 1

执行分配时也会发现此问题

for u in U:
    a_matrix[u, Y[u], :] = my_matrix[Y[u], X[u], :]

你知道我可以如何优雅地解决这个问题吗？

【问题讨论】：

您能否举一个最小的工作示例，可能是 3x3x2 矩阵或易于可视化的东西？
为了确保我理解，Y 是一个二维矩阵，Y[u] 是索引Y 维度的一维掩码。但是X 是一维矩阵，所以X[u] 是单个元素。换句话说，在工作 for 循环中，每次仅在所有三个维度中 Y 指示的行处将 1 添加到单个列。所以Y[u] 可能表示一些相同的行，这对应于将 1 多次添加到这些元素。对吗？
如果是这样，我认为最好的方法是先折叠Y 向量，然后计算一个求和矩阵以添加到整个事物中。
@AlexanderReynolds 正确。我已经用进一步的分配问题更新了这个问题，

标签： python arrays numpy matrix indexing

【解决方案1】：

简单地使用通常的 nd 数组形状的花式索引的直接方法不太可能解决您的问题。这就是我这么说的原因：Y 有布尔行，它告诉您在第一个维度上要采用哪些索引。所以Y[0] 和Y[1] 可能有不同数量的True 元素，因此Y 的行将沿第一维对具有不同长度的子数组进行切片。换句话说，您的数组形索引不能转换为矩形子数组。

但是，如果您考虑一下索引数组的含义，就会有出路。 Y 的行准确地告诉您要修改哪些元素。如果我们将所有索引杂乱无章地整理成一个庞大的一维花式索引集合，我们可以沿着我们想要索引的第一个维度精确定位每个 (x,y) 点。

请特别考虑以下示例（顺便说一下，您的问题严重缺失）：

A = np.arange(4*3*2).reshape(4,3,2)
Y = np.array([[True,False,False,True],
              [True,True,True,False],
              [True,False,False,True]])
X = np.array([2,1,2])

A 是形状 (4,3,2)，Y 是形状 (3,4)（第一行和最后一行故意相同），X 是形状 (3,)`（以及第一个和最后一个元素故意相同）。让我们将布尔索引转换为线性索引的集合：

U,inds = Y.nonzero()
#U: array([0, 0, 1, 1, 1, 2, 2])
#inds: array([0, 3, 0, 1, 2, 0, 3])

如您所见，U 是Y 中每个True 元素的行索引。这些是给出Y 的行和X 的元素之间对应关系的索引。第二个数组inds 是沿第一个维度的实际线性索引（对于给定的行）。

我们几乎完成了，我们只需将inds 的元素与来自X 的相应索引配对，用于第二维。这实际上很简单：我们只需要将X 与U 索引。

所以总而言之，以下两个是针对同一问题的等效循环和花式索引解决方案：

B = A.copy()
for u in range(X.size):
    A[Y[u],X[u],:] += 1
U,inds = Y.nonzero()
np.add.at(B,(inds,X[U]),1)

A 用循环修改，B 用np.add.at 修改。我们可以看到两者是相等的：

>>> (A == B).all()
True

如果你看一下这个例子，你会发现我故意复制了第一组和第三组索引。这表明np.add.at 与这些花哨的索引一起正常工作，并完成了在输入中多次出现的累积索引。（打印B并与A的初始值比较可以看到最后的项目增加了两次。）

【讨论】：