python - 为什么一个代码(matmul)比另一个(Python)快

Question

该问题涉及基本矩阵运算。在下面的代码中，c1 本质上等于 c2。然而，第一种计算方式比第二种方式快得多。其实一开始我以为第一种方式需要分配一个比a矩阵大两倍的b矩阵，因此可能会更慢。结果恰恰相反。为什么？

import time
import numpy as np
a = np.random.rand(20000,100)+np.random.rand(20000,100)*1j

tic = time.time()
b = np.vstack((a.real,a.imag))
c1 = b.T @ b
t1 = time.time()-tic

tic = time.time()
c2 = a.real.T @ a.real+a.imag.T@a.imag
t2 = time.time()-tic

print('t1=%f. t2=%f.'%(t1,t2))

一个示例结果是

t1=0.037965. t2=4.375873.

Answer 1

a.real和 a.imag就地访问，而 np.vstack创建一个新副本。方式@运算符 ( matmul() ) 处理 a.real和 a.imag需要更长的时间。为了使其更快，您可以创建每个副本，然后将其传递给 @或使用 np.dot(a.real.T, a.real)和 np.dot(a.imag.T, a.imag) (我不确定 BLAS 中的每个实现)。

对于大型矩阵，以下代码中的第一种方法应该仍然稍微快一点:

a = np.random.rand(20000,100)+np.random.rand(20000,100)*1j

tic = time.time()
b = np.vstack((a.real,a.imag))
c1 = b.T @ b
t1 = time.time()-tic

tic = time.time()
b = a.real.copy()
c = a.imag.copy()
c2 = b.T @ b + c.T @ c
t2 = time.time()-tic

print('t1=%f. t2=%f.'%(t1,t2))  

输出:

t1=0.031620. t2=0.021769.

编辑:深入研究:

让我们快速浏览一下各种内存布局:

a = np.random.rand(20000,100)+np.random.rand(20000,100)*1j
print('a.flags\n', a.flags)
print('a.real flags\n', a.real.flags)


a.flags
C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : True
...
a.real flags
C_CONTIGUOUS : False
F_CONTIGUOUS : False
OWNDATA : False

所以 a是 C_CONTIGUOUS和 a.real不是。我不知道如何 @实现计算，但我的猜测是缓存技巧和步幅以及展开的循环不同。我将把它留给专家来解释。现在， array.copy()默认为 C_CONTIGUOUS (小心: np.copy() 默认情况下不是 C_CONTIGUOUS。)这就是为什么上面的第二种方法和第一种方法一样快(其中 b 也是 C_CONTIGUOUS )。

@George C:我认为第一个稍微快一点，因为 np.vstack创建一个新的 C_CONTIGUOUS可以在一个地方利用缓存技巧的对象，而在第二种方法中 a.real.T @ a.real 的输出和 a.imag.T@a.imag位于不同的内存位置，需要额外的努力来计算。 Here is a link to more explanation.
免责声明:欢迎任何了解 NumPy 实现细节的专家编辑这篇文章。