'for' 循环与 MATLAB 中的矢量化答案

【问题标题】：'for' loop vs vectorization in MATLAB'for' 循环与 MATLAB 中的矢量化
【发布时间】：2011-11-26 00:43:05
【问题描述】：

我在 MATLAB 中进行了一些编程，并且按照建议，我一直在尝试使用矢量化。但最后程序很慢。所以我发现在一个地方使用循环时代码明显更快（下面的示例）。

我想知道我是否误解了什么或做错了什么，因为在这种情况下性能很重要，我不想继续猜测向量化或循环是否会更快。

% data initialization

k = 8;
n = 2^k+1;
h = 1/(n-1);
cw = 0.1;

iter = 10000;

uloc = zeros(n);
fploc = uloc;
uloc(2:end-1,2:end-1) = 1;
vloc = uloc;
ploc = ones(n);

uloc2 = zeros(n);
fploc2 = uloc2;
uloc2(2:end-1,2:end-1) = 1;
vloc2 = uloc2;
ploc2 = ones(n);

%%%%%%%%%%%%%%%%%%%%%%
% vectorized version %
%%%%%%%%%%%%%%%%%%%%%%
tic
for it=1:iter
    il=2:4;
    jl=2:4;
    fploc(il,jl) = h/6*(-uloc(il-1,jl-1) + uloc(il-1,jl)...
        -2*uloc(il,jl-1)+2*uloc(il,jl+1)...
        -uloc(il+1,jl) + uloc(il+1,jl+1)...
        ...
        -vloc(il-1,jl-1) - 2*vloc(il-1,jl)...
        +vloc(il,jl-1) - vloc(il,jl+1)...
        + 2*vloc(il+1,jl) + vloc(il+1,jl+1))...
        ...
        +cw*h^2*(-ploc(il-1,jl)-ploc(il,jl-1)+4*ploc(il,jl)...
        -ploc(il+1,jl)-ploc(il,jl+1));
end
toc


%%%%%%%%%%%%%%%%%%%%%%
%    loop version    %
%%%%%%%%%%%%%%%%%%%%%%
tic
for it=1:iter
    for il=2:4
        for jl=2:4
            fploc2(il,jl) = h/6*(-uloc2(il-1,jl-1) + uloc2(il-1,jl)...
                -2*uloc2(il,jl-1)+2*uloc2(il,jl+1)...
                -uloc2(il+1,jl) + uloc2(il+1,jl+1)...
                ...
                -vloc2(il-1,jl-1) - 2*vloc2(il-1,jl)...
                +vloc2(il,jl-1) - vloc2(il,jl+1)...
                + 2*vloc2(il+1,jl) + vloc2(il+1,jl+1))...
                ...
                +cw*h^2*(-ploc2(il-1,jl)-ploc2(il,jl-1)+4*ploc2(il,jl)...
                -ploc2(il+1,jl)-ploc2(il,jl+1));
        end
    end
end
toc

【问题讨论】：

标签： performance matlab for-loop vectorization

【解决方案1】：

我没有仔细阅读您的代码，但最新版本的 Matlab 中的 JIT 编译器已经改进到您所面临的情况很常见的程度 - 循环可能比矢量化代码更快。很难事先知道哪个会更快，因此最好的方法是以最自然的方式编写代码，对其进行分析，然后如果存在瓶颈，请尝试从循环切换到矢量化（或其他方式）。

【讨论】：

【解决方案2】：

MATLAB 的即时编译器 (JIT) 在过去几年中得到了显着改进。即使你是对的，通常应该对代码进行矢量化，但根据我的经验，这仅适用于某些操作和函数，并且还取决于你的函数正在处理多少数据。

找出最有效的方法是profile your MATLAB code，无论是否使用矢量化。

【讨论】：

当您说“这取决于您拥有多少数据”时，您介意更具体地说明您的意思吗？你的意思是循环通常在更大的数据集上表现更差？

【解决方案3】：

也许一个由几个元素组成的矩阵并不是一个很好的向量化效率测试。最后，它取决于应用程序的效果。

此外，通常矢量化代码看起来更好（更符合底层模型），但在很多情况下它不会，最终会损害实现。你所做的很棒，因为现在你知道什么最适合你。

【讨论】：

【解决方案4】：

我不会称之为矢量化。

您似乎正在执行某种过滤操作。这种过滤器的真正矢量化版本是原始数据乘以过滤器矩阵（即代表整个 for 循环的一个矩阵）。

这些矩阵的问题是它们非常稀疏（对角线周围只有几个非零元素），以至于使用它们几乎没有效率。您可以使用 sparse 命令，但即便如此，这种优雅的表示法可能并不能证明所需的额外内存是合理的。

Matlab 过去不擅长 for 循环，因为即使循环计数器等仍被视为复杂矩阵，因此在每次迭代时都会评估对此类矩阵的所有检查。我的猜测是，在你的 for 循环中，每次应用滤波器系数时，仍然会执行所有这些检查。

也许 matlab 函数 filter 和 filter2 在这里有用？您也可以阅读这篇文章：Improving MATLAB Matrix Construction Code : Or, code Vectorization for begginers

【讨论】：

【解决方案5】：

一种可能的解释是启动开销。如果在幕后创建了临时矩阵，请为内存分配做好准备。另外，我猜 MATLAB 不能推断出你的矩阵很小，所以会有循环开销。所以你的矢量化版本可能会以类似的代码结束

double* tmp=(double*)malloc(n*sizeof(double));
for(size_t k=0;k<N;++k)
    {
//  Do stuff with elements
    }
free(tmp);

将此与已知数量的操作进行比较：

double temp[2];
temp[0]=...;
temp[1]=...;

因此，与每次计算的工作负载相比，当 malloc-loopcounter-free 时间较长时，JIT 可能会更快。

【讨论】：