Vectorizing a loop that calculates cosd on every element?

Hi. I have radiation data for every minute of a day that is spectrally resolved so each minute contains several wavelengths. I want to project this radiation on a tilted plane with a nested for loop that looks like this:

for minute=1:1440 %minute during day
     for beta = 1:length(Tilt) %tilt of surface
         for lambda = 1:length(DirectNormal) %wavelength of incoming light

TiltedDirect(minute,beta,lambda)=DirectNormal(lambda,minute)*cosd(AOFINCEQ15(minute,beta));

           end
       end
  end

And what I want to do is go through all minutes of the day and all the tilts (minute and beta) on an element by element basis. BUT I want to perform the same cosine calculation on all wavelengths within that minute and tilt, so I don’t really need to go through all the lambda values per element as I am now. And since I have over 2000 lambdas the loop is incredibly slow. I know very little about vectorization but I assume this would be an ideal situation to vectorize by performing the cosine on all lambda values at the same time. But I can’t figure out how to write that in the code.

Could someone advise me?

Thanks!

Answer 1

[David Young]

You don't need vectorization to avoid computing the cosine on every iteration through the lambda loop - just put the calculation outside:

TiltedDirect = zeros(1440, length(Tilt), size(DirectNormal,1));  
for minute = 1:1440 %minute during day
     for beta = 1:length(Tilt) %tilt of surface
         c = cosd(AOFINCEQ15(minute,beta));
         for lambda = 1:size(DirectNormal,1) %wavelength of incoming light
           TiltedDirect(minute,beta,lambda) = DirectNormal(lambda,minute) * c;
         end
     end
end

Note that I've also changed length to size - you can't guarantee that length will give the size of the right dimensions. I've also preallocated TiltedDirect , which will help efficiency.

By the way, are you sure you want AOFINCEQ15(minute,beta) and not AOFINCEQ15(minute,Tilt(beta))?

If you'd still like to vectorize the inner loop, you can do this

TiltedDirect = zeros(1440, length(Tilt), size(DirectNormal,1));  
for minute = 1:1440 %minute during day
     for beta = 1:length(Tilt) %tilt of surface
         c = cosd(AOFINCEQ15(minute,beta));
         TiltedDirect(minute,beta,:) = DirectNormal(:, minute) * c;
       end
  end

[PetterS]

Brilliant! Thank you so much! Looping though this with the vectorised method you suggested now takes about 30 seconds, whereas before it took hours!!

However, the CPU load during this is still only around 30%. Ideally I would like to perform around 130 000 of these loops, and that would still take over one month to complete. Is there any way I can speed this up? I have heard about parfor instead of ordinary for loop, but it doesn’t seem to work in nested situations.

And if I write AOFINCEQ15(minute,Tilt(beta)) I get an error since the first element in Tilt is 0, which means it’s trying to access AOFINCEQ15(1,0)

[Jan]

The main advantage is the pre-allocation here.

Is AOFINCEQ15 a function or an array?

[David Young]

Just to be clear - I have no real reason to believe that there's anything wrong with AOFINCEQ15(minute,beta). My query was just a suggestion of something to check because the contents of the array Tilt are not actually used in the code you showed.

[PetterS]

AOFINCEQ15 is a 1440x181 matrix, and I agree with you – it is strange that I’m not actually using the values within Tilt. The Tilt variable is simply this: “Tilt=(0:1:180)';”, equally spaced values between 0 and 180. But AOFINCEQ15 has been calculated using Tilt so think AOFINCEQ15(:,1) corresponds to a tilt of 0°, so if I’m not too confused here I guess that would mean that beta=1 is the same as Tilt 0°.

[David Young]

Yes, that makes sense.