当我无法对计算进行矢量化时，如何加快 Matlab 嵌套的 for 循环？

Question

我有三个大小相同的大 3D 数组 [41*141*12403]，在 Matlab 代码中命名为 alpha、beta 和 ni。从它们中，我需要计算另一个具有相同大小的 3D 数组，该数组是通过使用每个元素的值结合无限和和定积分计算的计算从原始矩阵中逐元素获得的。因此，必须使用多个嵌套循环来进行此计算似乎是不可避免的。代码现在已经运行了几个小时（！），它仍然处于外循环的第一次迭代中（需要执行 41 次！！根据我的计算，这样程序将必须运行超过两年！！！）。我不知道如何优化代码。请帮我 ！！

我使用的代码：

    z_len=size(KELDYSH_PARAM_r_z_t,1);   % 41 rows
    r_len=size(KELDYSH_PARAM_r_z_t,2);   % 141 columns   
    t_len=size(KELDYSH_PARAM_r_z_t,3);   % 12403 slices

    sumRes=zeros(z_len,r_len,t_len);

    for z_ind=1:z_len
        z_ind     % in order to track the advancement of the calculation
        for r_ind=1:r_len
            for t_ind=1:t_len
                sumCurrent=0;
                sumPrevious=inf;
                s=0;

                while abs(sumPrevious-sumCurrent)>1e-6
                    kapa=kapa_0+s;    %some scalar
                    x_of_w=(beta(z_ind,r_ind,t_ind).*(kapa-ni...
                       (z_ind,r_ind,t_ind))).^0.5;               
                    sumPrevious=sumCurrent;
                    sumCurrent=sumCurrent+exp(-alpha(z_ind,r_ind,t_ind).* ...
                        (kapa-ni(z_ind,r_ind,t_ind))).*(x_of_w.^(2*abs(m)+1)/2).* ...
                            w_m_integral(x_of_w,m);
                    s=s+1;
                end

                sumRes(z_ind,r_ind,t_ind)=sumCurrent;
            end
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function  res=w_m_integral(x_of_w,m)

    res=quad(@integrandFun,0,1,1e-6);

    function y=integrandFun(t)
            y=exp(-x_of_w^2*t).*t.^(abs(m))./((1-t).^0.5);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Answer 1

选项 1 - 更多矢量化

您正在使用的模型非常复杂，并未解释所有术语，但某些部分仍可进一步矢量化。您的alpha、beta 和ni 矩阵大概是静态的和预先计算的？您的 s 值是一个标量，kapa 可以是任何一个，因此您也可以一次性预先计算 x_of_w 矩阵。这将给你一个非常小的加速，虽然你会花费内存来获得它——这些天7100万点是可行的，但需要大量的硬件。对你的 41 行中的每一行执行一次，可以巧妙地减轻负担。

剩下的就是积分本身。 quad 函数不接受向量输入 - 这将是一场噩梦，不是吗？ - integral 也没有，Mathworks 建议您改用它。但是，如果您的积分限制在每种情况下都相同，那么为什么不采用老式的方式进行积分呢？计算被积函数在 1 处的值的矩阵，计算被积函数在 0 处的值的另一个矩阵，然后取差。

然后您可以编写一个循环来计算整个输入空间的积分，然后测试所有矩阵元素的收敛性。制作一个掩码，记下那些没有收敛的掩码，并重新计算那些增加了s 的掩码。重复直到所有都收敛（或者你达到了迭代的阈值）。

选项 2 - 并行化

过去，matlab 的向量化操作比循环快得多。我现在找不到它的来源，但我想我已经读到它最近使用for 循环变得更快了，所以根据你可用的资源，你可能会通过并行化你当前的代码来获得更好的结果有。这也需要进行一些重构 - 最大的问题是在将数据复制到工作人员时的开销（您可以通过将输入分成块并只提供相关的输入来解决）和 parfor 循环不允许您可以使用某些变量，通常是覆盖整个空间的变量。再次切碎它们会有所帮助。

但是，如果您有 2 年的运行时间，我猜您至少需要 100 倍，这意味着集群！如果您在大学或其他地方，您可能能够在 500 核集群上获得几天时间，那就去吧...

如果您可以将积分写成封闭形式，那么它可能适合 GPU 计算。这些东西可以非常快速地完成某些类别的计算，但您必须能够并行化工作并将实际计算减少到主要包括加法和乘法的基本计算。 CUDA 库已经做了很多工作，matlab has an interface to them 所以请阅读这些内容。

选项 3 - 缩小范围

最后，如果上述两种方法都不能带来足够的加速，那么您可能必须缩小计算范围。尽可能多地修剪输入空间，并可能接受较低的收敛阈值。如果您知道在最里面的while 循环（其中有s 计数器的循环）内倾向于需要多少次迭代，那么可能会发现减少收敛标准会减少您需要的迭代次数，这可以加快它了。分析器可以帮助您了解您的时间都花在了哪些地方。

但底线是 7100 万个点需要一些时间来计算。到目前为止，您只能优化计算，对于这种规模的问题，您可能不得不投入硬件。