Thread Subject: Want Help Vectorizing a For Loop, #3

Dear Matlab community,

I'm hoping that someone with considerable expertise in programming for efficiency can help me vectorize the following for loops.

%#####################################################################
% Initialize variables
beta = linspace(-1000,1000,500);
omega = linspace(0 ,2*pi ,500);
alpha = linspace(-1000,1000,500);
Spectrum = zeros(500 ,500 ,500);
%#####################################################################


% Begin For Loops
%#####################################################################
for indexOmega = 1:500
     for indexBeta = 1:500
          for indexAlpha = 1:500
Spectrum(indexOmega,indexBeta,indexAlpha) = (abs((mcgama(1/2,i*omega(indexOmega))*cos((1/2+i*omega(indexOmega))*atan(beta(indexBeta)*alpha(indexAlpha))))/((alpha(indexAlpha))^2+(beta(indexBeta))^2)^(1/4+i*omega(indexOmega)/2)))^2;
          end
     end
end
%#####################################################################

It will currently take about 192 minutes total or 23 seconds per indexOmega loop.

Thanks,

Eric Diaz

Dear Matlab community,

I just wanted to add that I know this is my 3rd vectorization help request and it may appear like I am just being lazy. However, I really just don't get exactly how to do it. It's a pretty difficult concept to get hold of right away. I'm hoping that with time, I will get better at it (maybe by the 7th or 8th time it usually takes me). Anyway, thanks to everyone in the past for their help and thanks in advance to anyone who helps with this problem.

Eric

"Eric Diaz" <eric.diaz@gmail.com> wrote in message <hnc7em$n7f$1@fred.mathworks.com>...
> Dear Matlab community,
>
> I'm hoping that someone with considerable expertise in programming for efficiency can help me vectorize the following for loops.
>
> %#####################################################################
> % Initialize variables
> beta = linspace(-1000,1000,500);
> omega = linspace(0 ,2*pi ,500);
> alpha = linspace(-1000,1000,500);
> Spectrum = zeros(500 ,500 ,500);
> %#####################################################################
>
>
> % Begin For Loops
> %#####################################################################
> for indexOmega = 1:500
> for indexBeta = 1:500
> for indexAlpha = 1:500
> Spectrum(indexOmega,indexBeta,indexAlpha) = (abs((mcgama(1/2,i*omega(indexOmega))*cos((1/2+i*omega(indexOmega))*atan(beta(indexBeta)*alpha(indexAlpha))))/((alpha(indexAlpha))^2+(beta(indexBeta))^2)^(1/4+i*omega(indexOmega)/2)))^2;
> end
> end
> end
> %#####################################################################
>
> It will currently take about 192 minutes total or 23 seconds per indexOmega loop.
>
> Thanks,
>
> Eric Diaz

Eric Diaz wrote:

> I'm hoping that someone with considerable expertise in programming for
> efficiency can help me vectorize the following for loops.
>
> %#####################################################################
> % Initialize variables
> beta = linspace(-1000,1000,500);
> omega = linspace(0 ,2*pi ,500);
> alpha = linspace(-1000,1000,500);
> Spectrum = zeros(500 ,500 ,500);
> %#####################################################################
>
>
> % Begin For Loops
> %#####################################################################
> for indexOmega = 1:500
> for indexBeta = 1:500
> for indexAlpha = 1:500
> Spectrum(indexOmega,indexBeta,indexAlpha) =
> (abs((mcgama(1/2,i*omega(indexOmega))*cos((1/2+i*omega(indexOmega))*atan(beta(indexBeta)*alpha(indexAlpha))))/((alpha(indexAlpha))^2+(beta(indexBeta))^2)^(1/4+i*omega(indexOmega)/2)))^2;
>
> end
> end
> end

Remove from each loop the part that is constant for the loop. For
example, the mcgama (or is it mcgamma?) call will have a result that is
independant of indexAlpha or indexBeta, so that part can be moved
outside of both of those loops:

