MATLAB 速度优化答案

【问题标题】：MATLAB Speed OptimisationMATLAB 速度优化
【发布时间】：2013-10-04 10:09:52
【问题描述】：

有人可以帮忙吗？我是一个相当有经验的 Matlab 用户，但在加速下面的代码时遇到了麻烦。

使用 12 个核心，我能够在所有三个循环中一次运行的最快时间约为 200 秒。实际函数将被调用约 720 次，按此速率执行将需要 40 多个小时。根据 Matlab 分析器，大部分 CPU 时间都花在了指数函数调用上。我已经设法使用 gpuArray 大大加快了速度，然后在 Quadro 4000 显卡上运行 exp 调用，但这会阻止使用 parfor 循环，因为工作站只有一个显卡，这会抹杀任何收益。任何人都可以提供帮助，或者这段代码是否接近使用 Matlab 可以实现的最佳值？我用 openMP 编写了一个非常粗略的 c++ 实现，但收效甚微。

在此先感谢

function SPEEDtest_CPU

% Variable setup:
% - For testing I'll use random variables. These will actually be fed into 
%   the function for the real version of this code.
sy    = 320;
sx    = 100;
sz    = 32;
A     = complex(rand(sy,sx,sz),rand(sy,sx,sz));
B     = complex(rand(sy,sx,sz),rand(sy,sx,sz));
C     = rand(sy,sx);
D     = rand(sy*sx,1);
F     = zeros(sy,sx,sz);
x     = rand(sy*sx,1);  
y     = rand(sy*sx,1);
x_ind = (1:sx) - (sx / 2) - 1;
y_ind = (1:sy) - (sy / 2) - 1;


% MAIN LOOPS 
%  - In the real code this set of three loops will be called ~720 times!
%  - Using 12 cores, the fastest I have managed is ~200 seconds for one
%    call of this function.
tic
for z = 1 : sz
    A_slice = A(:,:,z);
    A_slice = A_slice(:);
    parfor cx = 1 : sx       
        for cy = 1 : sy       
            E = ( x .* x_ind(cx) ) + ( y .* y_ind(cy) ) + ( C(cy,cx) .* D );                                                          

            F(cy,cx,z) = (B(cy,cx,z) .* exp(-1i .* E))' * A_slice; 
        end       
    end   
end
toc

end

【问题讨论】：

标签： performance matlab parallel-processing

【解决方案1】：

需要考虑的一些事情：

你考虑过单打吗？

您能否对 cx、cy 部分进行矢量化，以便它们表示数组操作？

考虑更改浮点舍入或信号模式。

【讨论】：

Single可以给你减少内存。但是，基于此我认为它会阻止优化：ee.columbia.edu/~marios/matlab/accel_matlab.pdf，当然测试不会受到伤害。
@DennisJaheruddin 你错了。 Matlab 能够有效地将单精度指令转换为它们的 SSE 形式。我在很多场合都目睹了这种效果。您可以知道这种情况何时发生，因为您的 C++ 性能与您的 Matlab 相同。我在该文档中找不到对数值精度的引用？
感谢 Mikhail，转换为单精度将计算速度提高了约 10%。我觉得这很令人惊讶，因为我之前尝试过这个并且没有看到这样的收获。也许这与将较少的数据传递给多个内核有关？
@Mikhail 我的评论基于 pdf 第 2 页右下角的评论。也许我误解了这一点，或者它在这里不适用/不相关。我确实记得在某些情况下执行单次计算比双次计算要慢。
@DennisJaheruddin 哦，我明白了。无论如何，我是根据经验说话的，这篇论文是 2002 年的。

【解决方案2】：

如果您的数据是真实的（不复杂），如您的示例所示，您可以节省替换时间

(B(cy,cx,z) .* exp(-1i .* E))'

通过

(B(cy,cx,z) .* (cos(E)+1i*sin(E))).'

具体来说，在我的机器上(cos(x)+1i*sin(x)).' 花费的时间比exp(-1i .* x)' 少19%。

如果A 和B 很复杂：E 仍然是真实的，因此您可以在循环之外预先计算Bconj = conj(B)（根据您的数据大小，这大约需要 10 毫秒，并且只完成一次）然后替换

