arrayfun 可能比 matlab 中的显式循环慢得多。为啥？

Posted 2023-02-16

【中文标题】arrayfun 可能比 matlab 中的显式循环慢得多。为啥？

【英文标题】：arrayfun can be significantly slower than an explicit loop in matlab. Why?arrayfun 可能比 matlab 中的显式循环慢得多。为什么？

【发布时间】：2012-09-13 10:12:17

【问题描述】：

考虑以下arrayfun 的简单速度测试：

T = 4000;
N = 500;
x = randn(T, N);
Func1 = @(a) (3*a^2 + 2*a - 1);

tic
Soln1 = ones(T, N);
for t = 1:T
    for n = 1:N
        Soln1(t, n) = Func1(x(t, n));
    end
end
toc

tic
Soln2 = arrayfun(Func1, x);
toc

在我的机器上（Linux Mint 12 上的 Matlab 2011b），这个测试的输出是：

Elapsed time is 1.020689 seconds.
Elapsed time is 9.248388 seconds.

什么？！？ arrayfun 虽然公认是一个看起来更干净的解决方案，但速度要慢一个数量级。这是怎么回事？

此外，我对cellfun 进行了类似的测试，发现它比显式循环慢了大约 3 倍。同样，这个结果与我的预期相反。

我的问题是： 为什么arrayfun 和cellfun 这么慢？鉴于此，是否有任何充分的理由使用它们（除了让代码看起来不错）？

注意：我在这里说的是arrayfun 的标准版本，而不是并行处理工具箱中的GPU 版本。

编辑：为了清楚起见，我知道上面的 Func1 可以像 Oli 指出的那样被矢量化。我之所以选择它，是因为它会针对实际问题进行简单的速度测试。

编辑：按照 grunetta 的建议，我使用 feature accel off 重新进行了测试。结果是：

Elapsed time is 28.183422 seconds.
Elapsed time is 23.525251 seconds.

换句话说，很大一部分差异似乎在于 JIT 加速器在加速显式 for 循环方面比 arrayfun 做得更好。这对我来说似乎很奇怪，因为arrayfun 实际上提供了更多信息，即它的使用表明对Func1 的调用顺序无关紧要。另外，我注意到无论 JIT 加速器是打开还是关闭，我的系统都只使用一个 CPU...

【问题讨论】：

幸运的是，“标准解决方案”仍然是迄今为止最快的：tic; 3*x.^2+2*x-1; toc 经过的时间是 0.030662 秒。

@Oli 我想我应该预料到有人会指出这一点并使用无法矢量化的函数:-)

我很想看看当 JIT 加速器关闭时这个时间是如何变化的。执行命令“feature accel off”，然后重新运行您的测试。

@grungetta 有趣的建议。我已将结果与一些 cmets 添加到问题中。

让我把这个添加到相关问题列表中：What is the fastest way to perform arithmetic operations on each element of a cell array?

【参考方案1】：

您可以通过运行其他版本的代码来获得这个想法。考虑明确写出计算，而不是在循环中使用函数

tic
Soln3 = ones(T, N);
for t = 1:T
    for n = 1:N
        Soln3(t, n) = 3*x(t, n)^2 + 2*x(t, n) - 1;
    end
end
toc

在我的电脑上计算的时间：

Soln1  1.158446 seconds.
Soln2  10.392475 seconds.
Soln3  0.239023 seconds.
Oli    0.010672 seconds.

现在，虽然完全“矢量化”的解决方案显然是最快的，但您可以看到为每个 x 条目定义一个要调用的函数是一个巨大开销。只是明确地写出计算让我们加速了 5 倍。我想这表明 MATLABs JIT 编译器does not support inline functions。根据那里的gnovice的回答，实际上最好写一个普通函数而不是匿名函数。试试看。

下一步 - 移除（矢量化）内部循环：

tic
Soln4 = ones(T, N);
for t = 1:T
    Soln4(t, :) = 3*x(t, :).^2 + 2*x(t, :) - 1;
end
toc

Soln4  0.053926 seconds.

另一个因素 5 加速：这些语句中有一些内容说您应该避免 MATLAB 中的循环......或者真的有吗？那就看看这个吧

tic
Soln5 = ones(T, N);
for n = 1:N
    Soln5(:, n) = 3*x(:, n).^2 + 2*x(:, n) - 1;
end
toc

Soln5   0.013875 seconds.

更接近“完全”矢量化版本。 Matlab 按列存储矩阵。您应该始终（在可能的情况下）将您的计算构造为“按列”矢量化。

我们现在可以回到Soln3了。那里的循环顺序是“逐行”的。让我们改变它