for indexOmega = 1:500
   thisomega = i*omega(indexOmega);
   thisgamma = mcgama(1/2,thisomega);
   cosarg1 = 1/2 + thisomega;
   for indexBeta = 1:500
    thisbeta = beta(indexBeta)
    theseatans = atan(thisbeta .* alpha);
     for indexAlpha = 1:500
       Spectrum(indexOmega,indexBeta,indexAlpha) =
(abs((thisomega*cos(cosarg1*theseatans(indexAlpha))) /
etc

As you continue to work at it, you will find that you will be able to
vectorize all of the work done in the indexAlpha loop, and so will be
able to reduce it to two levels with an assignment to
Spectrum(indexOmega,indexBeta,:) . After that, you will be able to do
things like move all of the atan computations outside the indexBeta loop
by using techniques such as

ABatans = arrayfun(@(a) atan(beta .* a), alpha).';

then inside the indexBeta loop, BAatans(:,indexBeta)

(you'll probably have to play a bit with the order of operations and the
transpose to get the alpha as the first column -- remember, first column
is indexed fastest.)

Speaking of first column being indexed most quickly, you should use

Spectrum(indexAlpha,indexBeta,indexOmega) if your innermost loop is
indexAlpha.

Hi Walter,

Wow, thanks so much for your instructional explanation. I really appreciate learning how to do things by having someone kind of go through the thought processes behind steps. I can say that I definitely learned something this time that I will be able to use in the future.

I'll implement your suggestions and re-post a little bit later.

Thanks,

Eric Diaz

Walter Roberson <roberson@hushmail.com> wrote in message <hnceav$b9h$1@canopus.cc.umanitoba.ca>...
> Eric Diaz wrote:
>
> > I'm hoping that someone with considerable expertise in programming for
> > efficiency can help me vectorize the following for loops.
> >
> > %#####################################################################
> > % Initialize variables
> > beta = linspace(-1000,1000,500);
> > omega = linspace(0 ,2*pi ,500);
> > alpha = linspace(-1000,1000,500);
> > Spectrum = zeros(500 ,500 ,500);
> > %#####################################################################
> >
> >
> > % Begin For Loops
> > %#####################################################################
> > for indexOmega = 1:500
> > for indexBeta = 1:500
> > for indexAlpha = 1:500
> > Spectrum(indexOmega,indexBeta,indexAlpha) =
> > (abs((mcgama(1/2,i*omega(indexOmega))*cos((1/2+i*omega(indexOmega))*atan(beta(indexBeta)*alpha(indexAlpha))))/((alpha(indexAlpha))^2+(beta(indexBeta))^2)^(1/4+i*omega(indexOmega)/2)))^2;
> >
> > end
> > end
> > end
>
> Remove from each loop the part that is constant for the loop. For
> example, the mcgama (or is it mcgamma?) call will have a result that is
> independant of indexAlpha or indexBeta, so that part can be moved
> outside of both of those loops:
>
> for indexOmega = 1:500
> thisomega = i*omega(indexOmega);
> thisgamma = mcgama(1/2,thisomega);
> cosarg1 = 1/2 + thisomega;
> for indexBeta = 1:500
> thisbeta = beta(indexBeta)
> theseatans = atan(thisbeta .* alpha);
> for indexAlpha = 1:500
> Spectrum(indexOmega,indexBeta,indexAlpha) =
> (abs((thisomega*cos(cosarg1*theseatans(indexAlpha))) /
> etc
>
> As you continue to work at it, you will find that you will be able to
> vectorize all of the work done in the indexAlpha loop, and so will be
> able to reduce it to two levels with an assignment to
> Spectrum(indexOmega,indexBeta,:) . After that, you will be able to do
> things like move all of the atan computations outside the indexBeta loop
> by using techniques such as
>
> ABatans = arrayfun(@(a) atan(beta .* a), alpha).';
>
> then inside the indexBeta loop, BAatans(:,indexBeta)
>
> (you'll probably have to play a bit with the order of operations and the
> transpose to get the alpha as the first column -- remember, first column
> is indexed fastest.)
>
> Speaking of first column being indexed most quickly, you should use
>
> Spectrum(indexAlpha,indexBeta,indexOmega) if your innermost loop is
> indexAlpha.

Dear folks,

Walter's suggestions worked great! The average time of each outer loop went down to 0.208 seconds and 104sec total! That's >100x increase in speed and efficiency! Thanks Walter!

The only thing I couldn't get to work was that arrayfun thing. It kept on telling me that I had to tell it whether it could use scalar output or something.


Below is what I ended up doing.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize Variables
Spectrum = zeros(500,500,500);
beta = linspace(-1000 ,1000,500);
omega = linspace(0 ,2*pi,500);
alpha = linspace(-1000 ,1000,500);
alphaSquared = alpha.^2;
% theseBoverAaTans = arrayfun(@(a) atan(beta/a),alpha).';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for indexOmega = 1:500

thisOmega = 1i*omega(indexOmega);
thisGamma = mcgama(1/2,thisOmega);
cosarg1 = 1/2 + thisOmega;
denomExpArg = 1/4 + thisOmega/2;

for indexBeta = 1:500

thisBeta = beta(indexBeta);
thisBetaSquared = thisBeta^2;
theseBoverAaTans= atan(thisBeta./alpha);

Spectrum(:,indexBeta,indexOmega) = abs( (thisGamma.*cos(cosarg1.*theseBoverAaTans)) ./ ((alphaSquared+thisBetaSquared).^(denomExpArg)) ).^2;

end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Eric Diaz wrote:

> The only thing I couldn't get to work was that arrayfun thing. It kept
> on telling me that I had to tell it whether it could use scalar output
> or something.

Don't worry too much about that. Probably using

bsxfun(@(a,b) atan(b/a), alpha.', beta)

would be faster anyhow.


> Below is what I ended up doing.
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> % Initialize Variables
> Spectrum = zeros(500,500,500);
> beta = linspace(-1000 ,1000,500);
> omega = linspace(0 ,2*pi,500);
> alpha = linspace(-1000 ,1000,500);
> alphaSquared = alpha.^2;
> % theseBoverAaTans = arrayfun(@(a) atan(beta/a),alpha).';
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> for indexOmega = 1:500
>
> thisOmega = 1i*omega(indexOmega);
> thisGamma = mcgama(1/2,thisOmega);
> cosarg1 = 1/2 + thisOmega;
> denomExpArg = 1/4 + thisOmega/2;
>
> for indexBeta = 1:500
>
> thisBeta = beta(indexBeta);
> thisBetaSquared = thisBeta^2;
> theseBoverAaTans= atan(thisBeta./alpha);

I notice that in the previous version, you had atan(thisBeta .* alpha)
but in the present version you have changed it to division instead. As
alpha includes 0 in its range, you now have a division by 0 where you
didn't before. Generally that would be okay, as the division will
generally produce -inf or +inf and atan() of that is quite meaningful.
Unfortunately, thisBeta can also be 0, so you now have a 0/0 which is
NaN, and you are going to get one of those NaN at the same place for
every Omega (because NaN's "pollute" the rest of the calculation.) In
the previous version with the multiplication, this NaN could not happen.
Are you sure you want this situation?

On the other hand, because your range -1000 to 1000 is 2001 points, and
you are requesting that it be sampled at an even number of points (500),
your ranges will not include 0 exactly: the closest alpha and beta will
get to 0 are +/- 2.00400801603212 . Thus, the NaN described above will
not occur in practice... unless you happen to change your sizes to be
odd, or happen to nudge the upper or lower bound. I recommend adding a
comment warning about the potential consequences.


> Spectrum(:,indexBeta,indexOmega) = abs(
> (thisGamma.*cos(cosarg1.*theseBoverAaTans)) ./
> ((alphaSquared+thisBetaSquared).^(denomExpArg)) ).^2;
>
> end
> end
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Hey Walter,

Thanks for your warning. I actually was pretty greedy initially and I also didn't think up my limits properly. I've since cut back significantly on the size and ranges. This function is actually quite bizaar.

In addition, I fixed up my code because I wasn't getting results consistent with what I wanted. And, actually, it turns out that it's even 10x faster with the inside for loop. Not sure why, but it is.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all force
clear all force
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize variables
numEl = 40;
Spectrum = zeros(numEl,numEl,numEl,numEl);
beta = linspace(-1 ,1 ,numEl);
omega = linspace(0*pi ,2*pi ,numEl);
alpha = linspace(-1 ,1 ,numEl);
% theseBoverAaTans = arrayfun(@(a) atan(beta/a),alpha).';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Begin for loops
for indexOmega = 1:numEl
    tic
    thisOmega = 1i*omega(indexOmega);
    thisGamma = mcgama(1/2,thisOmega);
    cosarg1 = 1/2 + thisOmega;
    denomExpArg = 1/4 + thisOmega/2;
    
    
    for indexBeta = 1:numEl
        
        thisBeta = beta(indexBeta);
        thisBetaSquared = thisBeta^2;
        
        for indexAlpha= 1:numEl
            thisAlpha = alpha(indexAlpha);
            thisAlphaSquared = thisAlpha^2;
            thisBoverAaTans = atan(thisBeta/thisAlpha);
            
            Spectrum(:,indexAlpha,indexBeta,indexOmega) = (abs((thisGamma*cos(cosarg1*thisBoverAaTans))/((thisAlphaSquared+thisBetaSquared)^(denomExpArg))))^2;
            
        end
        
    end
    toc
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Eric Diaz wrote:

> In addition, I fixed up my code because I wasn't getting results
> consistent with what I wanted. And, actually, it turns out that it's
> even 10x faster with the inside for loop. Not sure why, but it is.
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> close all force
> clear all force
> clc
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> % Initialize variables
> numEl = 40;
> Spectrum = zeros(numEl,numEl,numEl,numEl);
> beta = linspace(-1 ,1 ,numEl);
> omega = linspace(0*pi ,2*pi ,numEl);
> alpha = linspace(-1 ,1 ,numEl);
> % theseBoverAaTans = arrayfun(@(a) atan(beta/a),alpha).';
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> % Begin for loops
> for indexOmega = 1:numEl
> tic
> thisOmega = 1i*omega(indexOmega);
> thisGamma = mcgama(1/2,thisOmega);
> cosarg1 = 1/2 + thisOmega;
> denomExpArg = 1/4 + thisOmega/2;
> for indexBeta = 1:numEl
> thisBeta = beta(indexBeta);
> thisBetaSquared = thisBeta^2;
> for indexAlpha= 1:numEl
> thisAlpha = alpha(indexAlpha);
> thisAlphaSquared = thisAlpha^2;
> thisBoverAaTans = atan(thisBeta/thisAlpha);
> Spectrum(:,indexAlpha,indexBeta,indexOmega) =
> (abs((thisGamma*cos(cosarg1*thisBoverAaTans))/((thisAlphaSquared+thisBetaSquared)^(denomExpArg))))^2;

If I have traced properly, all of the elements on the right hand side of
the assignment there are scalar variables, so the result will be a
scalar value, but you are assigning it to a vector. That's going to copy
the scalar value to each element of that vector. I think you should be
dropping the first dimension of the Spectrum array, and dropping the :
index there, if you are going to loop explicitly.

Note by the way that you can pre-calculate alpha.^2 outside of all of
the loops, and then just index that array instead of calculating
thisAlphaSquared inside the inner loop.
>
> end
> end toc
> end
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%