MATLAB中对向量进行索引是否效率低下？

Question

MATLAB中对向量进行索引是否效率低下？

arraysperformancematlabloopsvectorization

56

背景

我的问题是基于一些简单的观察而提出的，这些观察在一定程度上削弱了经验丰富的MATLAB用户通常持有/做出的信仰/假设：

MATLAB在内置函数和基本语言特性（如向量和矩阵索引）方面进行了非常好的优化。
在MATLAB中，循环速度较慢（尽管有JIT），如果算法可以用本机“向量化”的方式表达，则通常应避免使用循环。

结论：核心MATLAB功能高效，试图使用MATLAB代码超越它是难以实现的，如果不是不可能的话。

调查向量索引的性能

下面展示的示例代码是最基本的：我为所有向量条目分配标量值。首先，我分配一个空向量x：

tic; x = zeros(1e8,1); toc
Elapsed time is 0.260525 seconds.

如果我有一个向量x，我希望将所有元素都设为相同的值。在实际操作中，您可能会使用不同的方法，例如x = value*ones(1e8,1)，但是这里的重点是研究向量索引的性能。最简单的方法是编写：

tic; x(:) = 1; toc
Elapsed time is 0.094316 seconds.

我们称其为方法1（根据分配给x的值）。它似乎非常快（至少比内存分配快）。因为我这里所做的唯一一件事情就是操作内存，所以可以通过计算获得的有效内存带宽并将其与我的计算机的硬件内存带宽进行比较来估算此代码的效率：

eff_bandwidth = numel(x) * 8 bytes per double * 2 / time

在上述示例中，我乘以2，因为除非使用SSE流式处理，否则在内存中设置值需要同时从向量读取并写入到内存中。

eff_bandwidth(1) = 1e8*8*2/0.094316 = 17 Gb/s

我的电脑的STREAM基准内存带宽约为17.9 Gb/s，因此 - 在这种情况下MATLAB确实可以提供接近峰值性能！到目前为止，一切顺利。

如果您想将所有向量元素设置为某个值，则方法1很适合。但是，如果您想每隔step条目访问元素，则需要用例如1:step:end替换:。以下是与方法1的直接速度比较：

tic; x(1:end) = 2; toc
Elapsed time is 0.496476 seconds.

虽然你不会期望方法2表现出任何不同，但显然它造成了巨大的问题：没有理由的5倍减速。我怀疑在这种情况下MATLAB明确地分配了索引向量(1:end)。这在使用显式向量大小而不是end时得到了部分确认:

tic; x(1:1e8) = 3; toc
Elapsed time is 0.482083 seconds.

方法2和方法3的表现一样糟糕。

另一种可能性是显式地创建一个索引向量id并使用它来索引x。这样可以获得最灵活的索引能力。在我们的情况下：

tic;
id = 1:1e8; % colon(1,1e8);
x(id) = 4;
toc
Elapsed time is 1.208419 seconds.

这真是太厉害了 - 相比方法1，慢了12倍！我知道由于额外的内存使用（用于id），它应该比方法1表现更差，但为什么比方法2和3要差那么多呢？

让我们试着给循环一次机会 - 尽管听起来可能毫无希望。

tic;
for i=1:numel(x)
    x(i) = 5;
end
toc
Elapsed time is 0.788944 seconds.

一个大惊喜——循环比向量化方法快4倍，但仍比方法1、2和3慢。结果发现，在这种特殊情况下，你可以做得更好：

tic;
for i=1:1e8
    x(i) = 6;
end
toc
Elapsed time is 0.321246 seconds.

这项研究的最奇怪结果可能是：一个用MATLAB编写的循环明显优于本地向量索引。这肯定不应该是这样的。请注意，JIT编译的循环仍然比方法1几乎获得的理论峰值慢3倍。因此，仍有很大的改进空间。只是出乎意料（更适合的词语是“惊人”）的是通常的“向量化”索引（1:end）甚至更慢。

问题

