如何在 MATLAB 中编写向量化函数

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

在我之前的问题中：MATLAB 中的嵌套 for 循环非常慢（预分配） https://stackoverflow.com/q/7811239/518169我意识到使用向量化函数与 4 个嵌套循环相比运行时间相差 7 倍.

在该示例中，无需循环遍历 4 维数组的所有维度并计算每个向量的中值，而是直接调用中位数（stack，n）（其中 n 表示中值函数的工作维度）更干净、更快。

但中位数只是一个非常简单的例子我只是幸运它有这个实施尺寸参数.

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

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

你如何写一个函数my_median_nD它的作用类似于 MATLAB 的中值，通过采用一个 n 维数组和一个“工作尺寸”范围？

Update

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

% 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 倍还有，那为什么这么复杂.

Update 2

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

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

Update 2 （解决您更新的问题）

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

这就是为什么MATLAB 内部做了some循环得很好- 通常将它们编码为编译函数。

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

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

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

原答案 （关于如何构建自己的在数组上工作的快速函数)

实际上有很多内置函数采用维度参数：min(stack, [], n), max(stack, [], n), mean(stack, n), std(stack, [], n), median(stack,n), sum(stack, n)...以及其他内置函数（例如exp(), sin()自动处理整个数组的每个元素（即sin(stack)自动为您执行四个嵌套循环，如果stack是 4D），您可以构建许多您可能需要的功能，只需依赖现有的内置函数.

如果这对于特定情况还不够，您应该看看repmat http://www.mathworks.de/help/techdoc/ref/repmat.html, bsxfun http://www.mathworks.de/help/techdoc/ref/bsxfun.html, arrayfun http://www.mathworks.de/help/techdoc/ref/arrayfun.html and accumarray http://www.mathworks.de/help/techdoc/ref/accumarray.html这些是非常强大的函数，可以“以 MATLAB 方式”执行操作。只需在 SO 上搜索问题（或者更确切地说是答案）using https://stackoverflow.com/search?q=%5Bmatlab%5D%20repmat one https://stackoverflow.com/search?q=%5Bmatlab%5D%20bsxfun of https://stackoverflow.com/search?q=%5Bmatlab%5D%20arrayfun these https://stackoverflow.com/search?q=%5Bmatlab%5D%20accumarray，通过这种方式，我学到了很多关于 MATLAB 的优点。

As an example，假设您想实施p-norm http://en.wikipedia.org/wiki/Norm_%28mathematics%29#p-norm沿维度的堆栈n，你可以写

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

...您可以有效地重用“哪个维度的能力”sum.

Update

正如马克斯在评论中指出的那样，还可以看看冒号运算符 (:) http://www.mathworks.de/help/techdoc/ref/colon.html这是一个非常强大的工具，用于从数组中选择元素（甚至改变它的形状，这通常是通过reshape http://www.mathworks.de/help/techdoc/ref/reshape.html).

一般来说，看看这个部分数组运算 http://www.mathworks.de/help/techdoc/ref/f16-5872.html#f16-55525在帮助中 - 它包含repmat等人。上面提到了，而且还cumsum以及一些更晦涩的辅助函数，您应该将其用作构建块。