Vectorization question (trying to avoid for loops)

Dear all, I know that this must be quite a common question and I am sorry if the answer has already been posted somewhere, but after a good look all over this website and several other websites I still am not sure whether there is a way for me to avoid 'for' loops in my case. So here is part of my code which I am trying to vectorize:

G=zeros(t_max+1,9); %each column represents the proportion of people using something numbered 1 to 9. So: sum(G,2)=1
P=zeros(t_max+1,4); %each column represents the proportion of people using either item 1, 2, 3 or 4. So: sum(P,2)=1
SS=[1 1 2 1 1 2 3 3 4]; %each number corresponds to each column of P (item 1 to 4)
I=[a b c d e f  g h i]; %with a to i, 9 preallocated numbers (basically initial proportion values)
G(1,:)=[I];

for t=1:t_max
%G(t+1,:) is updated but not shown here 
 for ii=1:4 %columns of P
    for jj=1:9 %columns of G
       if SS(jj)==ii
          P(t+1,ii)=P(t+1,ii)+G(t+1,jj);
       end
       clear jj;
    end
    clear ii;
 end
end

Basically in plain english (well I will try my best explaining it), at each time step 't' I am updating the proportion in 'P' which is a sum of specific proportions in 'G' (I have not shown this but G(t+1,:) is calculated before the loop). Each proportion of P is updated depending on which column of P the numbers in 'SS' are specifying to. So if SS(:,6) specifies 2, we know that the proportion in the sixth column of G will be added to second column of P (cumulative sums). (I think that might be even less understandable so if you have any question just ask me).

Sorry for the really long post and hopefully someone will be able to tell me if I really need all these 'for' loops or if I can simplify it using vectorization.

Thank you!

Regards

Jonathan

[Adam]

I'm afraid I don't have time to consider it in enough detail though if someone else hasn't answered when I next have a spare moment I will take a look.

What caused me to rely though is to ask why do you have

clear jj

in the middle of a loop?! I have never tried it so I don't know what the effect is, but I cannot imagine it being a good effect. It would either mess up the loop entirely or do nothing if the loop variable is kept track of behind the scenes anyway (I'm not sure of the exact mechanics of how that works).

[dpb]

As is effectively it does nothing; the variable does disappear for the time between there and the next iteration but it then will be recreated and will have the proper value as the for expression is evaluated prior to the loop and not reevaluated.

That said, it is "poor practice" and should be avoided. Both of those clear statements should be deleted by OP simply as being unneeded and unwise.

[Jonathan E.]

I was trying to remove any unnecessary variables in the computation but I did not know how for loops work so was not sure whether it was really useful or not. Now I know so thanks. I will just add the clear argument after the loop. Thanks for your help!

Answer 1

[Kirby Fears]

I deleted the clear lines and switched the ( t_max+1) tendency to simply be t_max. I added t_max and I values up top so it runs on my machine. Here is the code with those changes:

t_max = 5; % needed a t_max value
G = zeros(t_max,9);
P = zeros(t_max,4);
SS = [1 1 2 1 1 2 3 3 4];
I = 1:9; % needed real I values
G(1,:) = I;

for t = 1:t_max,
    for ii = 1:4,
        for jj = 1:9,
            if SS(jj)==ii,
                P(t,ii) = P(t,ii) + G(t,jj);
            end
         end
    end
end

Then I simplified your nested loops to a single loop.

for s = 1:numel(SS),
    P(:,SS(s)) = P(:,SS(s)) + G(:,s);
end

This gives the same result as the nested loops, though G and P are mostly full of zeros still. Let us know if this is what you needed.

[Jonathan E.]

Dear Kirby Fears, Sorry for my late reply and thank you so much for your answer. It is exactly what I needed! That has improved the speed of my model already. It seems so simple now that you have outlined the way to do it... Thanks a lot.

[Star Strider]

The sincerest expression of appreciation here on MATLAB Answers is to Accept the Answer that most closely solves your problem.

Answer 2

[Jonathan E.]

Well actually if I may take a bit more of your time, I am also wondering how to simplify this one (I am really trying to figure it out but am a bit stuck to be honest). I have the same G matrix of sum 1 (to simplify it here I have considered it as an array). For example:

 G=[0.2 0.1 0.05 0.05 0.1 0.3 0.1 0.05 0.05]

I have a 9x9 W matrix (not shown here) which represents the probability for each 1 to 9 individuals to meet one each other.

I was wondering whether these for loops could be simplified:

 E=zeros(1,9);
 for ii=1:9
     for jj=1:9
          E(ii)=E(ii)+G(1,ii)*G(1,jj)*W(SS(ii),SS(jj));
     end
 end

Sorry if that seems obvious again...

Jonathan

[Kirby Fears]

Is SS the same as before? I would guess W is indexed with W(ii,jj), not W(SS(ii),SS(jj)). Please confirm and I can try to vectorize this.

[Kirby Fears]

Here is a simplification assuming W(ii,jj) was the proper indexing.

G = [0.2 0.1 0.05 0.05 0.1 0.3 0.1 0.05 0.05];
W = rand(numel(G));
E = zeros(size(G));
for ii = 1:numel(E),
    E(ii) = E(ii) + sum(G(ii)*G.*W(ii,:));
end

Here is the same thing with matrix multiplication instead.

G = [0.2 0.1 0.05 0.05 0.1 0.3 0.1 0.05 0.05];
W = rand(numel(G));
E = G.*(G*W');