【问题标题】：vectorize two nested for-loops in MATLAB在 MATLAB 中向量化两个嵌套的 for 循环
【发布时间】：2016-09-11 02:08:24
【问题描述】：

我有两个嵌套的 for 循环，用于格式化我加载的数据。循环具有以下结构：

data = magic(20000);
data = data(:,1:3);

for i=0:10
    for j=0:10
        data_tmp = data((1:100)+100*j+100*10*i,:);
        vx(:, i+1,j+1) = data_tmp(:,1);
        vy(:, i+1,j+1) = data_tmp(:,2);
        vz(:, i+1,j+1) = data_tmp(:,3);
    end
end

数组vx、vy 和vz 我会预先分配到它们所需的大小。但是，有没有办法对for-loops 进行矢量化以提高效率？由于第二个循环中的第一行data((1:100)+100*j+100*10*i,:)，我不相信这是这种情况，有没有更好的方法来做到这一点？

【问题讨论】：

标签： matlab for-loop vectorization

【解决方案1】：

原来你在循环中有重复的索引

i=k, j=10 和 i=k+1, j=0 k<10

例如，您阅读1:100 + 100*10 + 100*10*0，然后阅读1:100 + 100*0 + 100*10*1，它们是相同的。

使用重复索引重塑

如果这是您打算做的，那么矢量化需要多一步（索引生成）。

以下是我的建议（N=100, M=10 其中N 是data_tmp 的长度，M 是最大循环变量）

index = bsxfun(@plus,bsxfun(@plus,(1:N)',reshape(N*(0:M),1,1,M+1)),M*N*(0:M)); %index generation
vx = data(index);
vy = data(index + size(data,1));
vz = data(index + size(data,1)*2);

这不是那么理想，但它会起作用。

当我在笔记本电脑上进行测试时，它比使用预分配的原始代码快两倍。随着数据量的增加，差距越来越小。

无重复索引的重塑

如果不是，即，您想先在第 3 维方向上重塑每列，最后在第 2 维方向上），那么以下将起作用。

首先，这是我解释您的代码的方式

data = magic(20000);
data = data(:,1:3);

N = 100; M = 10;
for i=0:(M-1)
    for j=0:(M-1)
        data_tmp = data((1:N)+M*j+N*M*i,:);
        vx(:, i+1,j+1) = data_tmp(:,1);
        vy(:, i+1,j+1) = data_tmp(:,2);
        vz(:, i+1,j+1) = data_tmp(:,3);
    end
end

请注意，循环结束于 (M-1)。

以下是我的建议。

vx = permute(reshape(dat(1:N*M*M,1), N, M, M),[1,3,2]);
vy = permute(reshape(dat(1:N*M*M,2), N, M, M),[1,3,2]);
vz = permute(reshape(dat(1:N*M*M,3), N, M, M),[1,3,2]);

在我的笔记本电脑中，它比原始代码快 4 倍。随着尺寸的增加，差距接近 2。

【讨论】：

【解决方案2】：

以防万一你想坚持循环，这里有一个更快的方法：

data = randi(100,20000,3);
[vx,vy,vz] = deal(zeros(100,11,11));
[J,I] = ndgrid(1:11,1:11);
c = 1;
for k = 0:100:11000
    vx(:,I(c),J(c)) = data((1:100)+k,1);
    vy(:,I(c),J(c)) = data((1:100)+k,2);
    vz(:,I(c),J(c)) = data((1:100)+k,3);
    c = c+1;
end

我的猜测是 @Dohyun 答案中的 reshape 是您要寻找的（它比这快 x10，比您的代码快 x10000），但是对于下次使用循环时，这可能会很有用。

这是另一个不使用reshape 的选项，与reshape 版本类似：

[vx,vy,vz] = deal(zeros(100,10,11));
vx(:) = data(1:11000,1);
vy(:) = data(1:11000,2);
vz(:) = data(1:11000,3);
vx = permute(vx,[1 3 2]);
vy = permute(vy,[1 3 2]);
vz = permute(vz,[1 3 2]);

这个想法是在分配[vx,vy,vz] 时定义它们的形状。

【讨论】：