在MATLAB中，简单索引（方法2、3和4比方法1更慢）非常低效，还是我漏掉了什么？
为什么方法4比方法2和3慢那么多？
为什么使用1e8作为循环边界而不是numel(x)可以将代码加速2倍？

编辑阅读了Jonas的评论后，这是另一种使用逻辑索引的方法：

tic;
id = logical(ones(1, 1e8));
x(id) = 7;
toc
Elapsed time is 0.613363 seconds.

比方法4好得多。

为方便起见：

function test

tic; x = zeros(1,1e8); toc

tic; x(:) = 1; toc
tic; x(1:end) = 2; toc
tic; x(1:1e8) = 3; toc

tic;
id = 1:1e8; % colon(1,1e8);
x(id) = 4;
toc

tic;
for i=1:numel(x)
    x(i) = 5;
end
toc

tic;
for i=1:1e8
    x(i) = 6;
end
toc

end

- angainor

6

继续研究Matlab的神秘方法吧 :) - Andrey Rubshtein

1

非常有趣。我给你的比较添加了两个：使用 length(x) 而不是 numel(x) 会慢得多。使用 x = x*6; 和方法一一样快。 - R. Schifini

1

@angainor：true(1e8,1) 可能比逻辑转换更快。 - Jonas

2

@angainor 很好的问题 +1。很遗憾我没有答案，但我想指出一些关于Matlab循环的有趣事情。尝试在每个解决方案周围放置 for z = 1:1。在我的机器上（R2012b），对于解决方案1到4没有影响，但是对于5和6会导致一个数量级的减速。换句话说，JIT加速器仍然无法发现双重循环中的低效率，即使外部循环只有一个迭代并且没有任何作用！哦，还有一件事，也许更改您的 i 下标。从技术上讲，i 是一个内置函数。出于同样的原因，请避免使用 j。 :-) - Colin T Bowers

3

@AndrewJanke 是的，这很有趣。我想这就是专有软件的问题——发现和诊断这些奇怪的问题唯一的方法是通过间接实验。Matlab语法优美，拥有出色的IDE，但整个封闭源代码的事情有时会让人烦恼... - Colin T Bowers

显示剩余7条评论

3个回答

9

我对所有问题都没有答案，但我对方法2、3和4有一些精细的推测。关于方法2和3，似乎MATLAB为向量索引分配内存，并用值从1到1e8填充它。为了理解这一点，让我们看看发生了什么。默认情况下，MATLAB使用double作为其数据类型。分配索引数组需要与分配x相同的时间。

tic; x = zeros(1e8,1); toc
Elapsed time is 0.260525 seconds.

目前，索引数组仅包含零。像方法1一样以最佳方式分配值给x向量需要0.094316秒。现在，必须从存储器中读取索引向量以便用于索引。这是额外的0.094316/2秒。请注意，在x(:)=1向量中，x必须从存储器中读取和写入。因此，仅读取它需要一半的时间。假设x(1:end)=value只执行这些操作，则方法2和3的总时间应为

t = 0.260525+0.094316+0.094316/2 = 0.402

几乎正确，但并不完全正确。我只能推测，填充索引向量的操作可能是一个额外的步骤，并需要额外的0.094316秒。因此，t=0.4963，这或多或少符合方法2和方法3的时间。

这些只是猜测，但它们似乎证实了当使用本机向量索引时，MATLAB会显式地创建索引向量。个人认为这是一种性能缺陷。MATLAB的JIT编译器应该足够聪明，能够理解这种琐碎的构造并将其转换为对正确内部函数的调用。就目前而言，在今天的内存带宽受限架构下，索引执行速度约为20%的理论峰值。

因此，如果您关心性能，您将不得不将x(1:step:end)实现为MEX函数，类似于：

set_value(x, 1, step, 1e8, value);

现在，在MATLAB中这显然是非法的，因为你不允许在MEX文件中就地修改数组。

编辑关于方法4，可以尝试按以下方式分析各个步骤的性能：

