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

Posted 2023-03-11

【中文标题】如何在 MATLAB 中编写向量化函数

【英文标题】：How to write vectorized functions in MATLAB

【发布时间】：2011-12-10 11:35:30

【问题描述】：

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

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

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

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

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

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

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

更新

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

% 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

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

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

【问题讨论】：

在 Matlab 命令提示符下键入“open median”，看看 Mathworks 是如何做到的！然而，他们作弊 - sort(X, dim) 是内置的。

【参考方案1】：

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

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

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

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

% 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() 等内置函数会自动作用于整个数组的每个元素（即，如果 stack 是 4D，sin(stack) 会自动为您执行四个嵌套循环），您可以构建很多您可能只需要依赖现有的内置函数。

如果这对于特定情况还不够，您应该查看repmat、bsxfun、arrayfun 和accumarray，它们是非常强大的“MATLAB 方式”处理函数。只需在 SO 上搜索问题（或更确切地说是答案）usingoneofthese，我通过这种方式学到了很多关于 MATLAB 的优势。

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

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

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

更新

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

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

【讨论】：

还要看看矩阵整形，以及 : 运算符的许多用途。

我用 'mean' 而不是 'median' 做了另一个测试来使用没有排序的函数，结果更加疯狂。这样，内置函数实际上快了 160 倍。这是 0.12 秒对 19 秒！感谢您的回答和更新！

【参考方案2】：

矢量化

除了已经说过的之外，您还应该了解vectorization 涉及并行化，即对数据执行并发操作而不是顺序执行（想想 SIMD 指令），甚至在某些情况下利用线程和多处理器...

MEX 文件

现在虽然“解释与编译”这一点已经争论过了，但没有人提到您可以通过编写 MEX 文件来扩展 MATLAB，这些文件是用 C 编写的编译可执行文件，可以从内部作为普通函数直接调用MATLAB。这允许您使用 C 等较低级别的语言来实现性能关键部分。

列主要顺序

最后，在尝试优化某些代码时，请始终记住 MATLAB 以列优先顺序存储矩阵。与其他任意顺序相比，按该顺序访问元素可以产生显着的改进。

例如，在您之前的链接问题中，您正在计算沿某个维度的堆叠图像集的 median。现在，这些维度的排序顺序极大地影响了性能。插图：

%# sequence of 10 images
fPath = fullfile(matlabroot,'toolbox','images','imdemos');
files = dir( fullfile(fPath,'AT3_1m4_*.tif') );
files = strcat(fPath,filesep,files.name');      %'

I = imread( files1 );

%# stacked images along the 1st dimension: [numImages H W RGB]
stack1 = zeros([numel(files) size(I) 3], class(I));
for i=1:numel(files)
    I = imread( filesi );
    stack1(i,:,:,:) = repmat(I, [1 1 3]);   %# grayscale to RGB
end

%# stacked images along the 4th dimension: [H W RGB numImages]
stack4 = permute(stack1, [2 3 4 1]);

%# compute median image from each of these two stacks
tic, m1 = squeeze( median(stack1,1) ); toc
tic, m4 = median(stack4,4); toc
isequal(m1,m4)

时间差很大：

Elapsed time is 0.257551 seconds.     %# stack1
Elapsed time is 17.405075 seconds.    %# stack4

【参考方案3】：

您能否向我解释一下为什么这段代码与简单的嵌套循环相比如此有效？它和其他函数一样有嵌套循环。

嵌套循环的问题不在于嵌套循环本身。这是您在内部执行的操作。

每个函数调用（尤其是对非内置函数的调用）都会产生一点点开销；如果函数执行例如，更是如此无论输入大小如何，错误检查都需要相同的时间。因此，如果一个函数只有 1 毫秒的开销，如果你调用它 1000 次，你将浪费一秒钟。如果您可以调用它一次来执行矢量化计算，那么您只需支付一次开销。

此外，JIT compiler (pdf) 可以帮助矢量化 simple for 循环，例如，您只执行基本的算术运算。因此，在您的帖子中进行简单计算的循环会加快很多，而调用 median 的循环则不会。

【参考方案4】：

在这种情况下

M = median(A,dim) returns the median values for elements along the dimension of A specified by scalar dim

但是使用通用函数，您可以尝试使用 mat2cell 拆分数组（它可以处理 n-D 数组，而不仅仅是矩阵）并通过 cellfun 应用您的 my_median_1D 函数。下面我将使用median 为例来说明你得到了相同的结果，但是你可以将它传递给m 文件中定义的任何函数，或者使用@(args) 符号定义的匿名函数。

>> testarr = [[1 2 3]' [4 5 6]']

testarr =

     1     4
     2     5
     3     6

>> median(testarr,2)

ans =

    2.5000
    3.5000
    4.5000

>> shape = size(testarr)

shape =

     3     2

>> cellfun(@median,mat2cell(testarr,repmat(1,1,shape(1)),[shape(2)]))

ans =

    2.5000
    3.5000
    4.5000

【讨论】：

（注意mat2cell 调用的输出是一个行向量元胞数组。）