(B(cy,cx,z) .* exp(-1i .* E))'

通过

(Bconj(cy,cx,z) .* (cos(E)+1i*sin(E))).'

获得类似的收益。

【讨论】：

谢谢，虽然结果很有趣，但在这种情况下 A 和 B 是复杂的。我已更新问题以反映这一点。
@jack 查看更新的答案（第二部分，水平线之后）

【解决方案3】：

加速 MATLAB 代码的主要方法有两种； 预分配和矢量化。

您已预先分配好，但没有矢量化。为了最好地学习如何做到这一点，您需要很好地掌握线性代数和使用repmat 将向量扩展为多个维度。

矢量化可以带来多个数量级的加速，并且会以最佳方式使用内核（假设标志处于上升状态）。

你计算的数学表达式是什么，我可以帮忙吗？

【讨论】：

谢谢卢克。我已经尝试过以这种方式对代码进行矢量化，以删除内部循环，但没有看到加速。很可能是由于必须操纵产生的大矩阵尺寸。但也许这只是代表我的一个糟糕的实施。该函数实际上是一个用于重建磁共振图像的逆编码函数，并且与添加了相位项 (C(cy,cx).* D) 的离散傅里叶变换非常相似。
嗯...这是奇怪的行为。我被难住了，我害怕。祝你好运找到解决方案！

【解决方案4】：

您可以将x .* x_ind(cx) 移出最内层循环。我没有方便的 GPU 来测试时序，但您可以将代码分成三个部分，以便您使用 GPU 和 parfor

for z = 1 : sz
    E = zeros(sy*sx,sx,sy);
    A_slice = A(:,:,z);
    A_slice = A_slice(:);
    parfor cx = 1 : sx
        temp = ( x .* x_ind(cx) );       
        for cy = 1 : sy       
            E(:, cx, cy) = temp + ( y .* y_ind(cy) ) + ( C(cy,cx) .* D );                                                          
        end
    end
    temp = zeros(zeros(sy*sx,sx,sy));
    for cx = 1 : sx
        for cy = 1 : sy       
             % Ideally use your GPU magic here
             temp(:, cx, cy) = exp(-1i .* E(:, cx, cy)));
        end
    end
    parfor cx = 1 : sx
        for cy = 1 : sy       
            F(cy,cx,z) = (B(cy,cx,z) .* temp(:, cx, cy)' * A_slice; 
        end       
    end   
end

【讨论】：

感谢 Daniel，但不幸的是，使用这种方法会使 E 单精度约为 4GB，对于 GPU 的内存来说太大了。通过为 cx 循环中的每个 cx 创建一个临时变量来减小 E 的大小，可以解决这个内存问题。然而，过多的时间花在了 GPU 和 CPU 之间来回传递信息上，整个事情又变慢了。
@jack 那么你如何提高 GPU 速度？

【解决方案5】：

为了实现正确的并行化，您需要确保循环完全独立，因此检查在每次运行中不分配给 E 是否有帮助。

进一步尝试尽可能多地矢量化，一个简单的例子可能是：y.*y_ind(cy)

如果您只是一次为所有值创建适当的索引，则可以将其从最低循环中取出。

【讨论】：

【解决方案6】：

不确定它是否对速度有很大帮助 - 但由于 E 基本上是一个总和，也许您可以使用可以预先计算的 exp (i cx(A+1)x) = exp(i cx(A) x) * exp(i x) 和 exp(i x)。

这样您就不必每次迭代都计算 exp - 只需乘以，这应该会更快。

【讨论】：

【解决方案7】：

除了其他人在这里给出的其他好建议之外，A_slice 的乘法与cx,cy 循环无关，并且可以在它们之外进行，一旦两个循环都完成后乘以F。

同样，B*exp(...) 的共轭也可以在 cx,cy 循环之外，在乘以 A_slice 之前批量完成。

【讨论】：

【解决方案8】：

这一行： ( x .* x_ind(cx) ) + ( y .* y_ind(cy) ) + ( C(cy,cx) .* D );

是某种类型的卷积，不是吗？循环卷积在频域中要快得多，并且使用 FTT 优化了与频域的转换。

【讨论】：