矢量化代码比循环慢？MATLAB

Question

在我正在研究的问题中，有这样一部分代码，如下所示。定义部分只是向您展示数组的大小。下面我粘贴了矢量化版本--而且它的速度比以前慢了2倍。为什么会这样？我知道，如果向量化需要大量临时变量，就会发生这种情况，但是(看起来)这里不是这样的。

通常，除了parfor之外，我还能做些什么来加快代码的速度呢？

maxN = 100;  
levels = maxN+1;  
xElements = 101;  
umn = complex(zeros(levels, levels));  
umn2 = umn;  
bessels = ones(xElements, xElements, levels);    % 1.09 GB  
posMcontainer = ones(xElements, xElements, maxN);  

tic  
for j = 1 : xElements  
    for i = 1 : xElements  
        for n = 1 : 2 : maxN  
            nn = n + 1;  
            mm = 1;  
            for m = 1 : 2 : n  
                umn(nn, mm) = bessels(i, j, nn) * posMcontainer(i, j, m);  
                mm = mm + 1;  
            end  
        end  
    end  
end  
toc % 0.520594 seconds  


tic  
for j = 1 : xElements  
    for i = 1 : xElements  
        for n = 1 : 2 : maxN  
            nn = n + 1;  
            m = 1:2:n;  
            numOfEl = ceil(n/2);  
            umn2(nn, 1:numOfEl) = bessels(i, j, nn) * posMcontainer(i, j, m);  
        end  
    end  
end  
toc % 1.275926 seconds  

sum(sum(umn-umn2)) % veryfying, if all done right  

诚挚的问候,

亚历克斯

来自侧写员：

​

​

编辑：

在对@Jason answer的答复中，这个替代方案需要同样的时间：

for n = 1:2:maxN  
    nn(n) = n + 1;  
    numOfEl(n) = ceil(n/2);  
end  

for j = 1 : xElements  
    for i = 1 : xElements  
        for n = 1 : 2 : maxN  
            umn2(nn(n), 1:numOfEl(n)) = bessels(i, j, nn(n)) * posMcontainer(i, j, 1:2:n);  
        end  
    end  
end 

Edit2

答复@EBH：

重点是做以下工作：

parfor i = 1 : xElements  
    for j = 1 : xElements  
    umn = complex(zeros(levels, levels)); % cleaning  
    for n = 0:maxN
        mm = 1;
        for m = -n:2:n
            nn = n + 1; % for indexing

            if m < 0
                umn(nn, mm) = bessels(i, j, nn) * negMcontainer(i, j, abs(m));
            end

            if m > 0
                umn(nn, mm) = bessels(i, j, nn) * posMcontainer(i, j, m);
            end

            if m == 0
                umn(nn, mm) = bessels(i, j, nn);
            end

            mm = mm + 1; % for indexing
        end % m
    end % n
    beta1 = sum(sum(Aj1.*umn));
    betaSumSq1(i, j) = abs(beta1).^2;

    beta2 = sum(sum(Aj2.*umn));
    betaSumSq2(i, j) = abs(beta2).^2;
    end % j
end % i

我尽可能地加快了速度。您所编写的只是最后一个bessels和posMcontainer值，因此它不会产生相同的结果。在实际代码中，这两个容器不是用1填充的，而是用一些预先计算的值填充的。

Answer 1

在编辑之后，我可以看到umn只是另一个计算的临时变量。它仍然可以大部分是可矢量化的：

betaSumSq1 = zeros(xElements); % preallocating
betaSumSq2 = zeros(xElements); % preallocating
% an index matrix to fetch the right values from negMcontainer and
% posMcontainer:
indmat = tril(repmat([0 1;1 0],ceil((maxN+1)/2),floor(levels/2)));
indmat(end,:) = [];
% an index matrix to fetch the values in correct order for umn:
b_ind = repmat([1;0],ceil((maxN+1)/2),1);
b_ind(end) = [];
tempind = logical([fliplr(indmat) b_ind indmat+triu(ones(size(indmat)))]);

% permute the arrays to prevent squeeze:
PM = permute(posMcontainer,[3 1 2]);
NM = permute(negMcontainer,[3 1 2]);
B = permute(bessels,[3 1 2]);

for k = 1 : maxN+1 % third dim
    for jj = 1 : xElements % columns
        b = B(:,jj,k); % get one vector of B

        % perform b*NM for every row of NM*indmat, than flip the result:
        neg = fliplr(bsxfun(@times,bsxfun(@times,indmat,NM(:,jj,k).'),b));

        % perform b*PM for every row of PM*indmat:
        pos = bsxfun(@times,bsxfun(@times,indmat,PM(:,jj,k).'),b);

        temp = [neg mod(1:levels,2).'.*b pos].'; % concat neg and pos
        % assign them to the right place in umn:
        umn = reshape(temp(tempind.'),[levels levels]).';

        beta1 = Aj1.*umn;
        betaSumSq1(jj,k) = abs(sum(beta1(:))).^2;
        beta2 = Aj2.*umn;
        betaSumSq2(jj,k) = abs(sum(beta2(:))).^2;
    end
end

这使运行时间从~95秒减少到3秒(都没有parfor)，因此在几乎97%的中得到了改进。

Answer 2

我怀疑这是内存分配。您正在3深循环中重新分配m数组。尝试重新排列代码：

tic
for n = 1 : 2 : maxN
    nn = n + 1;
    m = 1:2:n;
    numOfEl = ceil(n/2);
    for j = 1 : xElements
        for i = 1 : xElements
            umn2(nn, 1:numOfEl) = bessels(i, j, nn) * posMcontainer(i, j, m);
        end
    end
end
toc % 1.275926 seconds

我在Igor中尝试了这一点，这是一种类似的语言，但是有不同的优化。因此，直接翻译的时间不像Matlab那样(在Igor中矢量化的速度略快)。但是重新排列循环确实加快了矢量化形式的速度。

在代码的第二部分，即设置umn2，在循环中，您有：

nn = n + 1;  
m = 1:2:n;  
numOfEl = ceil(n/2);  

这3行不需要来自i和j循环的任何输入，它们只使用n循环。因此，重新排序循环，例如i和j在n循环中，将意味着这3行的执行次数比xElements^2 (100^2)少了一倍。我怀疑需要时间的是m = 1:2:n行，因为这是分配数组。