How can you avoid for loops by using vectorization?

I am fairly new to matlab. I am used to coding in C that is why I use a lot for loops but recently I encounter a problem. I am working on assignment for school where I have to run monte carlo simulations for different diversity schemes. I got the simulations to work correctly. However, I am trying to implement this because it takes a very long time to run simulation since I am using a lot for loops. So, I was wondering if someone can help me out with this to avoid for loops and start using vectorization instead, since I know matlab will run faster. below is my code:

    stepSize = 1;
      M = (1:stepSize:4);
  rho = (1:stepSize:18);
  S = 1;

disp('Computing...');drawnow;
tic;

for i = 1:length(M)

     for j = 1:length(rho)

              Sum_MRC = 0;

            for k = 0:(M(i)-1)

               C = nchoosek((M(i)-1+k),k);
               gamma = sqrt(rho(j)/(1+rho(j)));
               Sum_MRC = Sum_MRC + (C * ((1+gamma)/2)^k);

            end

          MRC_Theoretical(i,j) = ((1-gamma)/2)^M(i) * Sum_MRC; 
          Upper_bound(i,j)= 0.5*(1/(1+rho(j))^M(i));

          Error_Counter = 0;
          Total_Counter = 0;
        while(Error_Counter<100)
              hr = randn(1,M(i));
              hi = randn(1,M(i));
              vr = randn(1,M(i));
              vi = randn(1,M(i));
              Hm = (hr + (1i*hi));
              Vm = (vr + (1i*vi));

                  Xm = ((sqrt(rho(j))*Hm*S + Vm)/sqrt(2)); %Construction of thefading model channel used.
              MRC_Decision = 0;
              for k = 1:M(i)
                  MRC_Decision = MRC_Decision + (Xm(k)*conj(Hm(k)));
                  %Estimating the decision made by the MRC diversity technique.
              end
              if (real(MRC_Decision) < 0)
                  Error_Counter = Error_Counter + 1; %Incrementing theerror counter if the MRC decision is not equal to thee one that is sent.
              end

          Total_Counter = Total_Counter + 1;
        end
          Error_Probability(i,j) = Error_Counter/Total_Counter; %Calculating the probability of error based on the Monte Carlo simulations.
    end
  end
  toc;

%Plot for the Monte-Carlo simualtion vs the theoretical values calculated.
figure(1)
semilogy(10*log10(rho),Error_Probability(3,:),'k*:',10*log10(rho),MRC_Theoretical(3,:),'r*:',10*log10(rho),Upper_bound(1,:),'b*:')
xlabel('Signal to Noise ratio');
ylabel('Average probability of error');
legend('Eb for Monte-Carlofor M=3 ','Theoretical Eb','Upperbound','Location','SouthWest');
title('Plot for MRC diversity comparing with Monte Carlo simulation','FontWeight','bold');

%Plot for the all the Theoretical values calculated for all the antennas

figure(2)
semilogy(10*log10(rho),MRC_Theoretical(1,:),'k*:',10*log10(rho),MRC_Theoretical(2,:),'r*:',10*log10(rho),MRC_Theoretical(3,:),'b*:',10*log10(rho),MRC_Theoretical(4,:),'g*:')
xlabel('Signal to Noise ratio');
ylabel('Average probability of error');
legend('Eb for M=1','Eb for M=2','Eb for M=3','Eb for M=4','Location','SouthWest');
title('Plot for MRC diversity for different M values','FontWeight','bold');

Answer 1

[Image Analyst]

It's not the for loops that's causing the problem. Your i loop goes 4 times and your j loop goes 18 times so that's 72 times. The k loop is anywhere from 1 to 4 times. So you have like a hundred iterations. You computer can do tens of millions of for loop iterations in a small fraction of a second. It will do a hundred iterations faster than a greased cheetah on steroids. The problem lies elsewhere in your calculations. But it's situations like this where for loops get improperly blamed for being slow when they have nothing to do with it. Look, I did 100000000, yes 100 million, for loop iterations and this was the result:

Elapsed time is 0.229573 seconds.

