numpy 为什么复共轭如此缓慢?

mspsb9vt  于 2023-05-17  发布在  其他
关注(0)|答案(3)|浏览(124)

它和复数乘法是一样的,这让我很困惑:

import numpy as np

def op0(x):
    return np.conj(x)
    
def op1(x0, x1):
    return x0 * x1
    
def op2(x0, x1):
    x0[:] = x1

for N in (50, 500, 5000):
    print(f"\nshape = ({N}, {N})")
    x0 = np.random.randn(N, N) + 1j*np.random.randn(N, N)
    x1 = np.random.randn(N, N) + 1j*np.random.randn(N, N)

    %timeit op0(x0)
    %timeit op1(x0, x1)
    %timeit op2(x0, x1)
shape = (50, 50)
3.55 µs ± 143 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
4.85 µs ± 261 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
1.85 µs ± 116 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

shape = (500, 500)
1.52 ms ± 60.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
1.96 ms ± 133 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
299 µs ± 50.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

shape = (5000, 5000)
163 ms ± 4.4 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
185 ms ± 11.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
39.8 ms ± 399 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

为什么翻转x.imag的标志这么贵?当然,在低层次上,它比几个乘法和加法((a + j*b)*(c + j*d))容易得多?
Windows 10 x64、numpy 1.23.5、Python 3.10.4

ndasle7k

ndasle7k1#

在记忆中来回移动需要花费大部分时间。
我想出了类似于莎拉梅塞尔的想法(在上面的评论中)。仅复制输入数组就需要花费大量时间:

def op0(x):
    return np.conj(x)

def op3(x):
    return x.copy()

这些操作之间的时间差约为20%。

就地共轭

正如杰罗姆在评论中提到的,使用out参数:

def op0B(x):
    np.conj(x, out=x)

对于N=5000,在我的机器上快3倍。
对于小的N,改进是可忽略的。有趣的是,对于N=500,加速甚至更好(5-6x)。

iqjalb3h

iqjalb3h2#

测试更充分,完整的结果在底部。主要成果:

op4 = lambda x:      np.conj(x, out=x)
op5 = lambda x0, x1: np.multiply(x0, x1, out=x0)

def op6(x0, x1): 
    x0[:] = x1

得分

shape = (5000, 5000)
4: 43.5 ms ± 2.63 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
5: 71.7 ms ± 4.15 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
6: 43.3 ms ± 3.27 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

NumPy中的operation(x, out=x)使其就位,如果错了请纠正我。将out=out与预分配的out一起使用会更慢(见下文)。所以...我想是被说服了
不确定我是否应该对复数乘法的速度印象深刻,或者对基本的数组东西的速度感到失望。

操作分解

S = N**2

  • op4S读取,S写入--S符号翻转
  • op54*S读取,2*S写入--4*S乘法,4*S加法
  • op62*S读取,2*S写入

op4op6S读取,S写入与S符号翻转相当-似乎公平。
op5op6:额外的2*S读取,4*S乘法,4*S加法甚至不会使计算时间加倍。写作是最贵的吗?
如果我的分析出错了就告诉我。

完整结果

import numpy as np

op0 = lambda x:           np.conj(x)
op1 = lambda x0, x1:      x0 * x1
op2 = lambda x, out:      np.conj(x, out=out)
op3 = lambda x0, x1, out: np.multiply(x0, x1, out=out)
op4 = lambda x:           np.conj(x, out=x)
op5 = lambda x0, x1:      np.multiply(x0, x1, out=x0)

def op6(x0, x1):
    x0[:] = x1

for N in (50, 500, 5000):
    print(f"\nshape = ({N}, {N})")
    out = np.zeros((N, N), dtype='complex128')
    x0  = np.ones((N, N),  dtype='complex128')
    x1 = x0.copy()

    %timeit op0(x0)
    %timeit op1(x0, x1)
    %timeit op2(x0, out)
    %timeit op3(x0, x1, out)
    %timeit op4(x0)
    %timeit op5(x0, x1)
    %timeit op6(x0, x1)
