正确地使用pyfftw来加速numpy

Question

我正试图从Matlab跳到numpy，但我迫切需要fft的速度，现在我知道pyfftw，但我不知道我是否正确地使用了它。我的方法是

import numpy as np
import pyfftw
import timeit

pyfftw.interfaces.cache.enable()

def wrapper(func, *args):
    def wrapped():
        return func(*args)
    return wrapped

def my_fft(v):
    global a
    global fft_object
    a[:] = v
    return fft_object()

def init_cond(X):
    return my_fft(2.*np.cosh(X)**(-2))

def init_cond_py(X):
    return np.fft.fft(2.*np.cosh(X)**(-2))

K = 2**16
Llx = 10.
KT = 2*K
dx = Llx/np.float64(K)
X = np.arange(-Llx,Llx,dx)

global a
global b
global fft_object
a = pyfftw.n_byte_align_empty(KT, 16, 'complex128')
b = pyfftw.n_byte_align_empty(KT, 16, 'complex128')
fft_object = pyfftw.FFTW(a,b)

wrapped = wrapper(init_cond, X)
print min(timeit.repeat(wrapped,repeat=100,number=1))

wrapped_two = wrapper(init_cond_py, X)
print min(timeit.repeat(wrapped_two,repeat=100,number=1))

我知道这里有构建器函数，也有通过pyfftw对fft和numpy fft调用的标准接口。不过，这些都表现得非常缓慢。通过首先创建一个fft_object实例，然后在全局使用它，我已经能够获得与numpy的快速或稍快的速度。

话虽如此，我工作的前提是智慧被暗藏起来。这是真的吗？我需要把这个说清楚吗？如果是的话，最好的方法是什么？

而且，我认为时间是完全不透明的。我用得好吗？它是否像我所说的那样储存智慧？提前感谢你能给予的任何帮助。

Answer 1

在交互式(ipython)会话中，我认为下面是您想要做的事情(ipython很好地处理了时间问题)：

In [1]: import numpy as np

In [2]: import pyfftw

In [3]: K = 2**16

In [4]: Llx = 10.

In [5]: KT = 2*K

In [6]: dx = Llx/np.float64(K)

In [7]: X = np.arange(-Llx,Llx,dx)

In [8]: a = pyfftw.n_byte_align_empty(KT, 16, 'complex128')

In [9]: b = pyfftw.n_byte_align_empty(KT, 16, 'complex128')

In [10]: fft_object = pyfftw.FFTW(a,b)

In [11]: a[:] = 2.*np.cosh(X)**(-2)

In [12]: timeit np.fft.fft(a)
100 loops, best of 3: 4.96 ms per loop

In [13]: timeit fft_object(a)
100 loops, best of 3: 1.56 ms per loop

In [14]: np.allclose(fft_object(a), np.fft.fft(a))
Out[14]: True

你读过教程吗？你不明白什么？

我建议使用建筑商界面来构造FFTW对象。播放各种设置，最重要的是线程的数量。

智慧不是默认存储的。你需要自己把它拿出来。

您的所有globals都是不必要的--您想要更改的对象是可变的，所以您可以很好地处理它们。fft_object总是指出同样的事情，所以这不是一个全球性的问题。理想情况下，您只是不希望ii上的循环。我建议如何构造数组，以便您可以在一个调用中完成所有操作。

编辑:编辑:我只粗略地看了一下你的代码就写了下面的一段，很明显，它是一个递归的更新，矢量化并不是一种明显的方法，没有一些严重的狡猾。我在底层对您的实现有一些评论，不过我怀疑您的问题是对如何最好地使用Python (或Matlab)这样的语言进行数值处理是一个更根本的误解。核心原则是尽可能地向矢量化。我的意思是把您的python调用滚到尽可能少的地方。不幸的是，我不知道如何使用您的示例(虽然我只考虑了2分钟)。如果这仍然失败，请考虑cython --尽管要确保您真的想走这条路(也就是说，您已经用尽了其他选项)。

关于全球:不要这样做。如果您想要创建一个带有状态的对象，可以使用一个类(这就是它们的用途)，或者在您的情况下使用一个闭包。全局几乎从来都不是您想要的(我认为我在编写python的所有文章中至少都有一个是合法的，在pyfftw的缓存代码中)。我建议阅读这么好的问题。Matlab是一门蹩脚的语言--这其中的一个原因是它糟糕的范围界定设施，这往往会导致坏习惯。

只有在要全局修改引用时才需要全局引用。我建议阅读更多关于Python作用域规则和python中哪些变量真的是的内容。

FFTW对象随身携带所需的所有数组，因此不需要单独传递它们。在设置或返回值时，使用调用接口几乎不需要任何开销(特别是如果您禁用了正常化)--如果您处于优化的水平，我强烈怀疑您已经达到了极限(我要提醒您，对于许多非常小的FFTs来说，这可能不是完全正确的，但此时您需要重新考虑向FFTW调用向量化的算法)。如果您发现每次更新数组(使用调用接口)时存在大量开销，那么这是一个错误，您应该按此方式提交它(我会很惊讶)。

总之，不要担心每次调用都会更新数组。这几乎肯定不是您的瓶颈，不过要确保您知道规范化，如果您愿意，可以禁用它(与update_arrays()和execute()方法的原始访问相比，它可能会稍微放慢速度)。

您的代码没有使用缓存。只有在使用interfaces代码时才使用缓存，并减少了在内部创建新FFTW对象时的Python开销。由于您正在自己处理FFTW对象，所以没有理由使用缓存。

对于获取FFTW对象，builders代码是一个约束较小的接口。我现在几乎总是使用构建器(从头创建FFTW对象要方便得多)。您想直接创建FFTW对象的情况非常少见，我很想知道它们是什么。

关于算法实现的评论:我不熟悉您正在实现的算法。然而，我对你目前是如何写的有一些评论。您正在计算每个循环上的nl_eval(wp)，但据我所知，这与前一个循环中的nl_eval(w)完全相同，因此您不需要计算它两次(但这附带了一个警告:当您在任何地方都有全局值时，很难看到发生了什么，所以我可能遗漏了一些东西)。

不要费心使用my_fft或my_ifft中的副本。只需执行fft_object(u) (在我的机器上为前向情况做2.29ms和1.67ms)。内部数组更新例程使复制变得不必要。而且，正如您编写的那样，您要复制两次：c[:]的意思是“复制到数组c"，而要复制到c的数组是v.copy()，即v的一个副本(总共两个拷贝)。

更明智的方法(可能也是必要的)是将输出复制到保持数组中(因为这样可以避免调用FFTW对象的临时结果)，不过要确保保持数组是正确对齐的。我相信您已经注意到这一点很重要，但是复制输出更容易理解。

你可以把你所有的头皮都移到一起。wn计算中的3可移动到nl_eval的my_fft中。您还可以将它与ifft中的正常化常量结合起来(并在pyfftw中关闭它)。

看看基本数组操作的numexpr。它可以提供相当多的速度比香草泡菜。

不管怎么说，从这一切中拿出你想要的。毫无疑问，我错过了一些东西，或者说了一些不正确的话，所以请尽可能谦卑地接受它。与Matlab相比，花点时间研究Python是如何运行的是值得的(事实上，忘记后者就好了)。