tic;
id = 1:1e8; % colon(1,1e8);
toc
tic
x(id) = 4;
toc

Elapsed time is 0.475243 seconds.
Elapsed time is 0.763450 seconds.

第一步，分配和填充索引向量的值所需时间与仅使用方法2和3相同。看起来这太多了 - 它应该最多需要分配内存和设置值的时间 (0.260525秒+0.094316秒=0.3548秒)，因此还有额外的开销约为0.12秒，我无法理解。第二部分 (x(id) = 4) 看起来也非常低效：它应该花费设置 x 值和读取 id 向量的时间 (0.094316秒+0.094316/2秒=0.1415秒) 加上对 id 值进行一些错误检查。用C编程，这两个步骤需要：

create id                              0.214259
x(id) = 4                              0.219768

代码使用检查一个double索引实际上代表一个整数，并且它适合x的大小:

tic();
id  = malloc(sizeof(double)*n);
for(i=0; i<n; i++) id[i] = i;
toc("create id");

tic();
for(i=0; i<n; i++) {
  long iid = (long)id[i];
  if(iid>=0 && iid<n && (double)iid==id[i]){
    x[iid] = 4;
  } else break;
}
toc("x(id) = 4");

第二步所需的时间比预期的0.1415秒要长，这是由于需要对id值进行错误检查。我认为开销过大-也许可以写得更好。尽管如此，所需的时间为0.4340秒，而不是1.208419秒。MATLAB在幕后做了什么-我不知道。也许有必要这样做，我只是看不到它。

当然，使用doubles作为索引引入了两个额外的开销级别：

double的大小是uint32的两倍。请记住，内存带宽是限制因素。

doubles需要转换为整数才能进行索引

可以使用整数索引在MATLAB中编写方法4：

tic;
id = uint32(1):1e8;
toc
tic
x(id) = 8;
toc

Elapsed time is 0.327704 seconds.
Elapsed time is 0.561121 seconds.

通过使用整数作为向量索引，性能明显提高了30％，这证明了我们应该使用整数作为向量的索引。然而，开销仍然存在。

从我目前的看法来看，在MATLAB框架内工作时，我们无法做任何事情来改进这种情况，我们必须等待Mathworks解决这些问题。

- angainor

@Jonas 我还没有这样做。我想在SO上问一下，看看我是否忽略了什么。我也希望有人从Mathworks回答 - 这里有很多来自MW的人。性能 bug 不是通常的 bug。代码可以工作并给出正确的结果。我预计 MW 的答案将是这是 Matlab 的设计特性，尽管我认为描述的一些问题至少可以通过改进 JIT 来消除。索引向量不需要显式创建。那只是一个简单的循环。或者它们不是，开销来自其他地方？ - angainor

TMW的人很少谈论实现细节，更不会在非正式场合谈及此事。不幸的是，自从他们开始了一个类似于Matlab Answers的SO克隆版后，他们似乎已经很少出现在SO上了。当您报告性能问题时，您可能会得到一些通用的答案，比如“我们正在处理它”，但只要他们将问题转发给开发人员，就有人知道我们知道还有工作要做。 - Jonas

@Jonas 我已经在 Matlab Central 上提交了这个问题。也许会引起一些工程师的兴趣。 - angainor

6

简单记录一下，在经过8年的发展中，MATLAB的性能特征已经发生了很大变化。

这是在R2017a（即OP帖子发布后5年）的情况：

Elapsed time is 0.000079 seconds.    % x = zeros(1,1e8);
Elapsed time is 0.101134 seconds.    % x(:) = 1;
Elapsed time is 0.578200 seconds.    % x(1:end) = 2;
Elapsed time is 0.569791 seconds.    % x(1:1e8) = 3;
Elapsed time is 1.602526 seconds.    % id = 1:1e8; x(id) = 4;
Elapsed time is 0.373966 seconds.    % for i=1:numel(x), x(i) = 5; end
Elapsed time is 0.374775 seconds.    % for i=1:1e8, x(i) = 6; end

