为什么在一个大矩阵 $M$ 上的 NumPy 减法比将 $M$ 分成较小的矩阵然后减法时要慢？

Question

我正在编写一些代码，其中我有几个矩阵，并想从每个矩阵的每一行中减去一个向量 $v$（然后对结果做一些其他的事情）。当我使用 NumPy 并希望尽可能地“向量化”时，我认为我可以通过将所有矩阵存储为一个大（“连接”）矩阵并从中减去 $v$ 来加快运行时间。问题是我的代码在这个假定的优化之后运行速度变慢了。事实上，在某些情况下，分解矩阵并分别减去会明显更快（参见下面的代码示例）。

你能告诉我是什么原因造成的吗？天真地，我会假设这两种方法都需要相同数量的基本减法运算，而大矩阵方法更快，因为我们避免使用纯 Python 循环分别遍历所有矩阵。

最初，我认为减速可能是由于初始化了一个更大的矩阵来存储减法的结果。为了测试这一点，我在我的测试函数之外初始化了一个大矩阵并将它传递给 np.subtract 命令。然后我认为广播可能导致性能缓慢，所以我手动将向量广播成与大矩阵相同的形状，然后减去得到的广播矩阵。两种尝试都未能使大矩阵方法具有竞争力。

我制作了以下 MWE 来展示这个问题。

导入 NumPy 和一个计时器：

import numpy as np
from timeit import default_timer as timer

然后我有一些参数可以控制矩阵的大小和数量。

n = 100  # width of matrix
m = 500  # height of matrix
k = 100  # number of matrices
M = 100  # upper bound on entries
reps = 100  # repetitions for timings

我们可以生成如下的测试矩阵列表。大矩阵只是列表中所有矩阵的串联。我们从矩阵中减去的向量是随机生成的。

list_of_matrices = [np.random.randint(0, M+1, size=(m,n)) for _ in range(k)]
large_matrix = np.row_stack(list_of_matrices)
vector = np.random.randint(0, M+1, size=n)

这是我用来评估减法速度的三个函数。第一个从列表中的每个矩阵中减去向量，第二个从（连接的）大矩阵中减去向量，最后一个函数试图通过预初始化输出矩阵和广播向量来加速后一种方法。

def list_compute(list_of_matrices, vector):
    for j in range(k):
        np.subtract(list_of_matrices[j], vector)

def array_compute(bidlists, vector):
    np.subtract(large_matrix, vector_matrix, out=pre_allocated)

pre_allocated = np.empty(shape=large_matrix.shape)
vector_matrix = np.broadcast_to(vector, shape=large_matrix.shape)
def faster_array_compute(large_matrix, vector_matrix, out_matrix):
    np.subtract(large_matrix, vector_matrix, out=out_matrix)

我通过运行对这三个函数进行基准测试

start = timer()
for _ in range(reps):
    list_compute(list_of_matrices, vector)
print timer() - start

start = timer()
for _ in range(reps):
    array_compute(large_matrix, vector)
print timer() - start

start = timer()
for _ in range(reps):
    faster_array_compute(large_matrix, vector_matrix, pre_allocated)
print timer() - start

对于上述参数，我得到了时间

0.539432048798
1.12959504128
1.10976290703

天真地，与多个矩阵方法相比，我希望大矩阵方法更快或至少具有竞争力。我希望有人能给我一些见解，了解为什么不是这样，以及如何加快我的代码速度！

Answer 1

变量pre_allocated 的类型是float8。输入矩阵是 int。你有一个隐式转换。尝试将预分配修改为：

pre_allocated = np.empty_like(large_matrix)

更改之前，我的机器上的执行时间是：

0.6756095182868318
1.2262537249271794
1.250292605883855

改动后：

0.6776479894965846
0.6468182835551346
0.6538956945388001

所有情况下的性能都相似。这些测量值存在很大差异。甚至可以观察到第一个是最快的。

由于预分配，似乎没有任何收益。

请注意，分配非常快，因为它只保留地址空间。 RAM 仅在实际访问事件时消耗。缓冲区为 20MiB，因此 CPU 上的 L3 缓存更大。执行时间将主要由页面错误和缓存的重新填充决定。此外，对于第一种情况，内存在被释放后重新分配。对于内存分配器，该资源可能是“热的”。因此，您不能直接将解决方案 A 与其他解决方案进行比较。

修改第一种情况下的“action”行，保留实际结果：

        np.subtract(list_of_matrices[j], vector, out=pre_allocated[m*j:m*(j+1)])

然后向量化操作的收益变得更加明显：

0.8738251849091547
0.678185239557866
0.6830777283598941