performance - 如何在 MATLAB 中编写矢量化函数

Question

我刚刚开始学习 MATLAB，我发现很难理解循环与向量化函数的性能因素。

在我之前的问题中:Nested for loops extremely slow in MATLAB (preallocated)我意识到使用矢量化函数与 4 个嵌套循环相比，运行时间相差 7 倍。

在该示例中，与其循环遍历 4 维数组的所有维度并计算每个向量的中值，不如调用 median(stack, n) 更清晰、更快，其中 n 表示中值函数的工作维度。

但中位数只是一个非常简单的示例，我很幸运它实现了这个维度参数。

我的问题是，您如何自己编写一个与实现此维度范围的函数一样高效的函数？

例如，您有一个函数 my_median_1D，它仅适用于一维向量并返回一个数字。

如何通过采用 n 维数组和“工作维度” 参数编写函数my_median_nD，其作用类似于 MATLAB 的中位数？

更新

我找到了计算高维中位数的代码

% In all other cases, use linear indexing to determine exact location
% of medians.  Use linear indices to extract medians, then reshape at
% end to appropriate size.
cumSize = cumprod(s);
total = cumSize(end);            % Equivalent to NUMEL(x)
numMedians = total / nCompare;

numConseq = cumSize(dim - 1);    % Number of consecutive indices
increment = cumSize(dim);        % Gap between runs of indices
ixMedians = 1;

y = repmat(x(1),numMedians,1);   % Preallocate appropriate type

% Nested FOR loop tracks down medians by their indices.
for seqIndex = 1:increment:total
  for consIndex = half*numConseq:(half+1)*numConseq-1
    absIndex = seqIndex + consIndex;
    y(ixMedians) = x(absIndex);
    ixMedians = ixMedians + 1;
  end
end

% Average in second value if n is even
if 2*half == nCompare
  ixMedians = 1;
  for seqIndex = 1:increment:total
    for consIndex = (half-1)*numConseq:half*numConseq-1
      absIndex = seqIndex + consIndex;
      y(ixMedians) = meanof(x(absIndex),y(ixMedians));
      ixMedians = ixMedians + 1;
    end
  end
end

% Check last indices for NaN
ixMedians = 1;
for seqIndex = 1:increment:total
  for consIndex = (nCompare-1)*numConseq:nCompare*numConseq-1
    absIndex = seqIndex + consIndex;
    if isnan(x(absIndex))
      y(ixMedians) = NaN;
    end
    ixMedians = ixMedians + 1;
  end
end

你能给我解释一下为什么这段代码比简单的嵌套循环更有效吗？它与其他函数一样具有嵌套循环。

我不明白怎么会快 7 倍，而且为什么这么复杂。

更新 2

我意识到使用中位数并不是一个很好的例子，因为它本身就是一个复杂的函数，需要对数组进行排序或其他巧妙的技巧。我用 mean 重新做了测试，结果更加疯狂:19 秒对 0.12 秒。这意味着求和的内置方式比嵌套循环快 160 倍。

我真的很难理解一种行业领先的语言如何根据编程风格有如此极端的性能差异，但我看到了以下答案中提到的要点。

Answer 1

更新 2 (解决您更新后的问题)

MATLAB 经过优化，可以很好地处理数组。一旦您习惯了它，实际上只需键入一行并让 MATLAB 自行完成完整的 4D 循环工作而无需担心它，这实际上非常好。 MATLAB 通常用于原型(prototype)设计/一次性计算，因此为编码人员节省时间并放弃 C[++|#] 的一些灵 active 是有意义的。

这就是为什么 MATLAB 在内部某些 循环非常好 - 通常将它们编码为编译函数。

您提供的代码片段并不真正包含执行主要工作的相关代码行，即

% Sort along given dimension
x = sort(x,dim);

换句话说，您显示的代码只需要通过现在排序的多维数组 x 中的正确索引访问中值(这不会花费太多时间)。访问所有数组元素的实际工作是由 sort 完成的，它是一个内置(即编译和高度优化的)函数。

原始答案 (关于如何在数组上构建您自己的快速函数)

实际上有相当多的内置函数带有维度参数:min(stack, [], n), max(stack, [], n), mean(stack, n), std(stack, [], n), 中位数(stack,n), sum(stack, n)... 以及其他内置函数，如 exp()、sin() 自 Action 用于每个元素你的整个数组(即如果 stack 是 4D，sin(stack) 自动为你做四个嵌套循环)，你可以建立很多你的函数可能只需要依赖现有的内置组件。

如果这对于特定情况还不够，您应该查看 repmat , bsxfun , arrayfun和 accumarray这些功能非常强大，可以“以 MATLAB 方式”做事。只需在 SO 上搜索问题(或答案)using one of these , 我通过这种方式学到了很多关于 MATLAB 的优点。

作为示例，假设您想实现 p-norm沿维度 n 的堆栈，你可以写

function result=pnorm(stack, p, n)
result=sum(stack.^p,n)^(1/p);

...您可以有效地重用 sum 的“which-dimension-capability”。

更新

正如 Max 在评论中指出的那样，也可以查看 colon operator (:)这是一个非常强大的工具，用于从数组中选择元素(或者甚至改变它的形状，这通常用 reshape 完成)。

一般来说，看看Array Operations 部分在帮助中 - 它包含 repmat 等。上面提到的，还有 cumsum 和一些更晦涩的辅助函数，您应该将它们用作构建 block 。