主页 / coding / 问题 / 79248192
Accepted
julian2000P
Asked: 2024-12-03 23:21:47 +0800 CST

为什么我使用 numpy 的矩阵乘法这么慢？

  • 772

我正在尝试在 numpy 中将两个维数相当大的矩阵相乘。请参阅下面的 3 种方法。我随机实现了这 3 个矩阵来展示我的问题。第一个矩阵，即Y1[:,:,0]首先是一个更大的 3d 数组的一部分。第二个是.copy()这个矩阵的，第三个是它自己的矩阵。

为什么第一次乘法比后两次乘法慢这么多？

import numpy as np
from time import time

Y1 = np.random.uniform(-1, 1, (5000, 1093, 201))
Y2 = Y1[:,:,0].copy()
Y3 = np.random.uniform(-1, 1, (5000, 1093))

W = np.random.uniform(-1, 1, (1093, 30))

# method 1
START = time()
Y1[:,:,0].dot(W)
END = time()
print(f"Method 1 : {END - START}")

# method 2
START = time()
Y2.dot(W)
END = time()
print(f"Method 2 : {END - START}")

# method 3
START = time()
Y3.dot(W)
END = time()
print(f"Method 3 : {END - START}")

输出时间分别大约为34、0.06、0.06秒。

我看到了区别：虽然最后两个矩阵是“真正的”二维数组，但第一个矩阵是我更大的三维数组的一部分。

子集化Y1[:,:,0]导致速度如此缓慢吗？另外，我注意到为矩阵 Y2 创建 Y1 的副本也相当慢。

毕竟，我得到了这个 3d 数组，并且必须重复计算 Y1 切片与（可能不同的）矩阵 W 的矩阵乘积。有没有更好/更快的方法来做到这一点？

提前致谢！

numpy

3 个回答

  1. Best Answer

    这是一个缓存问题。如果你研究成本差异与第三轴的大小，你首先会发现一个线性关系（k=1 => 没有差异，k=2，方法 1 的成本是两倍，k=3，方法 1 的成本是三倍，等等），上限是最大值（k=20 或 k=30 不会真正改变情况）

    该最大上限取决于其他轴的大小

    问题是，矩阵乘法（以及基本上任何数组运算）经常迭代进行。因此，内存中的数据是依次读取的。

    第一次读取数据会花费一些时间，因为内存很慢。但是当您访问内存中的数据时，会读取整行（例如 64 或 128 字节）并将其存储到缓存中。如果下一个操作使用矩阵中的下一个数字，并且该数字恰好位于内存中前一个数字的旁边，则它很可能属于同一缓存行。并且没有必要在内存中读取它，我们将它放在（更快的）缓存中。

    这有点过于简单了。而且它如何应用于矩阵乘法并不那么明显，因为矩阵乘法不是那么连续的。但基本上，在内存中使用彼此接近的数据越多，速度就越快。人们经常忽略这一点，认为这是一种黑客优化，以节省一些额外的纳秒。但效果可能很大。

    对于非常少量的数据，完全适合缓存（几千字节），以及足够复杂的算法，可以多次读取它（即使只是矩阵乘法也符合条件），它并没有真正显示出来，因为每个数据都会在几个计算步骤后进入缓存。

    但是如果缓存中装不下所有内容，那么数据之间的间隔越大，您重复使用缓存行的可能性就越小。而且您越需要重新读取内存中的数据。以至于读取内存是最大的成本。

    所以你的问题是

    Y1=np.random.uniform(-1, 1, (5000, 1093, 201))

    的每个数据Y1[:,:,0]至少相隔 201x8 = 1608 个字节。因此，缓存（对于Y1—— 它仍然用于W，但所有方法对此都相同）是无用的：没有机会快速访问 的值，Y1[:,:,0]因为我们已经在内存中读取了一个接近它的值：它们彼此相距很远。

    另一种方法可以让你相信这是你的问题，如果需要的话，这可能是解决方案，只要看看如果

    Y1=np.random.uniform(-1, 1, (201, 5000, 1093)).transpose(1,2,0)

    Y1和你的形状完全相同。同样的 5000x1093x201 数组。并且你保留了Y1[:,:,0]形状为 5000x1093 的相同子数据。

    Y1从纯粹的“数学”角度来看，我的和你的之间的唯一区别是看不见的；唯一的区别在于数据在物理内存中的确切存储位置。

    在我的 Y1 中，Y1[i,j,0]远离Y1[i+1,j,0]，并且远离Y1[i,j+1,0]（但接近Y1[i,j,1]，但这对你的情况没有帮助）。你可以通过观看来看到它Y1.strides

    Y1.strides
# (1757544, 1608, 8)

    它告诉你每个轴上两个连续值之间的字节数。你会发现它比所有轴上的典型缓存大小都大，但最后一个轴除外，最后一个轴是你不使用的

    虽然我的Y1

    (8744, 8, 43720000)

    当然，问题是，当你将问题简化为单个缓慢的部分时，你可能会得出结论，你应该Y1像我一样编码。

    但我认为您的代码不会分配 201 个数字，也不会使用其他 200 个数字。换句话说，您的代码中未显示的其他部分可能使用了第 3 个轴。

    因此，也许通过按最佳顺序对代码的这一部分进行排序而获得的提升Y1必须通过代码其他部分更慢的计算来付出。

    最后一点：在进行这种计算时，重要的是避免只运行一次。因为，这也是缓存的原因。第一个算法有偏见，因为它必须读取 W，而其他两个算法可能会发现它已经在缓存中等待（可能不是你的情况，因为你的数据太大了。但对于较小的数据，你会得出结论，无论第一种方法是什么，第一种方法都更慢，因为它是付出将数据加载到缓存中的代价的方法

    • 1

  2. 如果要比较各种方法的性能，则需要以原子方式考虑操作。 我会考虑两种情况：

    1. 切片复制点
    Y2 = Y1[:,:,0]
# get slice time
Y2 = Y2.copy()
# get copy time
Y2 = Y2.dot(W)
# get dot time
    1. 切片点
    Y2 = Y1[:,:,0]
# get slice time
Y2 = Y2.dot(W)
# get dot time

    这会告诉您时间是花在切片、复制还是点上。我怀疑对于大型数组来说，复制是开销最大的部分。它还会告诉您dot整个数组和切片的执行情况是否不同。

    一旦你知道了瓶颈在哪里，你就可以精准地找出问题所在，从而加快这一部分的速度。

    • 0

  3. 您可以使用爱因斯坦求和来稍微加快这一过程。

    np.einsum('ijk,jl->ilk', Y1, W)

    其中ijk是 的形状Y1，jl是 的形状W。这会产生一个维度为 的数组(5000, 30, 201)。在我的 Macbook 上，此操作需要 157 秒，这比切片 201 次要快得多。

    • 0

相关问题