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

6

我正在学习MATLAB,发现难以理解循环和矢量化函数的性能因素

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

在这个例子中,不是遍历4维数组的所有维度并为每个向量计算中位数,而是用median(stack, n)调用方式更加简洁快速,其中n表示median函数的工作维度。

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

我的问题是如何编写一个函数,使其像一个具有此维度范围实现的函数一样高效工作

例如,您有一个仅针对1-D向量并返回数字的函数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秒。这意味着内置求和方法比嵌套循环快160倍

对于我来说,真的很难理解为什么行业领先的语言在编程风格上具有如此极端的性能差异,但我看到了以下答案中提到的观点。


在 Matlab 命令提示符处键入“open median”并查看 Mathworks 如何执行!不过他们作弊了 - sort(X,dim) 是内置的。 - Max
4个回答

6

更新2(针对您更新的问题)

MATLAB针对数组进行了优化,因此一旦您习惯了它,只需输入一行代码,MATLAB就可以自己完成完整的4D循环操作,无需担心。MATLAB通常用于原型设计/一次性计算,因此节省编码者的时间并放弃C[++|#]的某些灵活性是有意义的。

这就是为什么MATLAB在内部实现了某些循环操作 - 通常通过将其编码为编译函数来实现。

您提供的代码片段实际上并不包含执行主要工作的相关代码行:

% 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)自动为您执行四个嵌套循环),您可以依靠现有的内置函数构建出许多可能需要的函数。
如果这些对于某个特定情况不够用的话,您可以查看 repmatbsxfunarrayfunaccumarray 这些非常强大的函数,以“MATLAB方式”完成任务。只需在SO上搜索问题(或者更确切地说是答案) 使用 其中之一 解决,我通过这种方法学到了很多关于MATLAB的优点。
作为一个例子,假设你想要在维度 n 上实现堆栈的 p-norm,你可以这样写。
function result=pnorm(stack, p, n)
result=sum(stack.^p,n)^(1/p);

...在这里,您可以有效地重复使用sum的“哪个维度能力”。

更新

正如Max在评论中指出的那样,还要看一下冒号操作符(:),它是从数组中选择元素(甚至可以改变其形状,这通常是用reshape来完成的)的非常强大的工具。

总的来说,请查看帮助中的数组操作部分 - 它包含了上面提到的repmat等函数,以及cumsum和一些更为晦涩的辅助函数,您应该将它们用作构建块。


1
还要看矩阵重塑和 : 运算符的许多用途。 - Max
我进行了另一项测试,使用“均值”而不是“中位数”来使用一个无需排序的函数,结果更加疯狂。这种方式内置函数实际上快了160倍。这是0.12秒对19秒!感谢您的答案和更新! - hyperknot

5

向量化

除了已经提到的内容,你还应该了解向量化涉及并行化,即对数据进行并发操作,而不是顺序执行(考虑SIMD指令),有些情况下甚至可以利用线程和多处理器...

MEX文件

尽管“解释性 vs. 编译性”这一点已经有人争论过了,但没有人提到你可以通过编写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( files{1} );

%# 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( files{i} );
    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

5

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

嵌套循环本身并不是问题所在,而是您在内部执行的操作。

每次函数调用(特别是非内置函数)都会产生一些开销;如果函数执行例如错误检查等操作,无论输入大小如何,它所需的时间相同,那么开销就更大。因此,如果一个函数只有1毫秒的开销,如果您调用它1000次,您将浪费1秒钟。如果您可以调用它一次来执行矢量化计算,则只需支付一次开销。

此外,JIT编译器(pdf)可以帮助矢量化简单的for循环,例如仅执行基本算术运算的循环。因此,在您的帖子中进行简单计算的循环速度大大加快,而调用median的循环则没有。


2
在这种情况下。
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调用的输出是一组行向量的单元数组。) - hatmatrix

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接