python - 为什么 NumPy 在一个大矩阵 $M$ 上的减法比将 $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

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

看来预分配没有带来任何 yield 。

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

修改第一种情况中的“action”行以保留实际结果:

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

然后，矢量化操作的 yield 变得更加明显:

0.8738251849091547
0.678185239557866
0.6830777283598941