如何使用数学将 NumPy 函数的性能提高 50%

一、说明

        2D 傅里叶变换是本世纪最重要的计算机科学算法之一。它已在我们的日常生活中得到应用，从Instagram过滤器到MP3文件的处理。

        普通用户最常用的实现，有时甚至是在不知不觉中，是 NumPy 的改编。然而，尽管它很受欢迎，但他们的算法并不是最有效的。通过一些简单的操作和 2015 年的一篇文章，我们在性能上击败了 NumPy 算法 30-60%。当前实现的核心问题是一个简单的事实，即它最初是从性能弱算法派生的。

二、NumPy实现的算法

        从本质上讲，NumPy实现的算法将常规的一维FFT依次应用于二维，这显然不能成为最优解。

        相反，在2015年，两位俄罗斯科学家提出了他们的算法版本，将一维蝶蝶变换的想法应用于二维信号。我们通过添加我们的想法有效地实现了他们的基本算法概念。

        在构建了本文中的朴素算法后，我们继续进行优化，如下所示：

void _fft2d( /* Square matrix of size N */ ) {
  // base case {
  if (N == 1) return;
  // } base case 
  int n = N >> 1;
    
  /* pseudo code {
  
  ...Creating 4 temprorary matrices here...
    
  // X(x, y, i, j)
  // x, y -- indexing over temporary submatricies
  // i, j -- indexing over rows, columns in each submatrix
 
  _fft2d(&X(0, 0), root * root, ...);
  _fft2d(&X(0, 1), root * root, ...);
  _fft2d(&X(1, 0), root * root, ...);
  _fft2d(&X(1, 1), root * root, ...);
  } pseudo code */
  for (int i = 0; i  
任何递归算法都可以通过增加基本情况的大小来增强。以下是我们的处理方式：
 
void _fft2d( /* Square matrix of size N */ ) {
  // base case {
  if (N == 1) return;
  if (N == 2) {
#define Y(y, x) (V[(y)*rowsize + (x)])
    auto x00 = Y(0, 0);
    auto x10 = Y(1, 0);
    auto x01 = Y(0, 1);
    auto x11 = Y(1, 1);
    Y(0, 0) = x00 + x10 + x01 + x11;
    Y(0, 1) = x00 + x10 - x01 - x11;
    Y(1, 0) = x00 - x10 + x01 - x11;
    Y(1, 1) = x00 - x10 - x01 + x11;
    return;
  }
  // } base case 
  // ...
}
 
        进一步的逻辑步骤是消除在每个递归步骤中创建四个不必要的临时子矩阵，而支持单个子矩阵。为此，我们使用了 algorithmica 文章中的概念，并将其修改为二维矩阵。此功能还有助于我们减少不必要的分配并增加缓存命中次数。
 
// Computing values for rev_bits[n]
auto revbits = [](size_t *v, size_t n) {
  int lg_n = log2(n);
  forn(i, n) {
    int revi = 0;
    forn(l, lg_n) revi |= ((i >> l) & 1) > 1;
#define X(y, x, i, j) (V[((y)*n + (i)) * rowsize + ((x)*n) + j])
#define params n, rowsize, W, (step  M) {
    int n = N >> 1;
#define Y(y, i, j) (V[((y)*n + (i)) * rowsize + j])
#define params n, M, rowsize, W, (step_i