不规则广播视图 numpy答案

【问题标题】：broadcast views irregularly numpy不规则广播视图 numpy
【发布时间】：2018-06-01 07:47:16
【问题描述】：

假设我想要一个大小为(n,m) 的numpy 数组，其中n 非常大，但有很多重复，即。 0:n1 相同，n1:n2 相同等（与n2%n1!=0，即不是定期间隔）。有没有办法在查看整个数组的同时为每个重复项只存储一组值？

示例：

unique_values = np.array([[1,1,1], [2,2,2] ,[3,3,3]]) #these are the values i want to store in memory
index_mapping = np.array([0,0,1,1,1,2,2]) # a mapping between index of array above, with array below

unique_values_view = np.array([[1,1,1],[1,1,1],[2,2,2],[2,2,2],[2,2,2], [3,3,3],[3,3,3]]) #this is how I want the view to look like for broadcasting reasons

我打算将数组（视图）乘以其他大小为(1,m) 的数组，然后取这个乘积的点积：

other_array1 = np.arange(unique_values.shape[1]).reshape(1,-1) # (1,m)
other_array2 = 2*np.ones((unique_values.shape[1],1)) # (m,1)
output = np.dot(unique_values_view * other_array1, other_array2).squeeze()

输出是一个长度为n的一维数组。

【问题讨论】：

您打算如何使用输出？仅供参考：第一阶段的视图，不保证以后不会强制复制。
@Divakar 真的。我计划与另一个形状为(1,m) 的数组相乘，然后将点积与另一个数组一起存储。主要考虑因素是将数组安装到内存中。如果强制复制，我可以分块完成最后一步
或者可以添加一些更通用的最小样本，比如index_mapping 具有更大的数字范围，unique_values 具有随机数？
“相同”是指它们是整数吗？或者只是非常接近的浮动？（即你能做按位等价吗）
@M.T Yakym Pirozhenko 的解决方案似乎适用于列出的最小案例。但同样，我看不出通用案例是什么。

标签： python numpy array-broadcasting

【解决方案1】：

根据您的示例，您可以简单地将索引映射分解到最后：

output2 = np.dot(unique_values * other_array1, other_array2).squeeze()[index_mapping]

assert (output == output2).all()

【讨论】：

【解决方案2】：

您的表达式有两个重要的优化：

最后做索引
先将other_array1 与other_array2 相乘，然后再与 unique_values

让我们应用这些优化：

>>> output_pp = (unique_values @ (other_array1.ravel() * other_array2.ravel()))[index_mapping]

# check for correctness
>>> (output == output_pp).all()
True

# and compare it to @Yakym Pirozhenko's approach
>>> from timeit import timeit
>>> print("yp:", timeit("np.dot(unique_values * other_array1, other_array2).squeeze()[index_mapping]", globals=globals()))
yp: 3.9105667411349714
>>> print("pp:", timeit("(unique_values @ (other_array1.ravel() * other_array2.ravel()))[index_mapping]", globals=globals()))
pp: 2.2684884609188884

如果我们观察两件事，这些优化很容易发现：

(1) 如果A 是mxn-矩阵且b 是n-向量，则

A * b == A @ diag(b)
A.T * b[:, None] == diag(b) @ A.T

(2) 如果A 是一个mxn-矩阵并且I 是一个k- 整数向量 range(m)然后

A[I] == onehot(I) @ A

onehot可以定义为

def onehot(I, m, dtype=int):
    out = np.zeros((I.size, m), dtype=dtype)
    out[np.arange(I.size), I] = 1
    return out

使用这些事实并缩写uv、im、oa1和oa2我们可以写

uv[im] * oa1 @ oa2 == onehot(im) @ uv @ diag(oa1) @ oa2

上述优化现在只是为这些矩阵乘法选择最佳顺序的问题，即

onehot(im) @ (uv @ (diag(oa1) @ oa2))

在此基础上逆向使用 (1) 和 (2)，我们得到了本文开头的优化表达式。

【讨论】：