如何优化matlab中嵌套循环的运行速度

Question

我需要在matlab中处理四维矩阵的元素。其中

N

在100左右。下面的程序计算起来非常耗时，不知道有没有合适的方法可以简化一些？

for i1 = 1 : N 
    for j1 = 1 : N 

        t1 = conj(PF(i1, j1));

        if abs(t1) ~= 0
            for i2 = 1 : N 
                for j2 = 1 : N 
                
                    t2 = PF(i2, j2);
                    mumual_intensity(i2, j2, i1, j1) = 
mumual_intensity(i2, j2, i1, j1) * t1 * t2;
                end
            end
        end
    end
end

我知道可以将索引网格修改为向量计算，但我自己未能成功实现这一点。

Answer 1

首先，您应该始终对循环进行排序，以便第一个数组维度变化最快（即内部循环）。这就是数据在内存中的存储方式，按内存顺序访问数据可以充分利用缓存。因此，按如下方式重新排序循环应该会显着加快代码速度：

for j1 = 1 : N  % last dimension is outer loop
    for i1 = 1 : N  % second to last dimension

        t1 = conj(PF(i1, j1));

        if abs(t1) ~= 0
            for j2 = 1 : N  % second dimension
                for i2 = 1 : N  % first dimension is inner loop
                
                    t2 = PF(i2, j2);
                    mumual_intensity(i2, j2, i1, j1) = mumual_intensity(i2, j2, i1, j1) * t1 * t2;
                end
            end
        end
    end
end

假设

PF

是一个数组，我们可以轻松地向量化乘法。

.*

是逐元素乘法。通过在末尾添加两个维度并沿这些维度复制数据，

mumual_intensity .* PF

将自动将

t2

扩展为与

mumual_intensity

相同的大小。也就是说，每个

mumual_intensity(i2, j2, :, :)

都会乘以

PF(i2, j2)

。

对于其他乘法（通过

t1

），我们必须首先在开始处添加两个维度，我们可以使用

shiftdim(PF,-2)

来完成。因此，循环可以替换为：

mumual_intensity = mumual_intensity .* conj(shiftdim(PF,-2)) .* PF;

在装有 MATLAB R2023b 的 M1 Mac 上，使用

N=100

和随机双数组，我看到原始代码约为 0.3 秒，重新排序的循环约为 0.2 秒，矢量化代码约为 0.08 秒。差异不是很大，但对于足够大的

N

，任何增益都是好的。