MATLAB 速度优化

Question

有人能帮忙吗？我是一名经验丰富的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

- jack

Answer 1

一些需要考虑的事情：

你有没有考虑使用单值？

能否将cx、cy部分向量化，使其成为数组操作？

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

- Mikhail

---

单个变量可以减少内存占用。然而，基于这一点，我认为它会阻止优化：http://www.ee.columbia.edu/~marios/matlab/accel_matlab.pdf ，当然，进行测试不会有任何损失。 - Dennis Jaheruddin

@DennisJaheruddin，你错了。Matlab能够有效地将单精度指令转换为它们的SSE形式。我在许多场合都见证过这种效果。当发生这种情况时，你可以通过C++性能与Matlab相同来判断。我在那份文件中找不到有关数值精度的参考资料？ - Mikhail

谢谢Mikhail，将计算转换为单精度使计算速度提高了约10％。我觉得这很令人惊讶，因为我以前尝试过这样做，没有看到这样的收益。也许这与向多个核心传递更少的数据有关？ - jack

@Mikhail 我的评论是基于PDF第2页右下角的评论。也许我误解了这个评论，或者它在这里不适用/相关。我记得有些情况下执行单个计算比双精度计算慢。 - Dennis Jaheruddin

@DennisJaheruddin 哦，我明白了。无论如何，我是根据经验说话的，而这篇论文是2002年的。 - Mikhail

@jack 你应该看到一个更大的增长，可能是30%。检查一下所有的计算是否都转换成了单精度。杂乱的 0.5*myarray 会破坏它。复数指数是数学密集型的，很容易比内存访问时间更长。此外，减少内存占用可以帮助提高带宽。最后，许多GPU在双精度和单精度之间的性能差异达到20倍。 - Mikhail

Answer 2

如果您的数据是实数（不是复数），就像您的示例一样，您可以节省时间替换

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

by

(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))'

by

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

为了获得类似的收益。

- Luis Mendo

---

谢谢，虽然这是一个有趣的结果，但在这种情况下A和B是复杂的。我已经更新了问题以反映这一点。 - jack

@jack 请查看更新后的答案（水平线之后的第二部分）。 - Luis Mendo

Answer 3

有两种主要方法可以加速MATLAB代码：预分配和向量化。

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

向量化可以导致多个数量级的加速，并且会最优地使用核心（前提是标志已启用）。

你正在计算什么数学表达式？我可以帮忙吗？

- Luke

---

谢谢Luke。我已经尝试过向量化代码，以去除内部循环，但没有看到加速。最可能是由于需要操作的结果较大矩阵大小所致。但也有可能只是我的实现不好。该函数实际上是用于重建磁共振图像的反编码函数，与具有附加相位项的离散傅里叶变换非常相似（C（cy，cx）.* D）。 - jack

嗯...这是奇怪的行为。我很困惑，恐怕无法提供帮助。祝你好运找到解决方案！ - Luke

Answer 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

- StrongBad

---

谢谢Daniel的建议，但不幸的是使用这种方法会使单精度下的E达到4GB，超出了GPU内存的容量。通过为cx循环中的每个cx创建一个临时变量来减小E的大小可以解决这个内存问题。然而，这样做会花费太多时间在GPU和CPU之间传递信息，整个过程又会变慢。 - jack

@jack，你是如何提高你的GPU速度的？ - StrongBad

Answer 5

这一行代码：( x .* x_ind(cx) ) + ( y .* y_ind(cy) ) + ( C(cy,cx) .* D ); 是某种卷积，对吧？在频域中，循环卷积要快得多，并且通过使用FFT来优化频域的转换。

- GLK

Answer 6

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

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

- Yair Altman

Answer 7

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

此外，请尽可能向量化，一个简单的例子可以是：y.*y_ind(cy) 如果您只是为所有值创建正确的索引，则可以将其从最低的循环中删除。

- Dennis Jaheruddin

Answer 8

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

这样，您就不必在每次迭代中评估exp，而只需进行乘法运算，这应该会更快。

- bdecaf