How to replace all for-loops ?

I had posted a function here for vectorization. Now this is the same function and even I tried myself to replace all for-loops like previous but I didn't succeed because in this the outer loop has enclosed all the for-loops and instead of pop being a vector, now it is a matrix due to which I failed to do so. The code is:

clear all; clc
u=[1  2 10 20]; 
dim=length(u);
pop=rand(10,dim);
C = size(pop,2);
P=C/2;
M=2*C;

e = zeros(size(pop,1),1);
for ii = 1:size(pop,1)
    % calculate xo
    b = pop(ii,:);
    xo=zeros(1,M);
    for k=1:M
        for i=1:P
            xo(1,k)=xo(1,k)+u(i)*exp(-1i*(k-1)*pi*cos(u(P+i)));
        end
    end
    
    % calculate xe
    xe=zeros(1,M);
    for k=1:M
        for i=1:P
            xe(1,k)=xe(1,k)+b(i)*exp(-1i*(k-1)*pi*cos(b(P+i)));
        end
    end
    
    abc=0.0;
    for m1=1:M
        abc=abc+(abs(xo(1,m1)-xe(1,m1))).^2; 
    end
    abc=abc/M; 
    
    e(ii)=abc
    
end

[Walter Roberson]

Do the computations involve complex values? I am wondering why you square the absolute value of the difference, instead of squaring the difference? The result would be the same unless there are imaginary parts.

[Sadiq Akbar]

Thanks a lot dear Walter Roberson for your prompt and kind response. Yes the values of xo and xe both are complex values.

[Eric Delgado]

I don't know if it's possible to remove all loops, but you can easyly remove three of your loops.

u   = [1  2 10 20]; 
pop = rand(10, numel(u));

P = width(pop)/2;
M = width(pop)*2;
H = height(pop);

e = zeros(H,1);
for ii = 1:H
    b = pop(ii,:);

    xo=zeros(1,M);
    xe=zeros(1,M);
    abc = 0;

    for jj=1:M
        for kk=1:P
            xo(1,jj) = xo(1,jj) + u(kk) * exp(-1i*(jj-1) * pi * cos(u(P+kk)));
            xe(1,jj) = xe(1,jj) + b(kk) * exp(-1i*(jj-1) * pi * cos(b(P+kk)));
        end
        abc = abc+(abs(xo(1,jj)-xe(1,jj))).^2; 
    end
    
    e(ii) = abc/M;    
end
e

e = 10×1

5.8441 5.0435 5.2775 5.3064 5.0977 5.5546 5.2565 5.0821 5.6409 5.3807

[Sadiq Akbar]

Thanks a lot dear Eric Delgado for your kind response. Yes, it works but if you look at my question, I want to replace all for-loops using the concept of vectorization. But your code has still for-loops.

[Rik]

Why exactly do you want to replace all for loops?

[Sadiq Akbar]

Thanks for your kind response dear Rik. The reason for replacing all for-loops is to reduce its machine execution time.

[Rik]

For loops in Matlab are actually quite fast, as long as you properly pre-allocate any large array. The main exception is if there is a direct function to do something. Note that arrayfun and cellfun will just hide the loop and will therefore be even slower than a normal loop (since hiding the loop creates additional overhead).

You should note that there is an interdependence between iterations, so the chances of finding a direct function are not looking good.

[Dyuman Joshi]

To add on @Rik's comment, I vectorised a part of the code, and here is the comparison.

You can see that double for loop is faster than the vectorized code.

u   = [1  2 10 20]; 
pop = rand(10, numel(u));

P = width(pop)/2;
M = width(pop)*2;
H = height(pop);

e1=vector(P,M,H,pop,u);
e2=twofor(P,M,H,pop,u);

%comparing results from both methods
result_comparison=isequal(e1,e2)

result_comparison = logical
   1

tic
for i=1:1e3
    y=vector(P,M,H,pop,u);
end
t1=toc;
fprintf('time taken by vector method = %0.4f seconds', t1)

time taken by vector method = 0.1316 seconds

tic
for i=1:1e3
    z=twofor(P,M,H,pop,u);
end
t2=toc;
fprintf('time taken by double for loop method = %0.4f seconds', t2)

time taken by double for loop method = 0.0789 seconds

