FFTW 与 Matlab FFT

Question

我在 matlab central 上发布了这个，但没有得到任何回复，所以我想我会在这里重新发布。

我最近在 Matlab 中编写了一个在 for 循环中使用 FFT 的简单例程； FFT 在计算中占主导地位。出于实验目的，我在 mex 中编写了相同的例程，它调用了 FFTW 3.3 库。事实证明，对于非常大的数组，matlab 例程比 mex 例程运行得更快（大约快两倍）。 mex 例程使用智慧并执行相同的 FFT 计算。我也知道 matlab 使用 FFTW，但他们的版本是否可能稍微优化一些？我什至使用了 FFTW_EXHAUSTIVE 标志，对于大型数组，它的速度仍然比 MATLAB 对应的慢两倍。此外，我确保我使用的 matlab 是带有“-singleCompThread”标志的单线程，并且我使用的 mex 文件未处于调试模式。只是好奇是否是这种情况-或者是否有一些我不知道的 matlab 正在使用的优化。谢谢。

这是墨西哥的部分：

void class_cg_toeplitz::analysis() {
// This method computes CG iterations using FFTs
    // Check for wisdom
    if(fftw_import_wisdom_from_filename("cd.wis") == 0) {
        mexPrintf("wisdom not loaded.\n");
    } else {
        mexPrintf("wisdom loaded.\n");
    }

    // Set FFTW Plan - use interleaved FFTW
    fftw_plan plan_forward_d_buffer;    
    fftw_plan plan_forward_A_vec;       
    fftw_plan plan_backward_Ad_buffer;
    fftw_complex *A_vec_fft;
    fftw_complex *d_buffer_fft;
    A_vec_fft = fftw_alloc_complex(n);
    d_buffer_fft = fftw_alloc_complex(n);

    // CREATE MASTER PLAN - Do this on an empty vector as creating a plane 
    // with FFTW_MEASURE will erase the contents; 
    // Use d_buffer
    // This is somewhat dangerous because Ad_buffer is a vector; but it does not
    // get resized so &Ad_buffer[0] should work
    plan_forward_d_buffer = fftw_plan_dft_r2c_1d(d_buffer.size(),&d_buffer[0],d_buffer_fft,FFTW_EXHAUSTIVE);
    plan_forward_A_vec = fftw_plan_dft_r2c_1d(A_vec.height,A_vec.value,A_vec_fft,FFTW_WISDOM_ONLY);
    // A_vec_fft.*d_buffer_fft will overwrite d_buffer_fft
    plan_backward_Ad_buffer = fftw_plan_dft_c2r_1d(Ad_buffer.size(),d_buffer_fft,&Ad_buffer[0],FFTW_EXHAUSTIVE);

    // Get A_vec_fft
    fftw_execute(plan_forward_A_vec);

    // Find initial direction - this is the initial residual
    for (int i=0;i<n;i++) {
        d_buffer[i] = b.value[i];
        r_buffer[i] = b.value[i];
    }    

    // Start CG iterations
    norm_ro = norm(r_buffer);
    double fft_reduction = (double)Ad_buffer.size(); // Must divide by size of vector because inverse FFT does not do this
    while (norm(r_buffer)/norm_ro > relativeresidual_cutoff) {        
        // Find Ad - use fft
        fftw_execute(plan_forward_d_buffer);    
        // Get A_vec_fft.*fft(d) - A_vec_fft is only real, but d_buffer_fft
        // has complex elements; Overwrite d_buffer_fft        
        for (int i=0;i<n;i++) {
            d_buffer_fft[i][0] = d_buffer_fft[i][0]*A_vec_fft[i][0]/fft_reduction;
            d_buffer_fft[i][1] = d_buffer_fft[i][1]*A_vec_fft[i][0]/fft_reduction;
        }        
        fftw_execute(plan_backward_Ad_buffer); 

        // Calculate r'*r
        rtr_buffer = 0;
        for (int i=0;i<n;i++) {
            rtr_buffer = rtr_buffer + r_buffer[i]*r_buffer[i];
        }    

        // Calculate alpha
        alpha = 0;
        for (int i=0;i<n;i++) {
            alpha = alpha + d_buffer[i]*Ad_buffer[i];
        }    
        alpha = rtr_buffer/alpha;

        // Calculate new x
        for (int i=0;i<n;i++) {
            x[i] = x[i] + alpha*d_buffer[i];
        }   

        // Calculate new residual
        for (int i=0;i<n;i++) {
            r_buffer[i] = r_buffer[i] - alpha*Ad_buffer[i];
        }   

        // Calculate beta
        beta = 0;
        for (int i=0;i<n;i++) {
            beta = beta + r_buffer[i]*r_buffer[i];
        }  
        beta = beta/rtr_buffer;

        // Calculate new direction vector
        for (int i=0;i<n;i++) {
            d_buffer[i] = r_buffer[i] + beta*d_buffer[i];
        }  

        *total_counter = *total_counter+1;
        if(*total_counter >= iteration_cutoff) {
            // Set total_counter to -1, this indicates failure
            *total_counter = -1;
            break;
        }
    }

    // Store Wisdom
    fftw_export_wisdom_to_filename("cd.wis");

    // Free fft alloc'd memory and plans
    fftw_destroy_plan(plan_forward_d_buffer);
    fftw_destroy_plan(plan_forward_A_vec);
    fftw_destroy_plan(plan_backward_Ad_buffer);
    fftw_free(A_vec_fft);
    fftw_free(d_buffer_fft);
};

