从基于列的数组返回多数加权投票

Question

我有一个矩阵x，具有3 x 3维数和一个3,的向量w：

x = np.array([[1, 2, 1],
[3, 2 ,1],
[1, 2, 2]])
w = np.array([0.3, 0.4, 0.3])

我需要生成另一个向量y，即每行x的多数票。x的每一列都按相应的值加权，以w为单位。像这样：

对于y[0]，它应该寻找X[0] => [1, 2, 1]

• 值为 1 的列 = 第一和第三 [0,2]
• 值为 2 = 秒的列 [1]
• 值为 3 的列 = 无

对按 X 中的值分组的列的权重(以w为单位)求和：

• 值为 1 的列的权重总和：0.3 + 0.3 = 0.6
• 值为 2 的列的权重总和：0.4
• 值为 3 的列的权重总和：0

由于值为 1 的列的权重总和最高，因此y[0] = 1.等等。

Answer 1

如果你懂广播，你可以用numpy来做。缺点是，由于代码是矢量化的，因此您执行的计算比您需要的要多。如果w向量的大小非常大，这将很重要。

也许有人想出了一种更简单的方法来写它，但这就是我不会考虑太多的方式。

答案先说：

i = np.arange(3) + 1
m = (x.reshape((1,4,3)) == i.reshape((3,1,1)))
np.argmax(np.sum(m, axis=2).T*w, axis=1) + 1

现在分步解释...请注意，通常最好从零开始计数，但我遵循了您的惯例。

1. 我添加了一行，因此数组不对称(更容易检查形状)

   In [1]: x = np.array([[1, 2, 1],
   ...:               [3, 2 ,1],
   ...:               [1, 2, 2],
   ...:               [3, 1, 3]])
   ...:
   ...: w = np.array([0.3, 0.4, 0.3])
   

2. 第一步是让索引数组i。您的约定从 1 开始。

   In [2]: i = np.arange(3) + 1
   

3. 棘手的步骤：创建一个形状为 (3,4,3) 的数组，其中数组的第 i 个条目是 (4,3) 数组，所有条目均为 0 或 1。它是 1 当且仅当 x == i。这是通过向x和i添加维度来完成的，以便可以广播它们。该操作基本上比较了x和i的所有组合，因为x的所有维度都匹配 size=1 维度的i，反之亦然：

   In [3]: m = (x.reshape((1,4,3)) == i.reshape((3,1,1)))*1
   In [4]: m
   Out[4]:
   array([[[1, 0, 1],
   [0, 0, 1],
   [1, 0, 0],
   [0, 1, 0]],
   [[0, 1, 0],
   [0, 1, 0],
   [0, 1, 1],
   [0, 0, 0]],
   [[0, 0, 0],
   [1, 0, 0],
   [0, 0, 0],
   [1, 0, 1]]])
   

