可以进行矢量化吗？会有什么好处吗？

Question

这个功能在我的跑步中吃了很多时间.但是看到的是内置功能的大部分时间polyarea.这段代码可以进行矢量化以提升性能吗？

Profiler报告 -

  time   calls
                  1 function [S S_area] = Polygons_intersection_Compute_area(S)
                  2 % Guillaume JACQUENOT
                  3 % guillaume at jacquenot at gmail dot com
                  4 % 2007_10_08
                  5 % 2009_06_16
                  6 % Compute area of each polygon of in S.
                  7 % Results are stored as a field in S
                  8 
  0.50   51945    9 S_area = struct('A', {}); 
  0.20   51945   10 for i=1:numel(S) 
  0.28  103890   11     S(i).area = 0; 
  1.34  103890   12     S_area(i).A = zeros(1,numel(S(i).P)); 
  0.69  103890   13     for j=1:numel(S(i).P) 
  9.24  103890   14         S_area(i).A(j) = polyarea(S(i).P(j).x,S(i).P(j).y); 
  0.28  103890   15         S(i).area      = S(i).area + (1-2*S(i).P(j).hole) * S_area(i).A(j);         
  0.01  103890   16     end 
  0.08  103890   17 end 

Answer 1

我看到4个问题.我将按照潜在性能增益的增加顺序讨论它们.

首先:您使用i和j作为循环变量名称.这些也是Matlab中虚构单元的名称,这意味着Matlab将不得不花一些时间查找你的意思.事情是,如果循环不是JIT(它不是你的,我会得到它),它必须在每次迭代时都这样做.

第二:索引多维结构需要比你想象的更多的时间.在这方面,Multi-D结构有些臭名昭着,你最好避免对它们进行过多的索引操作.通常制作元素的简单副本,对该副本执行所有操作,然后将副本写回结构可以提高性能.

第三:你没有S_area以最有效的方式预先分配.您甚至不预先分配结构,而是在预分配时在第一个循环中增长S_area(i).A.这一切都可以改进(见下文).

第四:polyarea是不是内置的功能,所以这种双循环将不JIT'ed.如果你调用循环中的任何函数,你或Mathworks用M语言(而不是C)编写,JIT编译器将无法编译你的循环.这是迄今为止 JIT框架中最烦人(和可改进的)限制,而JIT的循环通常比非JIT的循环运行速度快100倍或更快.

唯一的解决方案通常是在循环体中"内联"非内置函数.在Matlab中意味着:将函数体的全部内容复制粘贴到循环中,并对该体中调用的所有非内置函数递归执行此操作.

所有上述内容都会导致您的代码版本:

% pre-allocate S_area
S_area(numel(S)).A = [];
As = cellfun(@(x) zeros(numel(x),1), {S.P}, 'UniformOutput', false);
[S_area.A] = deal(As{:});

% number of polygons for all S(ii)
numPolys = cellfun(@numel, {S.P});

% enter loop
for ii = 1:numel(S)
    % extract S(ii) only once
    Sii = S(ii);

    Sii.area = 0;
    Aii = S_area(ii).A;        
    for jj = 1:numPolys(ii)

        p = Sii.P(jj);  % extract polygon only once
        x = p.x; % and its x and y components
        y = p.y;            
        sz = size(p);

        % NOTE: core of polyarea. Note that all checks and flexibility, and 
        % therefore user-friendliness, is GONE. Very little has to go wrong 
        % here before a hard-to-understand error is issued. 
        Area = reshape(abs(sum( (x([2:sz(1) 1],:) - x(:,:)).* ...
            (y([2:sz(1) 1],:) + y(:,:)))/2),[1 sz(2:end)]);

        Aii(jj) = Area;
        Sii.area = Sii.area + Area*(1-2*p.hole);
    end

    % place copies back into the strucure
    S_area(ii).A = Aii;
    S(ii).area = Sii.area;

end

我无法尽可能正确地测试这个,所以如果你发现一些错误,请告诉我,我会尽力纠正它们.