shape = (50, 50)
0: 4.32 µs ± 160  ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
1: 4.81 µs ± 280  ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
2: 2.83 µs ± 28.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
3: 4.25 µs ± 293  ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
4: 2.69 µs ± 38.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
5: 3.98 µs ± 67.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
6: 1.78 µs ± 29.6 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

shape = (500, 500)
0: 1.43 ms ± 94.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
1: 1.48 ms ± 26.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
2: 401 µs  ± 32.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
3: 845 µs  ± 66.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
4: 252 µs  ± 11.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
5: 532 µs  ± 17.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
6: 224 µs  ± 15.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

shape = (5000, 5000)
0: 165 ms  ± 3.26 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
1: 190 ms  ± 2.82 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
2: 62.4 ms ± 2.52 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
3: 87.7 ms ± 3.41 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
4: 43.5 ms ± 2.63 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
5: 71.7 ms ± 4.15 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
6: 43.3 ms ± 3.27 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
pw9qyyiw

pw9qyyiw3#

以下是基于许多以前的评论的答案:
不确定我是否应该对复数乘法的速度印象深刻,或者对基本的数组东西的速度感到失望。
内存访问很慢,并且由于内存墙(开始于>20年前),将来会更慢。如果你想要快速的东西,你不需要使用你的DRAM。Numpy并不是专门为此设计的(这是相当可悲的)。您的硬件达到17.2 GiB/s。这不是很多,但也不是很糟糕。现代PC可以达到40-60 GiB/s(有些,如M1,甚至可以达到>60 GiB)。现代英特尔CPU L3缓存可以达到~400 GiB/s,所以更多。
GPU通常具有显著更高的内存吞吐量,但比率computational_speed/memory_bandwidth仍然像CPU一样高(实际上甚至更高)。CPU现在有相当大的缓存,而GPU通常没有。请注意,GPU计算可能会被某些API延迟(惰性计算),因此您应该在基准测试中关注这一点(例如,您可以打印值)。
op 5 vs op 6:额外的2S读取,4S乘法,4*S加法甚至不会使计算时间加倍。写作是最贵的吗?
这里所有的操作都应该是内存受限的。只有内存操作才重要。
op5op6慢,因为op5读取两个数组,而op6读取一个。需要从DRAM传输更多数据,因此需要更多时间。此外,写可能比读更昂贵,但这取决于编译的汇编代码(因此编译器,优化标志和实际源代码)和目标体系结构。内存性能是一个复杂的主题,比看起来要复杂得多(有关此方面的更多信息,请参阅下文)。
请注意,它是两个独立的数组并不重要。一个大阵列实际上分为两部分将具有相同的影响。与硬件没有太大区别。

一般注意事项

关于计时,现代内存/CPU非常复杂,因此通常不容易理解基于基准测试的情况(特别是没有汇编代码)。从DRAM的Angular 来看,写入和读取几乎同样快。由于DRAM的工作方式,混合它们会降低性能。现代x86-64 CPU缓存使用写分配策略,导致在执行写操作时执行读操作。这导致写入比读取慢2倍。也就是说,假设编译器生成了非时态指令,那么可以使用非时态指令来避免这个问题
编译器通常不会生成非临时指令,因为当数组适合该高速缓存时,它们可能会慢得多,在这种情况下,用于构建Numpy的编译器无法知道数组的大小(运行时定义)。它也无法知道CPU缓存的大小。内存复制倾向于使用memcpy基本函数,该函数经过优化,可针对目标平台的大型数组使用非临时指令。AFAIK,这样的指令在Numpy中不用于乘法/加法等操作。(至少在我的机器上不是)。
请注意,我提到“现代CPU”是因为Zen之前的AMD CPU使用了一种非常不同的方法。在这种情况下,其他CPU体系结构(如POWER)的行为也非常不同。在高性能方面,任何细节都很重要。这也是为什么我说这个主题很复杂。最好的方法是在您的特定机器上进行低级分析,或者列出您的硬件,确切的Numpy包(或使用的汇编代码),以避免考虑理论上可能发生的许多可能的事情。

相关问题