4. 现在，您沿行求和(即 axis=2)，以获得每个选择在每行x中出现的次数(请注意，当您将其与x进行比较时，结果被转置：

   In [5]: np.sum(m, axis=2)
   Out[5]:
   array([[2, 1, 1, 1],
   [1, 1, 2, 0],
   [0, 1, 0, 2]])
   

5. 我希望你已经可以看到这是怎么回事。您可以直接阅读：在x的第一行中，1出现两次，2出现一次。在第二行x，都出现一次，在x的第三行，1出现一次，2出现两次，以此类推。

6. 将其乘以权重：

   In [7]: np.sum(m, axis=2).T*w
   Out[7]: 
   array([[0.6, 0.4, 0. ],
   [0.3, 0.4, 0.3],
   [0.3, 0.8, 0. ],
   [0.3, 0. , 0.6]])
   

7. 获取行的最大值(添加一个以符合您的约定)：

   In [8]: np.argmax(np.sum(m, axis=2).T*w, axis=1) + 1
   Out[8]: array([1, 2, 2, 3])
   

特殊情况：领带

评论中提到了以下案例：

x = np.array([[2, 2, 4, 1]])
w = np.array([0.1, 0.2, 0.3, 0.4])

权重的总和为：

[0.1, 0.4, 0., 0.4]

所以在这种情况下没有赢家。从这个问题中不清楚在这种情况下会怎么做。一个人可以拿走所有，什么都不拿...人们可以在最后查找这些情况：

final_w = np.sum(m, axis=2).T*w
result = np.argmax(np.sum(m*w, axis=2), axis=0) + 1
special_cases = np.argwhere(np.sum(final_w == np.max(final_w), axis=1) > 1)

注意：为了提高可读性，我使用了reshape方法，但我经常使用np.expand_dims或np.newaxis。像这样：

i = np.arange(3) + 1
m = (x[np.newaxis] == i[:, np.newaxis, np.newaxis])
np.argmax(np.sum(m, axis=2).T*w, axis=1) + 1

另一种选择：您也可以使用某种编译的代码。例如，在这种情况下，numba非常易于使用。

Answer 2

这是一个非常疯狂的方法，它涉及排序和索引，而不是添加新维度。这有点像np.unique使用的基于排序的方法。

首先找到每行中的排序索引：

rows = np.repeat(np.arange(x.shape[0]), x.shape[1])  # [0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3]
cols = np.argsort(x, axis=1).ravel()                 # [0, 2, 1, 2, 1, 0, 0, 1, 2, 1, 0, 2]

现在，您可以为每个列创建一个排序元素数组，包括未加权和加权。前者将用于获取求和的索引，后者实际上将被求和。

u = x[rows, cols]                            # [1, 1, 2, 1, 2, 3, 1, 2, 2, 1, 3, 3]
v = np.broadcast_to(w, x.shape)[rows, cols]  # [0.3, 0.3, 0.4, 0.3, 0.4, 0.3, 0.3, 0.4, 0.3, 0.4, 0.3, 0.3]

可以在以下位置找到要应用np.add.reduce的断点：

row_breaks = np.diff(rows).astype(bool)            # [0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0]
col_breaks = np.diff(u).astype(bool)               # [0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0]
break_mask = row_breaks | col_breaks               # [0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0]
breaks = np.r_[0, np.flatnonzero(break_mask) + 1]  # [ 0,  2,  3,  4,  5,  6,  7,  9, 10]

现在，您拥有每行中相同数字的权重总和：

sums = np.add.reduceat(v, breaks)  # [0.6, 0.4, 0.3, 0.4, 0.3, 0.3, 0.7, 0.4, 0.6]

但是您需要将它们分解为对应于每行唯一元素数量的段：

unique_counts = np.add.reduceat(break_mask, np.arange(0, x.size, x.shape[1]))
unique_counts[-1] += 1  # The last segment will be missing from the mask: # [2, 3, 2, 2]
unique_rows = np.repeat(np.arange(x.shape[0]), unique_counts)  # [0, 0, 1, 1, 1, 2, 2, 3, 3]

现在，您可以对每个段进行排序以查找最大值：

indices = np.lexsort(np.stack((sums, unique_rows), axis=0))  # [1, 0, 2, 4, 3, 5, 6, 7, 8]

每次运行结束时的索引由下式给出：

max_inds = np.cumsum(unique_counts) - 1  # [1, 4, 6, 8]

因此，最大金额为：

sums[indices[max_inds]]  # [0.6, 0.4, 0.7, 0.6]

您可以解开索引内索引，以从每一行中获取正确的元素。请注意，max_inds，以及依赖于它的所有内容都与x.shape[1]的大小相同，正如预期的那样：

result = u[breaks[indices[max_ind]]]

此方法看起来不是很漂亮，但它可能比在数组上使用额外的维度更节省空间。此外，无论x中的数字如何，它都可以工作。请注意，我从未以任何方式减去任何内容或调整x。事实上，所有行都是独立处理的，并且在构造breaks时，最大元素与下一个元素的最小元素相同的巧合被row_breaks打破。

TL;博士

享受：

def weighted_vote(x, w):
rows = np.repeat(np.arange(x.shape[0]), x.shape[1])
cols = np.argsort(x, axis=1).ravel()
u = x[rows, cols]
v = np.broadcast_to(w, x.shape)[rows, cols]
row_breaks = np.diff(rows).astype(bool)
col_breaks = np.diff(u).astype(bool)
break_mask = row_breaks | col_breaks
breaks = np.r_[0, np.flatnonzero(break_mask) + 1]
sums = np.add.reduceat(v, breaks)
unique_counts = np.add.reduceat(break_mask, np.arange(0, x.size, x.shape[1]))
unique_counts[-1] += 1
unique_rows = np.repeat(np.arange(x.shape[0]), unique_counts)
indices = np.lexsort(np.stack((sums, unique_rows), axis=0))
max_inds = np.cumsum(unique_counts) - 1
return u[breaks[indices[max_inds]]]

基准

基准测试按以下格式运行：

rows = ...
cols = ...
x = np.random.randint(cols, size=(rows, cols)) + 1
w = np.random.rand(cols)
%timeit weighted_vote_MP(x, w)
%timeit weighted_vote_JG(x, w)
assert (weighted_vote_MP(x, w) == weighted_vote_JG(x, w)).all()

我对weighted_vote_JG使用了以下概括，并进行了适当的更正：

def weighted_vote_JG(x, w):
i = np.arange(w.size) + 1
m = (x[None, ...] == i.reshape(-1, 1, 1))
return np.argmax(np.sum(m * w, axis=2), axis=0) + 1

行数： 100， 列数： 10

MP: 440 µs ± 5.12 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
* JG: 153 µs ± 796 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

行数： 1000， 列： 10

MP: 2.53 ms ± 43.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
* JG: 1.03 ms ± 12 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

行数： 10000， 列数： 10

MP: 23.5 ms ± 200 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
* JG: 16.6 ms ± 67.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

行数： 100000， 列： 10

MP: 322 ms ± 3.11 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
* JG: 188 ms ± 858 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

行数： 100， 列数： 100

* MP: 3.31 ms ± 257 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
JG: 12.6 ms ± 244 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

行数： 1000， 列数： 100

* MP: 31 ms ± 159 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
JG: 134 ms ± 581 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

行数： 10000， 列数： 100

* MP: 417 ms ± 7.06 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
JG: 1.42 s ± 126 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

行数： 100000， 列： 100

* MP: 4.94 s ± 25.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
JG: MemoryError: Unable to allocate 7.45 GiB for an array with shape (100, 100000, 100) and data type float64

故事的寓意：对于少量的列和权重，扩展的解决方案更快。对于更多的列，请改用我的版本。