为什么 scipy.signal.correlate2d() 使用如此低效的算法？

Question

几天前，我启动了一个脚本，利用 scipy.signal.correlate2d 来计算 4096x4096 图像的二维自相关。确切的调用是

zauto = signal.correlate2d(image, image, mode='full', boundary='wrap')

三天后，它仍在运行，看不到尽头。我最终意识到它一定是在进行逐个元素的强力卷积，这个过程与 N^2 一起使用，因此 4096^4 = 281 万亿次乘法和加法。

与此同时，我最终发现，通过对图像进行二维 FFT，将其转换为二维功率谱，然后进行逆 FFT，可以获得所需的结果；IE，

image -= np.mean(image)            # remove constant bias
zfft = np.fft.fft2(image)          # take 2-D FFT (complex)
zpower = zfft*np.conjug(zfft)      # convert to power spectrum
zauto = np.real(np.ifft2(zpower))  # take inverse FFT
zauto /= zauto[0,0]                # normalize

上述几行只需不到一分钟即可完成。

我的问题：为什么 scipy.signal.correlate2d 不至少包含在可行的情况下使用更有效的算法的选项，而不是让用户发现它根本无法在更大的图像上使用的困难方式？

Answer 1

恭喜，您重新发现了卷积定理！

开玩笑吧，scipy 确实可以让您选择在信号域或傅立叶域中进行卷积，只是不使用您选择的显式 2D 方法。该signal.correlate函数将处理 N 维卷积，并且将尝试为您选择一种好的方法（参数method="auto"），或者您可以强制它使用您想要的方法。

但请注意，傅里叶方法并不总是最有效的。虽然基于 FFT 的卷积确实比直接方法具有更好的渐近复杂度，但常数因子更大。执行所涉及的 FFT 需要大量的簿记和设置，而直接方法非常简单。"auto"这就是该选项和功能的原因choose_conv_method：“最佳”方法取决于输入的大小。

至于为什么？我无法真正回答这个问题，只能说库作者可能认为大多数使用这些工具的人都会意识到其中的权衡。该scipy.signal模块提供了许多类似的方法，并且文档的几个部分明确了基于 FFT 的卷积并不是万能药。