That's way more than you're doing and it's really fast. That's just for the for loop itself, of course adding any code would slow it down, but it's the extra code that would slow it down, not the for loop iterating. If you drive your Ferrari from here to Miami and it takes 4 days but you stopped at every stop light and convenience store along the way, you shouldn't accuse the Ferrari of being slow.

[Mostafa Nakhaei]

Actually ForLoops in Matlab is not comparble with ForLoops in C as C make the ForLoops supper efficient. Avoid using ForLoops in Matlab and try using matrix and arrays.

For example, if you want to calculate an integral in Matlab you should make use of matrix multiplication and summing the matrix than using a for loop to calculate the integral.

Thanks

Answer 2

[Henrik]

I'm a bit tired, so here are just a few suggestions.

You could calculate C once and for all outside the loops, this will save time since it's independent of j.

gamma can be calculated for all j at once, also outside the loop:

gamma=sqrt(rho./(1+rho));

This

for k = 1:M(i)
                MRC_Decision = MRC_Decision + (Xm(k)*conj(Hm(k)));
                %Estimating the decision made by the MRC diversity technique.
            end

can also be written as

MRC_Decision=dot(Xm,conj(Hm))

Here's the start of some code that I think will work. I have not tested it, though.

Apart from this, maybe try running your code with the profiler on to see where most of the compution time is spent.

gamma=sqrt(rho./(1+rho));
for i = 1:length(M)

for k=0:(M(i)-1)
C(k)=nchoosek((M(i)-1+k),k);
end
for j = 1:length(rho)
Sum_MRC=dot(C,(1+gamma(j)/2).^k);

etc. Remember that gamma is now a vector, so replace gamma elsewhere in the code with gamma(j)

Answer 3

[SK]

The ratio of TotalCounter to ErrorCounter is extremely high.

For just j = 4, TotalCounter goes up to 456,513. I didn't check how big it is for j = 18.

So the timing will be extremely sensitive to the content of the while loop. You need to optimize that.

1. Generate the random numbers together.

2. You could vectorize

for k = 1:M(i)
  MRC_Decision = MRC_Decision + (Xm(k)*conj(Hm(k)));
end

but since M is small, it may not make too much of a difference.

3. Also, for a small gain, write instead of the if statement, write:

Error_Counter = Error_Counter + (real(MRC_Decision) < 0);

Implementing this:

sqrt_two = sqrt(2);
a = sqrt(rho(j))*S; % assuming S is scalar.

while(Error_Counter < 100)

    h = randn(4,M(i));
    Hm = (h(1,:) + (1i*h(2,:)));
    Vm = (h(3,:) + (1i*h(4,:)));

    % Construction of the fading model channel used.
    Xm = (a*Hm + Vm)/sqrt_two;

    %Estimating the decision made by the MRC diversity technique.
    MRC_Decision = sum(Xm.*conj(Hm));

    %Incrementing the error counter if required.
    Error_Counter = Error_Counter + (real(MRC_Decision) < 0);

    Total_Counter = Total_Counter + 1;

end

On my machine, for j = 4, the time reduces from 16 secs to 10 secs.

But there is probably more you can do. The section:

h = randn(4,M(i));
Hm = (h(1,:) + (1i*h(2,:)));
Vm = (h(3,:) + (1i*h(4,:)));

% Construction of the fading model channel used.
Xm = (a*Hm + Vm)/sqrt_two;

MRC_Decision = sum(Xm.*conj(Hm));

can be simplified by removing the complex number calls - avoid it - it just slows things down in Matlab. Also sqrt(2) doesn't do anything of consequence. So finally:

Ec = 0;
Tc = 0;
a = S*sqrt(rho(j));
while (Ec < 100)

    Tc = Tc + 1;

    Hm = randn(M(i), 2);
    Vm = randn(M(i), 2);

    % Construction of the fading model channel used.
    Xm = (a*Hm + Vm);
        
    %Estimating the decision made by the MRC diversity technique.
    MRC_Decision = sum( Xm(:,1).*Hm(:,1) + Xm(:,2).*Hm(:,2) );

    %Incrementing the error counter if required.
    Ec = Ec + (MRC_Decision < 0);

end

Takes around 2 seconds for j = 4 (down from 16).