function e = vector(P,M,H,pop,u)
e=zeros(H,1);
xo=zeros(1,M);
xe=zeros(1,M);
for ii = 1:H
    b=pop(ii,:);

    xo=(exp(-1i*(0:M-1)'*pi*cos(u(P+1:2*P)))*u(1:P)')';
    xe=(exp(-1i*(0:M-1)'*pi*cos(b(P+1:2*P)))*b(1:P)')';
    e(ii)=sum(abs(xo-xe).^2)/M;    
end
end

function e = twofor(P,M,H,pop,u)
e = zeros(H,1);
for ii = 1:H
    b = pop(ii,:);
    xo = zeros(1,M);
    xe = zeros(1,M);
    abc = 0;
    for jj=1:M
        for kk=1:P
            xo(1,jj) = xo(1,jj) + u(kk) * exp(-1i*(jj-1) * pi * cos(u(P+kk)));
            xe(1,jj) = xe(1,jj) + b(kk) * exp(-1i*(jj-1) * pi * cos(b(P+kk)));
        end
        abc = abc+(abs(xo(1,jj)-xe(1,jj))).^2; 
    end
    
    e(ii) = abc/M;    
end
end

[Sadiq Akbar]

Thank you very much dear Dyuman Joshi for your kind response. Yes you are right. But again thee is a problem. When I run your code i.e.

tic
for i=1:1e3
    y=vector(P,M,H,pop,u);
end
t1=toc


tic
for i=1:1e3
    z=twofor(P,M,H,pop,u);
end
t2=toc

Then the execution time of twofor() is less than vector(). But when I run it like the following:

tic
e1=vector(P,M,H,pop,u);
time1=toc

tic
e2=twofor(P,M,H,pop,u);
time2=toc

Then execution time of vector() is less than twofor(). Now what to do?

[Rik]

Using tic toc will only give a general idea of execution times. The inherent variation in execution times may be large enough in this case that the order may affect which version wins if you only call them once.

What you should actually use is the timeit function.

[Sadiq Akbar]

Thanks a lot dear Rik for your kind response. But I didn't heard about it before. How is it used in my case please?

[Dyuman Joshi]

As @Rik mentioned, tic toc will only give an approximate idea of execution time. Results from running tic-toc a single time will always vary. That is why I ran in individual for loops to get a better idea of their execution time.

If you want to see the stark difference between run times, increase the number of iterations. I have increased the count to 1e5. You can see the (obvious) difference in execution times now.

u   = [1  2 10 20];
pop = rand(10, numel(u));

P = width(pop)/2;
M = width(pop)*2;
H = height(pop);

e1=vector(P,M,H,pop,u);
e2=twofor(P,M,H,pop,u);

%comparing results from both methods
result_comparison=isequal(e1,e2)

result_comparison = logical
   1

tic
for i=1:1e5
    y=vector(P,M,H,pop,u);
end
t1=toc;
fprintf('time taken by vector method = %0.4f seconds', t1)

time taken by vector method = 9.4559 seconds

tic
for i=1:1e5
    z=twofor(P,M,H,pop,u);
end
t2=toc;
fprintf('time taken by double for loop method = %0.4f seconds', t2)

time taken by double for loop method = 5.7637 seconds

function e = vector(P,M,H,pop,u)
e=zeros(H,1);
xo=zeros(1,M);
xe=zeros(1,M);
for ii = 1:H
    b=pop(ii,:);

    xo=(exp(-1i*(0:M-1)'*pi*cos(u(P+1:2*P)))*u(1:P)')';
    xe=(exp(-1i*(0:M-1)'*pi*cos(b(P+1:2*P)))*b(1:P)')';
    e(ii)=sum(abs(xo-xe).^2)/M;    
end
end

function e = twofor(P,M,H,pop,u)
e = zeros(H,1);
for ii = 1:H
    b = pop(ii,:);
    xo = zeros(1,M);
    xe = zeros(1,M);
    abc = 0;
    for jj=1:M
        for kk=1:P
            xo(1,jj) = xo(1,jj) + u(kk) * exp(-1i*(jj-1) * pi * cos(u(P+kk)));
            xe(1,jj) = xe(1,jj) + b(kk) * exp(-1i*(jj-1) * pi * cos(b(P+kk)));
        end
        abc = abc+(abs(xo(1,jj)-xe(1,jj))).^2; 
    end
    
    e(ii) = abc/M;    
end
end

[Sadiq Akbar]

Thanks a lot for your kind respone dear Dyuman Joshi . Yes, you are right but how will we do it using timeit function as Rik pointed out above?

[Sadiq Akbar]

Thanks a lot for your kind respone dear Dyuman Joshi . Yes it's ok now.

[Sadiq Akbar]

I don't see any "Accept this answer" button. Where is that?

[Dyuman Joshi]

Because there is no answer, only comments.

Do you want to accept my comment as answer?

[Sadiq Akbar]

Thanks a lot dear Dyuman Joshi for your kind response. Yes, i want to accept your answer and give vote to others also but how because vote icon is not visible with others?

[Jan]

I've moved Dyuman Joshi's comment to the answers section, such that it can be accepted now.

Answer 1

[Dyuman Joshi]

That would be something like this -

e1=vector;
e2=twofor;

%comparing results from both methods
result_comparison=isequal(e1,e2)

result_comparison = logical
   1

fprintf('time taken by vector method = %e seconds', timeit(@vector))

time taken by vector method = 3.112200e-04 seconds

fprintf('time taken by double for loop method = %e seconds', timeit(@twofor))

time taken by double for loop method = 2.392200e-04 seconds

function e = vector
rng(1)
u   = [1  2 10 20];
pop = rand(10, numel(u));

P = width(pop)/2;
M = width(pop)*2;
H = height(pop);

e=zeros(H,1);
xo=zeros(1,M);
xe=zeros(1,M);
for ii = 1:H
    b=pop(ii,:);

    xo=(exp(-1i*(0:M-1)'*pi*cos(u(P+1:2*P)))*u(1:P)')';
    xe=(exp(-1i*(0:M-1)'*pi*cos(b(P+1:2*P)))*b(1:P)')';
    e(ii)=sum(abs(xo-xe).^2)/M;    
end
end

function e = twofor
rng(1)
u   = [1  2 10 20];
pop = rand(10, numel(u));

P = width(pop)/2;
M = width(pop)*2;
H = height(pop);

e = zeros(H,1);
for ii = 1:H
    b = pop(ii,:);
    xo = zeros(1,M);
    xe = zeros(1,M);
    abc = 0;
    for jj=1:M
        for kk=1:P
            xo(1,jj) = xo(1,jj) + u(kk) * exp(-1i*(jj-1) * pi * cos(u(P+kk)));
            xe(1,jj) = xe(1,jj) + b(kk) * exp(-1i*(jj-1) * pi * cos(b(P+kk)));
        end
        abc = abc+(abs(xo(1,jj)-xe(1,jj))).^2; 
    end
    
    e(ii) = abc/M;    
end
end