注意循环使用1:numel(x)比索引x(1:end)更快，似乎数组1:end仍在创建中，而对于循环则不是这样。现在，在MATLAB中不要矢量化更好！

（我添加了一个赋值x(:)=0，在任何计时区域之外分配矩阵，以实际分配内存，因为zeros只保留内存。）

在MATLAB R2020b（在线）（3年后），我看到以下时间：

Elapsed time is 0.000073 seconds.    % x = zeros(1,1e8);
Elapsed time is 0.084847 seconds.    % x(:) = 1;
Elapsed time is 0.084643 seconds.    % x(1:end) = 2;
Elapsed time is 0.085319 seconds.    % x(1:1e8) = 3;
Elapsed time is 1.393964 seconds.    % id = 1:1e8; x(id) = 4;
Elapsed time is 0.168394 seconds.    % for i=1:numel(x), x(i) = 5; end
Elapsed time is 0.169830 seconds.    % for i=1:1e8, x(i) = 6; end

x(1:end) 现在被优化成与 x(:) 相同，向量 1:end 不再被明确创建。

- Cris Luengo

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Jonas · Accepted Answer

当然，我只能推测。但是当我在启用JIT编译器和禁用JIT编译器的情况下运行您的测试时，我得到了以下结果：

 % with JIT   no JIT
    0.1677    0.0011 %# init
    0.0974    0.0936 %# #1 I added an assigment before this line to avoid issues with deferring
    0.4005    0.4028 %# #2
    0.4047    0.4005 %# #3
    1.1160    1.1180 %# #4
    0.8221   48.3239 %# #5 This is where "don't use loops in Matlab" comes from 
    0.3232   48.2197 %# #6
    0.5464   %# logical indexing

分割可以展示我们在哪里有速度增加：

% withoutJit./withJit
    0.0067 %# w/o JIT, the memory allocation is deferred
    0.9614 %# no JIT
    1.0057 %# no JIT
    0.9897 %# no JIT
    1.0018 %# no JIT
   58.7792 %# numel
  149.2010 %# no numel

初始化时的明显加速是因为在关闭JIT时，MATLAB似乎会延迟内存分配直到使用它，所以x=zeros(...)实际上并没有做任何事情。(感谢@angainor)。

方法1到4似乎没有从JIT中受益。我想#4可能由于subsref中的额外输入测试而变慢，以确保输入的格式正确。 numel的结果可能与编译器难以处理不确定数量的迭代有关，或者与检查循环边界是否正确有一些开销（尽管没有JIT测试表明只需要约0.1秒）。

令人惊讶的是，在我的机器上运行R2012b时，逻辑索引似乎比#4更慢。

我认为这再次表明了MathWorks在加速代码方面做出的巨大工作，而“不要使用循环”并不总是最佳选择（至少目前是这样）。然而，我发现向量化通常是一个好方法，因为（a）JIT在更复杂的循环上失败了，（b）学习向量化可以让你更好地理解Matlab。

结论：如果您想要速度，请使用分析器，并在切换Matlab版本时重新进行分析。如评论中@Adriaan所指出的，现在使用timeit()来测量执行速度可能更好。

参考以下略微修改的测试函数

function tt = speedTest

tt = zeros(8,1);

tic; x = zeros(1,1e8); tt(1)=toc;

x(:) = 2;

tic; x(:) = 1; tt(2)=toc;
tic; x(1:end) = 2; tt(3)=toc;
tic; x(1:1e8) = 3; tt(4)=toc;

tic;
id = 1:1e8; % colon(1,1e8);
x(id) = 4;
tt(5)=toc;

tic;
for i=1:numel(x)
    x(i) = 5;
end
tt(6)=toc;

tic;
for i=1:1e8
    x(i) = 6;
end
tt(7)=toc;

%# logical indexing
tic;
id = true(1e8,1));
x(id)=7;
tt(8)=toc;