这是 matlab 部分：

% Take FFT of A_vec.
A_vec_fft = fft(A_vec); % Take fft once

% Find initial direction - this is the initial residual 
x = zeros(n,1); % search direction
r = zeros(n,1); % residual
d = zeros(n+(n-2),1); % search direction; pad to allow FFT
for i = 1:n
    d(i) = b(i); 
    r(i) = b(i); 
end

% Enter CG iterations
total_counter = 0;
rtr_buffer = 0;
alpha = 0;
beta = 0;
Ad_buffer = zeros(n+(n-2),1); % This holds the product of A*d - calculate this once per iteration and using FFT; only 1:n is used
norm_ro = norm(r);

while(norm(r)/norm_ro > 10^-6)
    % Find Ad - use fft
    Ad_buffer = ifft(A_vec_fft.*fft(d)); 

    % Calculate rtr_buffer
    rtr_buffer = r'*r;

    % Calculate alpha    
    alpha = rtr_buffer/(d(1:n)'*Ad_buffer(1:n));

    % Calculate new x
    x = x + alpha*d(1:n);

    % Calculate new residual
    r = r - alpha*Ad_buffer(1:n);

    % Calculate beta
    beta = r'*r/(rtr_buffer);

    % Calculate new direction vector
    d(1:n) = r + beta*d(1:n);      

    % Update counter
    total_counter = total_counter+1; 
end

就时间而言，对于 N = 50000 和 b = 1:n，mex 大约需要 10.5 秒，matlab 大约需要 4.4 秒。我正在使用 R2011b。谢谢

Answer 1

由于我不了解 MATLAB FFT 实现的任何细节，因此仅提供一些观察而不是明确的答案：

• 根据您的代码，我可以看到速度差异的两种解释：
  • 速度差异可以通过 FFT 优化级别的差异来解释
  • MATLAB 中的 while 循环执行次数明显减少

我假设您已经研究了第二个问题并且迭代次数是可比较的。 （如果不是，这很可能是一些准确性问题，值得进一步调查。）

现在，关于 FFT 速度比较：

• 是的，理论上 FFTW 比其他高级 FFT 实现更快，但只有在将苹果与苹果进行比较时才有意义：这里是在更底层的程序集级别比较实现，其中不仅算法的选择，而且针对特定处理器和具有不同技能的软件开发人员的实际优化都在起作用
• 一年来，我在许多处理器上优化或审查了优化的 FFT 组装（我从事基准测试行业），出色的算法只是故事的一部分。有一些非常特定于您正在编码的架构的考虑因素（考虑延迟、指令调度、优化寄存器使用、内存中的数据排列、考虑分支采用/未采用延迟等）并且会产生差异与算法的选择同样重要。
• 在 N=500000 的情况下，我们还讨论了大内存缓冲区：又一扇通往更多优化的大门，可以快速针对您运行代码的平台进行调整：避免缓存未命中的能力不会算法决定了数据的流动方式以及软件开发人员可能使用了哪些优化来有效地将数据导入和导出内存。
• 虽然我不知道 MATLAB FFT 实现的细节，但我很确定 DSP 工程师大军已经（并且仍在）对其优化进行磨练，因为它是众多设计的关键。这很可能意味着 MATLAB 拥有合适的开发人员组合来生成更快的 FFT。

Answer 2

由于低级和特定于架构的优化，这是典型的性能提升。

Matlab 使用来自 Intel MKL（数学内核库）二进制文件 (mkl.dll) 的 FFT。这些是英特尔针对英特尔处理器优化（在汇编级别）的例程。即使在 AMD 上，它似乎也能带来不错的性能提升。

FFTW 看起来像一个没有优化的普通 c 库。因此使用 MKL 可以获得性能提升。

Answer 3

我在 MathWorks 网站 [1] 上发现了以下评论：

注意 2 的大幂：对于 FFT 维度，它是 2 的幂 2、2^14和2^22之间，MATLAB软件使用特殊预加载 内部数据库中的信息以优化 FFT 计算。 当 FTT 的维度是 2 的幂时，不进行调整， 除非你使用命令 fftw('wisdom', []) 清除数据库。

虽然它与 2 的幂有关，但它可能暗示 MATLAB 在将 FFTW 用于某些（大）数组大小时采用了自己的“特殊智慧”。考虑：2^16 = 65536。

[1] R2013b 文档可从http://www.mathworks.de/de/help/matlab/ref/fftw.html 获得（2013 年 10 月 29 日访问）

Answer 4

编辑： @wakjah 对此答案的回复是准确的：FFTW 确实支持通过其 Guru 接口拆分实数和虚数内存存储。因此，我关于黑客攻击的说法并不准确，但如果不使用 FFTW 的 Guru 界面，则可以很好地应用——默认情况下就是这种情况，所以还是要小心！

首先，很抱歉迟到了一年。我不相信您看到的速度提升来自 MKL 或其他优化。 FFTW 和 Matlab 之间有一些根本不同的地方，那就是复杂数据在内存中的存储方式。

在 Matlab 中，复数向量 X 的实部和虚部是单独的数组 Xre[i] 和 Xim[i]（在内存中是线性的，在分别对它们中的任何一个进行操作时效率很高）。

在FFTW中，实部和虚部默认交错为double[2]，即X[i][0]为实部，X[i][1]为虚部。

因此，要在 mex 文件中使用 FFTW 库，不能直接使用 Matlab 数组，而必须先分配新的内存，然后将来自 Matlab 的输入打包为 FFTW 格式，然后将来自 FFTW 的输出解包为 Matlab 格式。即

X = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
Y = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

然后

for (size_t i=0; i<N; ++i) {
    X[i][0] = Xre[i];
    X[i][1] = Xim[i];
}

然后

for (size_t i=0; i<N; ++i) {
    Yre[i] = Y[i][0];
    Yim[i] = Y[i][1];
}

因此，这需要 2 次内存分配 + 4 次读取 + 4 次写入 - 大小均为 N。这在处理大型问题时确实会影响速度。

我有一种预感，Mathworks 可能已经破解了 FFTW3 代码，使其能够直接以 Matlab 格式读取输入向量，从而避免了上述所有情况。

在这种情况下，只能分配 X 并将 X 用于 Y 以就地运行 FFTW（如 fftw_plan_*(N, X, X, ...) 而不是 fftw_plan_*(N, X, Y, ...)），因为它将被复制到 Yre 和 Yim Matlab 向量中，除非应用程序需要/受益于将 X 和 Y 分开。

编辑：实时查看运行Matlab的fft2()和我基于fftw3库的代码时的内存消耗，它表明Matlab只分配了一个额外的复杂数组（输出)，而我的代码需要两个这样的数组（*fftw_complex 缓冲区加上 Matlab 输出）。 Matlab 和 fftw 格式之间的就地转换是不可能的，因为 Matlab 的实数和虚数数组在内存中不是连续的。这表明 Mathworks 破解了 fftw3 库以使用 Matlab 格式读取/写入数据。

另一个对多次调用的优化是持久分配（使用mexMakeMemoryPersistent()）。我不确定 Matlab 实现是否也能做到这一点。

干杯。

附言作为旁注，Matlab 复数数据存储格式对于分别对实向量或虚向量进行操作更有效。在 FFTW 的格式中，您必须进行 ++2 内存读取。