Vectorized code slower than loops? MATLAB

Question

In the problem Im working on there is such a part of code, as shown below. The definition part is just to show you the sizes of arrays. Below I pasted vectorized version - and it is >2x slower. Why it happens so? I know that i happens if vectorization requiers large temporary variables, but (it seems) it is not true here.

And generally, what (other than parfor, with I already use) can I do to speed up this code?

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  

Best regards,
Alex

From the profiler:

Edit:

In reply to @Jason answer, this alternative takes the same time:

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:
In reply to @EBH :
The point is to do the following:

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

I speeded it up as much, as I was able to. What you have written is taking only the last bessels and posMcontainer values, so it does not produce the same result. In the real code, those two containers are filled not with 1, but with some precalculated values.

Answer 1

After your edit, I can see that umn is just a temporary variable for another calculation. It still can be mostly vectorizable:

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

This reduce running time from ~95 seconds to less 3 seconds (both without parfor), so it improves in almost 97%.

Answer 2

I would suspect it is memory allocation. You are re-allocating the m array in a 3 deep loop. try rearranging the code:

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

I was trying this in Igor pro, which a similar language, but with different optimizations. So the direct translations don't time the same way as Matlab (vectorized was slightly faster in Igor). But reordering the loops did speed up the vectorized form.

In your second part of the code, that is setting umn2, inside the loops, you have:

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

Those 3 lines don't require any input from the i and j loops, they only use the n loop. So reordering the loops such that i and j are inside the n loop will mean that those 3 lines are done xElements^2 (100^2) times less often. I suspect it is that m = 1:2:n line that takes time, since that is allocating an array.