很多嵌套循环 - 我怎样才能让它更快？

Question

我的这部分代码（如下所示）由于存在许多嵌套循环而需要很长时间才能运行。有没有办法可以消除这些嵌套循环以使其运行得更快？

for k = 1:numel(UpscaledZLen.pq)
    for j = 1:numel(UpscaledRowLen.pq)
        for i = 1:numel(UpscaledColLen.pq)

            iZ = 1;
            while iZ <= UpscaledZLen.pq(k)
                for ZLag = 1:iZ

                    ihRow = 1;
                    while ihRow <= UpscaledRowLen.pq(j)
                        for hRowLag = 1:ihRow

                            ihCol = 1;
                            while ihCol <= UpscaledColLen.pq(i)
                                for hColLag = 1:ihCol
                                    temp1(hColLag) = trapz(AnalGamma(hRowLag,...
                                        1:hColLag,...
                                        ZLag));
                                end
                                temp2(ihCol) = trapz(temp1);
                                ihCol = ihCol + 1;
                            end

                            temp3(hRowLag) = trapz(temp2);
                        end
                        temp4(ihRow) = trapz(temp3);
                        ihRow = ihRow + 1;
                    end

                    temp5(ZLag) = trapz(temp4);
                end
                temp6(iZ) = trapz(temp5);
                iZ = iZ + 1;
            end

            NormVariance_AnalCorrAvg(i, j, k) = (2/((UpscaledRowLen.pq(j)*...
                UpscaledColLen.pq(i)*UpscaledZLen.pq(k))^2))*trapz(temp6);
        end
    end
end

这段代码试图实现以下积分表达式：

EDIT (SSCCE): 举个简单的例子，我采用了以下变量大小，然后使用上面的代码查看运行这些特定变量大小需要多少时间： p>

nRow = 4;
nCol = 4;
nZ = 4;

RowLenScale.pq = 1:nRow;
ColLenScale.pq = 1:nCol;
UpscaledRowLen.pq = RowLenScale.pq(rem(RowLenScale.pq(end), RowLenScale.pq) == 0);
UpscaledColLen.pq = RowLenScale.pq(rem(ColLenScale.pq(end), ColLenScale.pq) == 0);

ZLenScale.pq = 1:nZ;
UpscaledZLen.pq = ZLenScale.pq(rem(ZLenScale.pq(end), ZLenScale.pq) == 0);

AnalGamma = rand(nRow, nCol, nZ);

对于这个例子，它只需要 0.321976 秒，而对于 nRow = 100;、nCol = 100; 和 nZ = 20; 的原始情况，它需要超过 24 小时并且仍在使用最外层循环索引运行在k = 4。

Answer 1

首先，您的代码有问题，因为临时变量没有初始化，并且上次迭代留下的值在下一次迭代中仍然存在。更正此问题并将while 替换为for 循环后，代码如下所示：

NormVariance_AnalCorrAvg = nan(numel(UpscaledColLen.pq), numel(UpscaledRowLen.pq), numel(UpscaledZLen.pq));
for k = 1 : numel(UpscaledZLen.pq)
    for j = 1 : numel(UpscaledRowLen.pq)
        for i = 1 : numel(UpscaledColLen.pq)

            temp6 = nan(1, UpscaledZLen.pq(k));
            for iZ = 1 : UpscaledZLen.pq(k)
                temp5 = nan(1, iZ);
                for ZLag = 1 : iZ
                    temp4 = nan(1, UpscaledRowLen.pq(j));
                    for ihRow = 1 : UpscaledRowLen.pq(j)
                        temp3 = nan(1, ihRow);
                        for hRowLag = 1 : ihRow
                            temp2 = nan(1, UpscaledColLen.pq(i));
                            for ihCol = 1 : UpscaledColLen.pq(i)
                                temp1 = nan(1, ihCol);
                                for hColLag = 1:ihCol
                                    temp1(hColLag) = trapz(AnalGamma(hRowLag,...
                                        1:hColLag,...
                                        ZLag));
                                end
                                temp2(ihCol) = trapz(temp1);
                            end
                            temp3(hRowLag) = trapz(temp2);
                        end
                        temp4(ihRow) = trapz(temp3);
                    end
                    temp5(ZLag) = trapz(temp4);
                end
                temp6(iZ) = trapz(temp5);
            end
            NormVariance_AnalCorrAvg(i, j, k) = (2/((UpscaledRowLen.pq(j)*...
                UpscaledColLen.pq(i)*UpscaledZLen.pq(k))^2))*trapz(temp6);

        end
    end
end

这里的初始化也会导致数组的预分配，这是加速循环的基本推荐方法之一。但是，在这种情况下，它没有效果。

为了对代码进行基准测试，我使用rand 使用您的 sn-p，但使用了

nRow = 20;
nCol = 20;
nZ = 2;

在我的电脑上，清理后的代码需要 25.5 秒才能完成。

---

如何让它更快？好旧的矢量化：

hColLag 上的最内层循环可以替换为 cumtrapz，将运行时间减少到 8.2 秒。

temp1 = cumtrapz(AnalGamma(hRowLag, 1 : ihCol, ZLag), 2);

