矩阵矩阵乘法的 CPU 时间答案

【问题标题】：CPU Time for matrix matrix multiplication矩阵矩阵乘法的 CPU 时间
【发布时间】：2019-08-09 19:32:42
【问题描述】：

我正在尝试决定是否应该同时或顺序处理几个相似但独立的问题（可能在不同的计算机上并行处理）。为了做出决定，我需要比较以下操作的 cpu 时间：

time_1 是计算 X(with shape (n,p)) @ b (with shape (p,1)) 的时间。
time_k 是计算X(形状为(n,p))@B(形状为(p,k))的时间。

其中 X、b 和 B 是随机矩阵。这两个操作之间的区别在于第二个矩阵的宽度。

天真地，我们期望 time_k = k x time_1。使用更快的矩阵乘法算法（Strassen 算法、Coppersmith–Winograd 算法），time_k 可能小于 k x time_1，但这些算法的复杂度仍然比我在实践中观察到的要大得多。因此我的问题是：如何解释这两种计算在 cpu 时间方面的巨大差异？

我使用的代码如下：

import time
import numpy as np
import matplotlib.pyplot as plt

p     = 100
width = np.concatenate([np.arange(1, 20), np.arange(20, 100, 10), np.arange(100, 4000, 100)]).astype(int)

mean_time = []
for nk, kk in enumerate(width):
    timings = []
    nb_tests = 10000 if kk <= 300 else 100
    for ni, ii in enumerate(range(nb_tests)):
        print('\r[', nk, '/', len(width), ', ',  ni, '/', nb_tests, ']', end = '')
        x     = np.random.randn(p).reshape((1, -1))
        coef  = np.random.randn(p, kk)
        d     = np.zeros((1, kk))
        start = time.time()
        d[:]  = x @ coef
        end   = time.time()
        timings.append(end - start)

    mean_time.append(np.mean(timings))

mean_time = np.array(mean_time)


fig, ax = plt.subplots(figsize =(14,8))
plt.plot(width, mean_time, label =  'mean(time\_k)')
plt.plot(width, width*mean_time[0], label = 'k*mean(time\_1)')
plt.legend()
plt.xlabel('k')
plt.ylabel('time (sec)')
plt.show()

【问题讨论】：

你的意思是什么？向我们展示您的困惑。

标签： python numpy matrix time multiplication

【解决方案1】：

你不仅仅是时间乘法运算。 time.time() 需要时间才能完成。

>>> print(time.time() - time.time())
-9.53674316406e-07

当乘以尝试次数 (10000) 时，实例的数量就会变得很大开销，对于 n=100，您实际上是将 1.000.000 次调用 time.time() 与 100 次常规 numpy 数组乘法进行比较。

为了快速进行基准测试，Python 提供了一个不存在此问题的专用模块：请参阅timeit

【讨论】：

我同意 timeit 是一个更好的计时工具，但我不认为这是造成差异的原因，特别是因为我计算了“nb_tests”试验的平均时间，所以偏移量不是总结。我尝试使用 timeit 并将内部循环替换为以下内容：mean_time.append(timeit.timeit(stmt = 'd[:] = x @ coef', setup = 'import numpy as np; p = 100; x = np.random.randn(p).reshape((1, -1)); coef = np.random.randn(p, {width_k}); d = np.zeros((1, {width_k}));'.format(width_k = kk),))。观察结果保持不变。
通过删除内部循环并向 timeit 添加一个数字参数，我在 4000 端获得了橙色曲线和蓝色曲线之间的 60:1 比率，这听起来很合理。毕竟逐位相乘比较慢。

【解决方案2】：

这个细节的原因很复杂。你知道当PC运行X @ b时，它会执行许多其他需要的指令，可能是load data from RAM to cache等等。换句话说，成本时间包含两部分——Cost_A 表示的 CPU 中的“实际计算指令”和Cost_B 表示的“其他需要的指令”。我有一个想法，只是我的猜测，是Cost_B 导致time_k << k x time_1。

由于 b 的形状很小（例如 1000 x 1），“其他需要的指令”相对花费的时间最多。因为 b 的形状很大（例如 1000 x 10000），所以它相对较小。下面的一组实验可以给出一个不太严格的证明。我们可以看到，当 b 的形状从 (1000 x 1) 增加到 (1000 x ) 时，成本时间增加非常缓慢。

import numpy as np
import time

X = np.random.random((1000, 1000))

b = np.random.random((1000, 1))
b3 = np.random.random((1000, 3))
b5 = np.random.random((1000, 5))
b7 = np.random.random((1000, 7))
b9 = np.random.random((1000, 9))
b10 = np.random.random((1000, 10))
b30 = np.random.random((1000, 30))
b60 = np.random.random((1000, 60))
b100 = np.random.random((1000, 100))
b1000 = np.random.random((1000, 1000))

def test_cost(X, b):
    begin = time.time()
    for i in range(100):
        _ = X @ b
    end = time.time()
    print((end-begin)/100.)

test_cost(X, b)
test_cost(X, b3)
test_cost(X, b5)
test_cost(X, b7)
test_cost(X, b9)
test_cost(X, b10)
test_cost(X, b30)
test_cost(X, b60) 
test_cost(X, b100)
test_cost(X, b1000)

output:
0.0003210139274597168
0.00040063619613647463
0.0002452659606933594
0.00026523590087890625
0.0002449488639831543
0.00024344682693481446
0.00040068864822387693
0.000691361427307129
0.0011700797080993653
0.009680757522583008

更多，我在linux中用pref做了一组实验。对于pref，Cost_B 可能更大。我有8个python文件，第一个如下。

import numpy as np
import time
def broken2():
    mtx = np.random.random((1, 1000))
    c = None
    c = mtx ** 2

broken2()

我已将输出处理到表 A，如下所示。

我做了一个简单的分析，我将邻居实验中操作数（likes，cache-misses）的误差除以time elapsed(seconds)的误差。然后，我得到下表B。从表中我们可以发现，随着b的形状增加，形状与成本时间之间的线性关系更加明显。而导致time_k << k x time_1的主要原因可能是cache misses（从RAM加载数据到缓存），首先稳定了。

【讨论】：