Numpy中的置换矩阵乘法

如何解决Numpy中的置换矩阵乘法

我有一个条目乱序的矩阵。例如，A = [a1,a2] 是具有分区 a1 和 a2 的典型矩阵（a1 和 a2 不一定是标量，但可以是子数组）。我的矩阵的顺序是 A = [a2,a1]。

我仍然想执行正常的矩阵乘法，就好像 A 的顺序正确一样。有没有办法在不复制的情况下做到这一点？

此用例是在循环缓冲区上进行乘法运算。我一直在尝试避免使用 numpy.roll 作为缓冲区。很多时候，循环缓冲区是碎片化的，但我只想相乘，因为我知道缓冲区碎片在哪里。

似乎无论我做什么，我最终都必须分配与最终矩阵乘积大小相等的元素，因为乘以 A 的任何分区仍将具有相同的列数，所以最终维度不会改变，因为 A1（A 中将乘以分区 1 的部分）将是 m x k1 并且分区将是 k1 x n（其中 k1 是分区中的行数）。

np.matmul 似乎没有添加到现有数组而不是覆盖数组的选项。这让我很困扰，因为不应该需要另一个分配，因为元素只是按不同的顺序排列。我想我可以相应地映射元素，但事实证明比这更困难。

示例：

data = np.empty((10,3))
buffer = CircularBuffer(data)

for i in range(10):
   data.push(i)

# buffer now reads [[0,0],[1,1,1],...]

data.push(-1)
data.push(-2)
# buffer overwrites oldest element
# reads [[-1,-1,-1],[-2,-2,-2],[2,2,2],...] 

A = np.arange(30).reshape(3,10)

# We want this to perform A @ [[2,...,[-1,-2]]
#                              buf_start                      buf_end
# Unfortunately,the buffer is out of order!
A @ buffer

显然缓冲区中有两个“分区”。从缓冲区的开始到存储元素的结尾 [[2,...] 以及从存储元素的开始到缓冲区的结尾 [[-1,-2]]。

我们知道它们应该采用的顺序（显示在代码中），但我们希望乘法按照 BUFFER 的顺序而不是原始数据本身的顺序表示，我想避免复制。

解决方法

发布的问题提到 A = [a1,a2] 其中 a1 和 a2 可能是子数组而不是标量。

因此，对于我的工作示例，我将 A.shape 定义为 (10,3,2)。 根据矩阵乘法的 numpy 规则（适用于 np.matmul() 以及 @ 运算符以及 np.dot()），对于 A @ buffer 是有效的， A 的最后一个维度必须与 buffer 的倒数第二个维度匹配。因此，我将 buffer.shape 定义为 (10,2,3)。 （注意与 A 的形状不同）

这是我对 A 和 buffer 的初始化：

A      = np.arange(60).reshape (10,3)
buffer = -1 * np.arange(60).reshape (10,2)

然后，对于这些形状，A @ buffer 的形状为 (10,2)。因此，为了保存矩阵乘法的结果，我将数组 my_result 定义为：

my_result = np.empty(A.shape[:-1]+buffer.shape[-1:]) # To hold the result of A @ buffer

现在，问题中提到buffer的行是乱序的。

让 k 成为“缓冲区碎片/分区”的索引（在 axis0 的 buffer 上）。 （因此，0 <= k <= 9）。

k = 2

那么，A 和 buffer 之间的矩阵乘法不是一次性执行，而是分两步完成，如下所示：

np.matmul(A[:k],buffer[-k:],out=my_result[:k])     # Fills a part of my_result
np.matmul(A[k:],buffer[:-k],out=my_result[k:])     # Fills the other part of my_result

测试：

my_result1 = np.matmul(A,np.roll(buffer,k,axis=0))
np.array_equal(my_result,my_result1)

输出：

正确

总而言之，我们在这里做的是：

• 分配一个未初始化的数组 my_result 来保存 matmul() 的结果
• 调用 matmul() 两次，每次，传递不同的 my_result 切片作为 out 的 matmul() 参数的值，并注意索引 {{1 }} 和 A 以正确的方式。

据我所知，这里不涉及冗余数据复制，除非您将 buffer 视为某种数据复制。这只是发生在两部分中的原始矩阵乘法。