Vectorize the code of a summation written with nested for loops

Hi everyone,

I want to implement the following inequality in an efficient form (it's part of an optimization problem I'm solving with the genetic algorithm from the Global Optimization Toolbox):

I still realised it straight with some for loops, like you can see in the code:

for k = 1:K
    summation = 0;
    for w = 1:W
        for m = 1:M
            summation = summation + a(m) * y(m,k,w);
        end
    end
    ineq(k) = summation - c(k);
end

And because Matlab is slow with loops, I want to vectorize my code. This is how far i got:

for w = 1:W
    summation(w,k) = a(m) * y(m,k,w);
end
summation = sum(summation);
ineq = summation - c;

Finally I really don't how to replace the last for loop. It would be great if someone has a solution for it or another idea how to implement the inequality in vectorized form.

[Jonas Kreich]

Unfortunately I'm stucking again with this summation (it's part of my fitness function):

My code with for loops looks like this:

summand = 0; 
for w = 1:W
    for k = 1:K
        for l = 1:K
            for m = 1:M
              for n = 1:M
                    summand = summand + dis(k,l) * trcst(m,n) * y(m,k,w) * y(n,l,w);
                end
            end
        end
    end
end

With W=3, K=10, M=21. So far it works fine, but now i also want to reformulate the code without the for loops and thats the point I'm stucking.

Answer 1

[Guillaume]

I'm assuming:

• a is a column vector. If not, reshape it into a column vector with

   a = a(:);  %or a = reshape(a, [], 1);

• c is a row vector. If not,

c = reshape(c, 1, []);

• M, K, W are the dimensions of y (i.e. [M, K, W] = size(y)), M is also the length of a, and K is also the length of c. If not, crop the relevant vector / matrix so it's the correct size

y = y(1:M, 1:K, 1:W); 
a = a(1:M);
c = c(1:K);

With these assumption fulfilled, it's a one liner using scalar expansion and sum with the correct dimension:

If using R2016b,

ineq = sum(sum(a .* y, 1), 3) - c

If using earlier versions,

ineq = sum(sum(bsxfun(@times, a, y), 1), 3) - c

[Jonas Kreich]

Wow, actually it's really simple :D

Thanks a lot Guillaume!

[Guillaume]

For your new problem, you basically have 5 loops, so you have to reshape your data along 5 dimensions. The only tricky bit is working out, which dimension each variable should be mapped to.

If you map w, k, l, m, n along the 1st, 2nd, ... 5th, respectively, then:

• dis should be a 1xKxL matrix. If it starts as a KxL matrix, then:

dis = permute(dis, [3 1 2]);

• tcrts should be a 1x1x1xMxN matrix. If it starts as a MxN matrix, then:

trcts = permute(trcts, [3 4 5 1 2]);

• the first y should be a WxKx1xM matrix, if it starts as a MxKxW matrix, then:

y1 = permute(y, [3 2 4 1]);

• the second y should be a Wx1xLx1xN matrix, if it starts as a LxNxW matrix, then:

y2 = permute(y, [3 4 1 5 2]);

At which point you can let scalar expansion do its business and:

summand = dis .* tcrts .* y1 .* y2;
summand = sum(summand(:));

Answer 2

[Jonas Kreich]

I nearly answered the new problem from the comment I added to my question by myself.

dis = reshape(dis', numel(dis), 1);
trcst = reshape(trcst', 1, numel(trcst));
for w = 1:3
    tmp(:,:,w) = kron(y(:,:,w), y(:,:,w));
    summand(w) = trcst * tmp(:,:,w) * dis;
end
summand = sum(summand);

Thats how far I got. But now I'm stucking again about the last small for loop.