How to speed up my code

I haven't been able to vectorize the code below. Below,

 lent=600, lens=950, K=400.

I have values for B, d1, d2, ExpSum2, d1star, and d2star.

The size of d1,d2,d1star,d2star are all 950x600x400.
The size of B is 1350x600.
The size of ExpSum2 is 400.
I have values for bs, and length(bs)=1550.

I have posted two sets of for loops below. They should do the same thing, but the first one ran faster for me.

I have also vectorized the multiplication B.*d2 and pulled it outside the for loop, but I couldn't vectorize

B(s+k,tt).*d1(s,tt,k).

I still had to run the for loop, and the run time increased. Right now it takes the program between 2 hours and 2hours and 20 minutes to run.

Is there a way to significantly speed this up? I was hoping I could vectorize this, but if not I may try a mex file. Thanks for your time.

for tt=1:lent
  for k=1:K 
    for s=1:lens  
      if bs(s)>.1
        J(s,tt,k)=0;
      else
        J(s,tt,k)=.99^k*(B(s+k,tt)*d1(s,tt,k)-B(s,tt)*d2(s,tt,k)*ExpSum2(k));
      end
      if J(s,tt,k)<0
        J(s,tt,k)=0;
      end
      if bs(s)>.1
        Je(s,tt,k)=0;
      else
        Je(s,tt,k)=.99^k*(B(s+k,tt)*d1star(s,tt,k)-B(s,tt)*d2star(s,tt,k)*ExpSum2(k));
      end
      if Je(s,tt,k)<0
        Je(s,tt,k)=0;
      end
    end
  end
end

for tt=1:lent
  for k=1:K 
    for s=1:lens 
      if bs(s)<=.1
        J(s,tt,k)=max(0,.99^k*(B(s+k,tt)*d1(s,tt,k)-B(s,tt)*d2(s,tt,k)*ExpSum2(k)));
        Je(s,tt,k)=max(0,.99^k*(B(s+k,tt)*d1star(s,tt,k)-B(s,tt)*d2star(s,tt,k)*ExpSum2(k)));
      end
    end
  end
end

[Adam]

J and JE don't appear to be being pre-sized unless you left that code off. Obviously if vectorisation is possible then pre-sizing will likely become irrelevant though.

[Brian Hogle]

Yes, I preallocated. Sorry, this is my first time asking a question, and I wasn't sure how much of the code to post. I'm also new to MATLAB. Thanks for reading.

[dpb]

If you preallocated, then J, JE will be zero everywhere they're not set to anything else, already so looks like you could rewrite the code like

if bs(s)>.1
  J(s,tt,k)=0;
else
  J(s,tt,k)=.99^k*(B(s+k,tt)*d1(s,tt,k)-B(s,tt)*d2(s,tt,k)*ExpSum2(k));
end

to

if bs(s)<=0.1
  J(s,tt,k)=.99^k*(B(s+k,tt)*d1(s,tt,k)-B(s,tt)*d2(s,tt,k)*ExpSum2(k));
end

that won't save a ton, but storing a zero on top of an existing one doesn't do anything useful... :)

After that, it appears that you could simply wait until done and then write

J(J<0)=0;

after the nested loop and save the if inside the loops as well; it doesn't appear that subsequent J are dependent on earlier values so leaving the <0 entry while computing doesn't appear to matter.

I've got other commitment; will try to take a look at the other Q?? re: vectorizing later...

[Stephen23]

"I wasn't sure how much of the code to post"

Post all of your code by uploading it as an attachment: click the paperclip button to do this.

[Brian Hogle]

Thanks, I've attached the full code. Dpb, you're right. You're comments make this much better. It looks cleaner too. Thank you.

[dpb]

I'd suggest kutting the sizes down and also attach a .mat file containing sufficient data to run a test case if expecting anybody to do more than just air-code/read what have.

Also, looks like a lot more code there than just this section; have you profiled to determine where the actual bottleneck is for certain?

[Brian Hogle]

Yes, I profiled. The big issue was the nested loops containing the J and Je. Also, the loops for d1, d2, d1star, and d2star take time (about 30 minutes) but the majority of the rest of the time is spent doing the calculation for J and Je. What I may end up doing is just using a smaller lens value.

I appreciate all your help.

Answer 1

[dpb]

Not surprising replacing with max increased run time...see if the following produces same result/makes any significant timing difference; you might want to cut the sizes down some for timing testing instead of the full-blown thing if is as long as say.

KK(1,1,:)=0.99.^[1:K];
for tt=1:lent
  for k=1:K 
    for s=1:lens  
      if bs(s)>=0.1
        J(s,tt,k)=B(s+k,tt)*d1(s,tt,k)-B(s,tt)*d2(s,tt,k);
      end
      if bs(s)<=0.1
        Je(s,tt,k)=B(s+k,tt)*d1star(s,tt,k)-B(s,tt)*d2star(s,tt,k);
      end
    end
  end
end
J=J.*KK.*ExpSum2;
JE=JE.*KK.*ExpSum2;
J(J<0)=0; Je(Je<0)=0;

ADDENDUM

Maybe my eyes fool me, but see if you can't replace the loop over s with

ix=bs>0.1;
J(ix,tt,k)=B(ix+k,tt)*d1(ix,tt,k)-B(ix,tt)*d2(ix,tt,k);
ix=~ix;
Je(ix,tt,k)=B(ix+k,tt)*d1star(ix,tt,k)-B(ix,tt)*d2star(ix,tt,k);

[Brian Hogle]

Dpb, I really like this answer. It only took a minute for when lens=200. So, I changed lens to 950, and it's been running for maybe 2 hours.

[dpb]

Did the refactoring make any difference of note on speed?

[Brian Hogle]

It made a big difference for smaller values of lens, but when lens=950, it ran a long time (I'd say a solid 2 hours), then I got and error:

"Array dimensions must match for binary array op.

Error in calcjrtKK (line 120) J=J.*KK.*ExpSum2;"

I will fix this.

[dpb]

Hmmm...which release are you running? That syntax requires a relatively new enhancement of auto-expansion. If you're on an earlier release may need bsxfun or explicit expansion of the one x/y position 3D elements.

I did check here on R2017b I had correct syntax but I forget which release was first introduced.

[Brian Hogle]

I'm using R2017b. I don't know. I'm confused because it ran successfully previously, but now it won't run anymore.

[dpb]

Let's see the specific code and error message in context.

The test I did was to create a little 3D array that mimics the J array and a K to multiply by with values that I could make sure got the expected answer:

>> K(1,1,:)=0.99.^[1:2];  % the K vector...
K(:,:,1) =
  0.9900
K(:,:,2) =
  0.9801
>> J          % the test J array...
J(:,:,1) =
   1     1
   1     1
J(:,:,2) =
   2     2
   2     2
>> J.*K             % check on implicit expansion and multiply
ans(:,:,1) =
  0.9900    0.9900
  0.9900    0.9900
ans(:,:,2) =
  1.9602    1.9602
  1.9602    1.9602
>>

That's what we expected...this is also R2017b here.

NB: It looks like I left one explicit step out above; you also have to turn the ExpSum2 vector into the version by third dimension--

ES(1,1,:)=ExpSum2;

so the implicit expansion occurs over the third dimension; perhaps you missed seeing that was needed and the dimensions don't matchup is the problem?