tic
Soln6 = ones(T, N);
for n = 1:N
    for t = 1:T
        Soln6(t, n) = 3*x(t, n)^2 + 2*x(t, n) - 1;
    end
end
toc

Soln6  0.201661 seconds.

更好，但仍然很糟糕。单循环 - 好。双循环 - 不好。我猜 MATLAB 在提高循环性能方面做了一些不错的工作，但循环开销仍然存在。如果你在里面有一些更重的工作，你不会注意到的。但由于此计算受内存带宽限制，您确实会看到循环开销。您会更清楚地看到在那里调用 Func1 的开销。

那么，arrayfun 是怎么回事？那里也没有内联函数，所以开销很大。但是为什么比双嵌套循环差那么多呢？实际上，使用 cellfun/arrayfun 的话题已经被广泛讨论过很多次（例如here、here、here 和here）。这些函数很慢，您不能将它们用于如此细粒度的计算。您可以使用它们来简化代码以及在单元格和数组之间进行花哨的转换。但是函数需要比你写的更重：

tic
Soln7 = arrayfun(@(a)(3*x(:,a).^2 + 2*x(:,a) - 1), 1:N, 'UniformOutput', false);
toc

Soln7  0.016786 seconds.

请注意，Soln7 现在是一个单元格.. 有时这很有用。代码性能现在已经相当不错了，如果需要cell作为输出，使用全向量化解后就不需要再转换矩阵了。

那么为什么 arrayfun 比简单的循环结构慢呢？不幸的是，我们不可能肯定地说，因为没有可用的源代码。你只能猜测，由于arrayfun是一个通用的函数，它处理各种不同的数据结构和参数，在简单的情况下不一定很快，你可以直接表示为循环嵌套。我们无法知道开销从何而来。是否可以通过更好的实现来避免开销？也许不吧。但不幸的是，我们唯一能做的就是研究性能，以确定哪些情况下效果很好，哪些情况下效果不佳。

更新由于这个测试的执行时间很短，为了得到可靠的结果，我现在添加了一个围绕测试的循环：

for i=1:1000
   % compute
end

下面给出了一些时间：

Soln5   8.192912 seconds.
Soln7  13.419675 seconds.
Oli     8.089113 seconds.

您看到，arrayfun 仍然很糟糕，但至少不比矢量化解决方案差三个数量级。另一方面，具有按列计算的单个循环与完全矢量化版本一样快......这一切都是在单个 CPU 上完成的。如果我切换到 2 个内核，Soln5 和 Soln7 的结果不会改变 - 在 Soln5 中，我必须使用 parfor 来使其并行化。忘记加速... Soln7 不并行运行，因为 arrayfun 不并行运行。另一方面，Olis 矢量化版本：

Oli  5.508085 seconds.

【讨论】：

很好的答案！和 matlab central 的链接都提供了非常有趣的阅读。非常感谢。

还有一个有趣的更新！这个答案一直在给出:-)

只是一个小评论；回到 MATLAB 6.5，cellfun 被实现为 MEX 文件（旁边有 C 源代码）。这实际上很简单。当然它只支持应用 6 个硬编码函数之一（你不能传递函数句柄，只能传递一个带有函数名称的字符串）

arrayfun + 函数句柄 = 慢！在繁重的代码中避免它们。

@Amro 如果 arrayfun 和 cellfun 如你所说的那样实现，那么很可能 MATLAB 不使用任何 SIMD 或 MIMD 命令。但它可能是，JIT 编译器知道如何做，就像每个好的编译器一样。这可以解释与“Soln5”和“Oli”的时差因素。

【参考方案2】：

那是因为！！！！

x = randn(T, N); 

不是gpuarray类型；

你需要做的就是

x = randn(T, N,'gpuArray');

【讨论】：

我认为您需要更仔细地阅读@angainor 的问题和出色的答案。它与gpuarray 没有任何关系。几乎可以肯定，这就是为什么这个答案被否决的原因。

@Colin - 我同意 angainor 的更彻底，但答案没有提到“gpuArray”。我认为'gpuArray'在这里是一个很好的贡献（如果它是正确的）。此外，“这里发生了什么？” 这个问题有点草率，所以我认为它为其他方法打开了大门，例如矢量化数据并将其发送到 GPU。我让这个答案骑起来，因为它可以为未来的访问者增加价值。如果我打错电话了，我深表歉意。

您还忘记了gpuarray 仅支持 nVidia 显卡这一事实。如果他们没有这样的硬件，那么您的建议（或缺乏）就毫无意义。 -1

另一方面，gpuarray是matlab矢量化编程的光剑。