Is this a valid vectorization of my loop?

I have a loop that I want to calculate that looks like this:

TiltedGroundreflected = zeros(1440, length(Tilt), size(UpScaledBotGlobal,2));  
for minute = 1:1440
    for beta = 1:length(Tilt) 
        for lambda=1:length(Template) 
TiltedGroundreflected(minute,beta,lambda)=GrassReflectance(lambda)*UpScaledBotGlobal(minute,lambda)*(1-cosd(Tilt(beta)))/2;
        end
    end
end

It takes extremely long time to run because of the third for loop (the lambda one), so I wanted to see if I could speed it up. Each lambda represents a wavelength and the loop gets a corresponding reflectance value of grass for each certain wavelength and multiplies with an incident wavelength coming from the sun stored in UpScaledBotGlobal(minute,lambda). So since the wavelength is changing each iteration Im not sure if it can be vectorised or it has to be slow? Anyway, I tried to get rid of the lambda and speed it up by rewriting it like this:

TiltedGroundreflected = zeros(1440, length(Tilt), size(UpScaledBotGlobal,2));  
for minute = 1:1440
    for beta = 1:length(Tilt) 
        c=GrassReflectance(:)'.*UpScaledBotGlobal(minute,:);
TiltedGroundreflected(minute,beta,:)=c*(1-cosd(Tilt(beta)))/2;      
    end
end

And the calculation did go at least 100 times faster, which is great. But since I’m new to vectorization and don’t really 100% know what I’m doing (especially with the “.*”, that symbol is basically just something I try when Matlab won’t accept normal multiplication) my question is – is this correct? Am I calculating the same thing here or should I have written it in a different way?

Answer 1

[Roger Stafford]

I think there is still inefficiency in your computation. For every different beta, you repeat the same computation

   c=GrassReflectance(:)'.*UpScaledBotGlobal(minute,:);

That should have been performed before entering the for beta loop. Also for every different minute (1440 of them) the same computation

   (1-cosd(Tilt(beta)))/2

is performed.

See if the following would be faster:

   G2 = reshape(GrassReflectance,1,1,[]);
   U2 = reshape(UpScaledBotGlobal,1440,1,[]);
   T2 = reshape((1-cosd(Tilt))/2,1,[],1);
   TiltedGroundreflected = bsxfun(@times,bsxfun(@times,G2,U2),T2);

(You could put all that in one (very long) command.)

[PetterS]

Good god that was fast!

But Ive never seen anything like that, I dont understand what you are doing in the code, you are not even looping though anything :S

Could you explain how this code works or is the method called anything particular that I can read about? I have a lot of similar loops that are taking a really long time to run, can this be implemented on all nested for loops? That would potentially revolutionize my life. Like this one:

TiltedDiffuse = zeros(1440, length(Tilt), size(UpScaledBotDiffuse,2));  
for minute = 1:1440
    for beta = 1:length(Tilt)  
TiltedDiffuse(minute,beta,:)=UpScaledBotDiffuse(minute,:)*(1+cosd(Tilt(beta)))/2;
    end
end

[Stephen23]

There is plenty of documentation and best-practice for writing fast MATLAB code. You might as well start with what MATLAB themselves advise:

http://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html

http://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

Basically you need to think about MATLAB as a high-level language that operates on matrices/arrays, and write your code to suit: vectorization is an important part of this concept, and this is what makes Roger Stafford's code many times faster than your own.

Writing MATLAB code as if it were a poor cousin of C and trying to use the same code style is not a good path to writing fast and neat MATLAB code: instead it is better to learn to avoid loops, learn to vectorize calculations, learn to use logical indexing and learn to read the documentation.

You can also find a few good tutorials and guides for writing fast code in MATLAB, such as this:

http://de.mathworks.com/matlabcentral/fileexchange/22943-guidelines-for-writing-clean-and-fast-code-in-matlab

[Roger Stafford]

It's easier to explain using your second example. This code:

   TiltedDiffuse = zeros(1440, length(Tilt), size(UpScaledBotDiffuse,2));  
   for minute = 1:1440
     for beta = 1:length(Tilt)  
       TiltedDiffuse(minute,beta,:) = ...
       UpScaledBotDiffuse(minute,:)*(1+cosd(Tilt(beta)))/2;
     end
   end

can be done by this code:

   TiltedDiffuse = bsxfun(@times,reshape(U,1440,1,[]),...
     reshape((1+cosd(T(beta)))/2,1,[],1));

The inefficiencies in the first code are of the same type as I mentioned earlier. For the 'minute' indices, all 1440 of them, you are recalculating each value of (1+cosd(Tilt(beta)))/2 1440 times. That's an enormous inefficiency, by a factor of 1440 to 1. If you were doing that cosine calculation by hand, it would soon become apparent that you were repeating yourself needlessly. Each of those calculations ought to be saved for multiple use instead of repeated. That is often what is accomplished by the use of the bsxfun function. The 'reshape' function allows you to deal with placing data along appropriate dimensions. I would suggest that you read about the 'bsxfun' and 'reshape' functions very carefully at:

   http://www.mathworks.com/help/matlab/ref/bsxfun.html

and

   http://www.mathworks.com/help/matlab/ref/reshape.html