ihCol 上的循环也是如此，将运行时间减少到 1.85 秒。

temp2 = cumtrapz(cumtrapz(AnalGamma(hRowLag, 1 : UpscaledColLen.pq(i), ZLag), 2), 2);

hRowLag 上的下一个外循环仅用于计算几个类似的trapzs。这不是必需的，因为trapz 是完全矢量化的；循环可以被一个调用替换，这将运行时间减少到 0.34 秒。

temp3 = trapz(cumtrapz(cumtrapz(AnalGamma(1 : ihRow, 1 : UpscaledColLen.pq(i), ZLag), 2), 2), 2);

ihRow 上的循环计算累积积分；使用cumtrapz 可以将运行时间缩短到大约 0.07 秒。

temp4 = cumtrapz(trapz(cumtrapz(cumtrapz(AnalGamma(1 : UpscaledRowLen.pq(j), 1 : UpscaledColLen.pq(i), ZLag), 2), 2), 2), 1);

对iZ 和ZLag 上的循环应用相同的逻辑，我们得到以下代码

NormVariance_AnalCorrAvg = nan(numel(UpscaledColLen.pq), numel(UpscaledRowLen.pq), numel(UpscaledZLen.pq));
for k = 1 : numel(UpscaledZLen.pq)
    for j = 1 : numel(UpscaledRowLen.pq)
        for i = 1 : numel(UpscaledColLen.pq)
            temp6 = cumtrapz(trapz(cumtrapz(trapz(cumtrapz(cumtrapz(...
                AnalGamma(1 : UpscaledRowLen.pq(j), 1 : UpscaledColLen.pq(i), 1 : UpscaledZLen.pq(k)), ...
                2), 2), 2), 1), 1), 3);
            NormVariance_AnalCorrAvg(i, j, k) = (2/((UpscaledRowLen.pq(j)*...
                UpscaledColLen.pq(i)*UpscaledZLen.pq(k))^2))*trapz(temp6);

        end
    end
end

运行时间为 0.043 秒，提高了约 600 倍。

编写所有这些嵌套的trapzs 和cumtrapzs 当然很容易出错。我使用了应用于保存到文件的结果的校验和，以确保它不会因我对代码的更改而改变。不过，您应该仔细检查以完全确定。

---

也许可以走得更远。前三个 for 循环再次做一些累积的事情，所以可能更多的cumtrapz 是有序的。 trapz 没有很有效地实现，因为它有很多错误检查代码。由于它所做的基本上是平均值的总和（请参阅trapezoidal rule），因此在整个网格上预先计算多维平均值将允许将这些调用转换为对sum 的调用。比梯形规则更简单的是rectangle method，您将计算函数的值，而不是在边缘，而是在网格单元的中心。然后当然你可以尝试分析部分地计算积分；并且您可以降低网格的分辨率。

希望这会有所帮助！

Answer 2

您可以采取一些措施来提高性能。首先，您可以用 for 循环替换 while 循环。其次，您可以用 sum 替换集成的 trapz。如果你用 sum 代替，而不是在每个循环中计算总和，而是在最后计算总和，你会得到一些改进。最后，您可以提取最终的乘法计算。

伽马函数的输入真的是（hRowLag，1：hColLag，ZLag）还是（hRowLag，hColLag，ZLag）。如果 gamma 函数可以矢量化，您可以构建一个网格来为函数提供数据，而不是循环。

temp = zeros(max(UpscaledZLen.pq),max(UpscaledZLen.pq),max(UpscaledRowLen.pq),max(UpscaledRowLen.pq),max(UpscaledColLen.pq),max(UpscaledColLen.pq));
my_integral = zeros(numel(UpscaledZLen.pq), numel(UpscaledRowLen.pq), numel(UpscaledColLen.pq));

for k = 1:numel(UpscaledZLen.pq)
    for j = 1:numel(UpscaledRowLen.pq)
        for i = 1:numel(UpscaledColLen.pq)
            for iZ = 1:UpscaledZLen.pq(k)
                for ZLag = 1:iZ
                    for ihRow = 1:UpscaledRowLen.pq(j)
                        for RLag = 1:ihRow
                            for ihCol = 1:UpscaledColLen.pq(i)
                                for CLag = 1:ihCol
                                    temp(iZ,ZLag,ihRow,RLag,ihCol,CLag,:) = AnalGamma(RLag,1:CLag,ZLag);
                                end
                            end
                        end
                    end
                end
            end
            my_integral(i,j,k) = sum(temp(:));
        end
    end
end

my_z = repmat(permute(UpscaledZLen.pq(:),[3,2,1]),[numel(UpscaledRowLen.pq),numel(UpscaledColLen.pq),1]);
my_r = repmat(permute(UpscaledRowLen.pq(:),[2,1]),[numel(UpscaledRowLen.pq),1,numel(UpscaledZLen.pq)]);
my_c = repmat(UpscaledColLen.pq(:),[1,numel(UpscaledColLen.pq),numel(UpscaledZLen.pq)]);

NormVariance_AnalCorrAvg = my_integral .* 2 ./ (my_z .* my_r .* my_c).^2;