维度如何影响 pyfftw 中的性能?
How do dimensions affect performance in pyfftw?
我正在尝试使用 FFT 和 pyfftw 实现 3d 卷积。我在 SO:
的另一个 post 中使用代码 posted 作为基础
class CustomFFTConvolution(object):
def __init__(self, A, B, threads=1):
shape = (np.array(A.shape) + np.array(B.shape))-1
#shape=np.array(A.shape) - np.array(B.shape)+1
if np.iscomplexobj(A) and np.iscomplexobj(B):
self.fft_A_obj = pyfftw.builders.fftn(
A, s=shape, threads=threads)
self.fft_B_obj = pyfftw.builders.fftn(
B, s=shape, threads=threads)
self.ifft_obj = pyfftw.builders.ifftn(
self.fft_A_obj.get_output_array(), s=shape,
threads=threads)
else:
self.fft_A_obj = pyfftw.builders.rfftn(
A, s=shape, threads=threads)
self.fft_B_obj = pyfftw.builders.rfftn(
B, s=shape, threads=threads)
self.ifft_obj = pyfftw.builders.irfftn(
self.fft_A_obj.get_output_array(), s=shape,
threads=threads)
def __call__(self, A, B):
s1=np.array(A.shape)
s2=np.array(B.shape)
fft_padded_A = self.fft_A_obj(A)
fft_padded_B = self.fft_B_obj(B)
ret= self.ifft_obj(fft_padded_A * fft_padded_B)
return self._centered(ret, s1 - s2 + 1)
def _centered(self,arr, newshape):
# Return the center newshape portion of the array.
newshape = np.asarray(newshape)
currshape = np.array(arr.shape)
startind = (currshape - newshape) // 2
endind = startind + newshape
myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
return arr[tuple(myslice)]
我的数据 A 的形状为 (931, 411, 806),我的过滤器 B 的形状为 (32, 32, 32)。如果我 运行 此代码在 24 核机器中使用 24 个线程,则操作需要 263 秒。
现在,如果我 运行 在同一台机器上进行相同的实验,但这次 A 的形状为 (806, 411, 931) 只是轴的交换 ,代码采用只有16秒。这是什么原因?
是否有获得最佳性能的经验法则?也许填充其中一个维度?
谢谢!
既然考虑了填充,填充的大小是否可以增加到偶数,或者小质数的倍数?选择偶数可以将挂钟时间除以3.
根据维度,某些 DFT 算法可能不可用或效率不高。
例如,执行 DFT 最有效的算法之一是 Cooley-Tuckey algorithm. It consist in dividing the DFT of a signal of composite size N=N1*N2 into N1 DTFs of size N2. As a consequence, it works better for composite sizes obtained by multiplying small prime factors (2, 3, 5, 7) for which dedicated efficient algorithms are provided in FFTW. From the documentation of FFTW:
For example, the standard FFTW distribution works most efficiently for arrays whose size can be factored into small primes (2, 3, 5, and 7), and otherwise it uses a slower general-purpose routine. If you need efficient transforms of other sizes, you can use FFTW’s code generator, which produces fast C programs (“codelets”) for any particular array size you may care about. For example, if you need transforms of size 513 = 19*33, you can customize FFTW to support the factor 19 efficiently.
您的填充尺寸具有高质因数:
931=>962=2*13*37
411=>442=2*13*17
806=>837=3*3*3*31
可以扩展填充以更接近具有小质数的数字,例如 980、448 和 864。然而,填充 3D 图像会导致内存占用量显着增加,以至于并非总是可行。
为什么改变维度的顺序会改变计算时间?
差异可能是由于输入数组是真实的。 因此,在一个维度上执行 R2C DFT,然后在第二个和第三个维度上执行 C2C 以计算 3D DFT。如果要变换的第一个维度的大小是偶数,那么R2C变换可以变成一半大小的复数DFT,如图here。此技巧不适用于奇数大小。因此,一些快速算法可能会随着 962 和 837 的翻转而变得可用。
这里是测试它的代码:
import pyfftw
import matplotlib.pyplot as plt
import multiprocessing
import numpy as np
from timeit import default_timer as timer
def listofgoodsizes():
listt=[]
p2=2
for i2 in range(11):
p3=1
for i3 in range(7):
p5=1
for i5 in range(2):
listt.append(p2*p3*p5)
p5*=5
p7=1
for i7 in range(2):
listt.append(p2*p3*p7)
p7*=7
p3*=3
p2*=2
listt.sort()
return listt
def getgoodfftwsize(n,listt):
for i in range(len(listt)):
if listt[i]>=n:
return listt[i]
return n
def timea3DR2CDFT(n,m,p):
bb = pyfftw.empty_aligned((n,m, p), dtype='float64')
bf= pyfftw.empty_aligned((n,m, (p/2+1)), dtype='complex128')
pyfftw.config.NUM_THREADS = 1 #multiprocessing.cpu_count()
fft_object_b = pyfftw.FFTW(bb, bf,axes=(0,1,2))
print n,m,p
start = timer()
fft_object_b(bb)
end = timer()
print end - start
#three prime numbers !
n=3*37
m=241
p=5*19
timea3DR2CDFT(n,m,p)
# to even size :
neven=2*((n+1)/2)
meven=2*((m+1)/2)
peven=2*((p+1)/2)
timea3DR2CDFT(neven,meven,peven)
#to nearest multiple of prime
listt=listofgoodsizes()
ngood=getgoodfftwsize(n,listt)
mgood=getgoodfftwsize(m,listt)
pgood=getgoodfftwsize(p,listt)
timea3DR2CDFT(ngood,mgood,pgood)
在我的电脑上,它打印:
111 241 95
0.180601119995
112 242 96
0.0560319423676
112 252 96
0.0564918518066
我正在尝试使用 FFT 和 pyfftw 实现 3d 卷积。我在 SO:
的另一个 post 中使用代码 posted 作为基础class CustomFFTConvolution(object):
def __init__(self, A, B, threads=1):
shape = (np.array(A.shape) + np.array(B.shape))-1
#shape=np.array(A.shape) - np.array(B.shape)+1
if np.iscomplexobj(A) and np.iscomplexobj(B):
self.fft_A_obj = pyfftw.builders.fftn(
A, s=shape, threads=threads)
self.fft_B_obj = pyfftw.builders.fftn(
B, s=shape, threads=threads)
self.ifft_obj = pyfftw.builders.ifftn(
self.fft_A_obj.get_output_array(), s=shape,
threads=threads)
else:
self.fft_A_obj = pyfftw.builders.rfftn(
A, s=shape, threads=threads)
self.fft_B_obj = pyfftw.builders.rfftn(
B, s=shape, threads=threads)
self.ifft_obj = pyfftw.builders.irfftn(
self.fft_A_obj.get_output_array(), s=shape,
threads=threads)
def __call__(self, A, B):
s1=np.array(A.shape)
s2=np.array(B.shape)
fft_padded_A = self.fft_A_obj(A)
fft_padded_B = self.fft_B_obj(B)
ret= self.ifft_obj(fft_padded_A * fft_padded_B)
return self._centered(ret, s1 - s2 + 1)
def _centered(self,arr, newshape):
# Return the center newshape portion of the array.
newshape = np.asarray(newshape)
currshape = np.array(arr.shape)
startind = (currshape - newshape) // 2
endind = startind + newshape
myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
return arr[tuple(myslice)]
我的数据 A 的形状为 (931, 411, 806),我的过滤器 B 的形状为 (32, 32, 32)。如果我 运行 此代码在 24 核机器中使用 24 个线程,则操作需要 263 秒。 现在,如果我 运行 在同一台机器上进行相同的实验,但这次 A 的形状为 (806, 411, 931) 只是轴的交换 ,代码采用只有16秒。这是什么原因? 是否有获得最佳性能的经验法则?也许填充其中一个维度? 谢谢!
既然考虑了填充,填充的大小是否可以增加到偶数,或者小质数的倍数?选择偶数可以将挂钟时间除以3.
根据维度,某些 DFT 算法可能不可用或效率不高。 例如,执行 DFT 最有效的算法之一是 Cooley-Tuckey algorithm. It consist in dividing the DFT of a signal of composite size N=N1*N2 into N1 DTFs of size N2. As a consequence, it works better for composite sizes obtained by multiplying small prime factors (2, 3, 5, 7) for which dedicated efficient algorithms are provided in FFTW. From the documentation of FFTW:
For example, the standard FFTW distribution works most efficiently for arrays whose size can be factored into small primes (2, 3, 5, and 7), and otherwise it uses a slower general-purpose routine. If you need efficient transforms of other sizes, you can use FFTW’s code generator, which produces fast C programs (“codelets”) for any particular array size you may care about. For example, if you need transforms of size 513 = 19*33, you can customize FFTW to support the factor 19 efficiently.
您的填充尺寸具有高质因数:
931=>962=2*13*37
411=>442=2*13*17
806=>837=3*3*3*31
可以扩展填充以更接近具有小质数的数字,例如 980、448 和 864。然而,填充 3D 图像会导致内存占用量显着增加,以至于并非总是可行。
为什么改变维度的顺序会改变计算时间? 差异可能是由于输入数组是真实的。 因此,在一个维度上执行 R2C DFT,然后在第二个和第三个维度上执行 C2C 以计算 3D DFT。如果要变换的第一个维度的大小是偶数,那么R2C变换可以变成一半大小的复数DFT,如图here。此技巧不适用于奇数大小。因此,一些快速算法可能会随着 962 和 837 的翻转而变得可用。
这里是测试它的代码:
import pyfftw
import matplotlib.pyplot as plt
import multiprocessing
import numpy as np
from timeit import default_timer as timer
def listofgoodsizes():
listt=[]
p2=2
for i2 in range(11):
p3=1
for i3 in range(7):
p5=1
for i5 in range(2):
listt.append(p2*p3*p5)
p5*=5
p7=1
for i7 in range(2):
listt.append(p2*p3*p7)
p7*=7
p3*=3
p2*=2
listt.sort()
return listt
def getgoodfftwsize(n,listt):
for i in range(len(listt)):
if listt[i]>=n:
return listt[i]
return n
def timea3DR2CDFT(n,m,p):
bb = pyfftw.empty_aligned((n,m, p), dtype='float64')
bf= pyfftw.empty_aligned((n,m, (p/2+1)), dtype='complex128')
pyfftw.config.NUM_THREADS = 1 #multiprocessing.cpu_count()
fft_object_b = pyfftw.FFTW(bb, bf,axes=(0,1,2))
print n,m,p
start = timer()
fft_object_b(bb)
end = timer()
print end - start
#three prime numbers !
n=3*37
m=241
p=5*19
timea3DR2CDFT(n,m,p)
# to even size :
neven=2*((n+1)/2)
meven=2*((m+1)/2)
peven=2*((p+1)/2)
timea3DR2CDFT(neven,meven,peven)
#to nearest multiple of prime
listt=listofgoodsizes()
ngood=getgoodfftwsize(n,listt)
mgood=getgoodfftwsize(m,listt)
pgood=getgoodfftwsize(p,listt)
timea3DR2CDFT(ngood,mgood,pgood)
在我的电脑上,它打印:
111 241 95
0.180601119995
112 242 96
0.0560319423676
112 252 96
0.0564918518066