Matlab：加速大型数组运算 - 矢量化？

Question

我希望显着加快以下涉及大型数组（特别是名为“proj”的数组）操作的代码。与代码的其余部分相比，“proj”本身的计算足够快 - 请参阅“变量”部分。然而，“proj”的操作似乎需要花费很多时间。

下面的代码有三个版本，都在各自的部分中，标题为“原始代码，嵌套 for 循环”、“在 for 循环之前使 proj 1d”和“矢量化尝试，使 proj 4d”。所有三种方法似乎都提供相似的速度（非常慢）。尽管代码的某些部分不是最优的，但时间概况显示，所有三个版本的代码的瓶颈都与“proj”上的操作有关，该操作占用了大约 > 90% 的计算时间（计算温度和结果，第 38-39、61-62、80-81 行）。

我还尝试在一些外部循环上使用 parfor，但改进甚微。将来，我将增加变量 D、A、B、C 和 xyz 的大小，因此节省时间在这里至关重要。

%% variables

D = 5;
A = 360/D;
B = 0.2:0.5:5.2;
C = 100;
xyz = rand(500,3);

qx = cos(deg2rad(D).*(1:A));
qy = sin(deg2rad(D).*(1:A));

Rij = zeros(size(xyz,1),size(xyz,1),3);
for i1 = 1:size(Rij,3)
    for i2 = 1:size(Rij,2)
        for i3 = 1:size(Rij,1)
            Rij(i3,i2,i1) = xyz(i2,i1)-xyz(i3,i1);
        end
    end
end

proj = zeros(size(xyz,1),size(xyz,1),A);
Ccount = zeros(A,C);
for i1 = 1:A
    proj(:,:,i1) = Rij(:,:,1)*qx(i1) + Rij(:,:,2)*qy(i1);
    [Ccount(i1,1:C),~,~] = histcounts(proj(:,:,i1),C);
end

%% original code, nested for loops

tic;

result1 = zeros(A,length(B));

for i1 = 1:A
    for i2 = 1:length(B)
        temp = zeros(size(xyz,1),size(xyz,1));
        for i3 = 1:C
            temp = temp + (1./size(xyz,1)).*Ccount(i1,i3).*...
                cos(B(i2).*proj(:,:,i1));
        end
        temp = sum(temp,[1 2],'omitnan');
        result1(i1,i2) = temp;
    end
end
timer1 = toc;
disp(['time: ' num2str(round(timer1,1)) ' s']);

%% making proj 1D before for loops

tic;

proj1d = proj(:);
result2 = zeros(A,length(B));

for i1 = 1:A
    for i2 = 1:length(B)
        temp = zeros(size(xyz,1)*size(xyz,1),1);
        for i3 = 1:C
            i4 = i1*size(xyz,1).^2;
            i5 = i4 - size(xyz,1).^2 + 1;
            temp = temp + (1./size(xyz,1)).*Ccount(i1,i3).*...
                cos(B(i2).*proj1d(i5:i4));
        end
        temp = sum(temp,[1 2],'omitnan');
        result2(i1,i2) = temp;
    end
end
timer1 = toc;
disp(['time: ' num2str(round(timer1,1)) ' s']);

%% vectorization attempt (making proj 4d)

tic;

proj4d = repmat(proj,1,1,1,length(B));
B4d = repmat((permute(B,[1 3 4 2])),size(xyz,1),size(xyz,1),A);
temp = zeros(A,length(B));

for i3 = 1:C
    temp = temp + 1./size(xyz,1).*Ccount(:,i3).*...
        squeeze(sum(cos(proj4d.*B4d),[1 2],'omitnan'));
end
result3 = temp;
timer1 = toc;
disp(['time: ' num2str(round(timer1,1)) ' s']);

Answer 1

在原始的“嵌套 for 循环”代码中，内部循环一遍又一遍地重新计算相同的内容。不依赖于内部循环变量

i3

的计算应移出内部循环：

tic;
result1 = zeros(A,length(B));
for i1 = 1:A
    for i2 = 1:length(B)
        temp = zeros(size(xyz,1),size(xyz,1));
        for i3 = 1:C
            temp = temp + (1./size(xyz,1)) .* Ccount(i1,i3) .* ...
                cos(B(i2).*proj(:,:,i1));
        end
        temp = sum(temp,[1 2],'omitnan');
        result1(i1,i2) = temp;
    end
end
toc

变成：

tic;
result2 = zeros(A,length(B));
for i1 = 1:A
    for i2 = 1:length(B)
        temp = zeros(size(xyz,1),size(xyz,1));
        for i3 = 1:C
            temp = temp + (1./size(xyz,1)) .* Ccount(i1,i3);
        end
        temp = temp .* cos(B(i2).*proj(:,:,i1));
        temp = sum(temp,[1 2],'omitnan');
        result2(i1,i2) = temp;
    end
end
toc

这将我机器上的执行时间从 54.5 秒减少到 3.2 秒。

我们现在注意到，我们在

temp

中累积的是一个标量，只有将结果与

proj(:,:,i1)

组合时，它才变成矩阵。所以让我们保留

temp

为标量：

tic;
result2 = zeros(A,length(B));
for i1 = 1:A
    for i2 = 1:length(B)
        temp = 0;
        for i3 = 1:C
            temp = temp + (1./size(xyz,1)) .* Ccount(i1,i3);
        end
        temp = temp .* cos(B(i2).*proj(:,:,i1));
        temp = sum(temp,[1 2],'omitnan');
        result2(i1,i2) = temp;
    end
end
toc

现在需要 0.69 秒。

但是除以

size(xyz,1)

也与

i3

无关，我们也可以将其分解出来。所以内部循环只是对

Ccount(i1,:)

求和，我们不需要循环：

tic;
result2 = zeros(A,length(B));
for i1 = 1:A
    for i2 = 1:length(B)
        temp = sum(Ccount(i1,:)) ./ size(xyz,1);
        temp = temp .* cos(B(i2).*proj(:,:,i1));
        temp = sum(temp,[1 2],'omitnan');
        result2(i1,i2) = temp;
    end
end
toc

最后的更改使代码更简单，所以它很好，但实际上并没有显着减少计算时间，因为其他计算更昂贵。

接下来，我们再次注意到内循环的第一行不依赖于内循环变量，并且可以在外部完成：

tic;
result2 = zeros(A,length(B));
for i1 = 1:A
    x = sum(Ccount(i1,:)) ./ size(xyz,1);
    for i2 = 1:length(B)
        temp = x .* cos(B(i2).*proj(:,:,i1));
        result2(i1,i2) = sum(temp,[1 2],'omitnan');
    end
end
toc

我会把它留在这里。我们现在是 0.63 秒。

result2

与

result1

不同，因为我们改变了计算顺序，所以我们得到了不同的舍入误差。差异应该接近双浮点值的精度，

eps(result1)

:

max(abs((result1-result2)./eps(result1)),[],'all')

返回2。