9

I have two functions, one of which calls the other. I need to do some calculations but I was unable to optimize my code to do it efficiently.

In search I came across the memoization technique, which improved my code significantly, but still not enough.

Below is my code:

Clear[function1,function2]

function1[i_,j_]:=function1[i,j]=
  Which[
      i >= 40, 1.0,
      j >= 40, 0.0,
      FractionalPart[j] >= 0.40, 0.0,
      True, 0.18*0.45*(1 - function2[j, i+1]) + 0.66*0.45*(1 - function2[j, i+2]) + 0.16*0.45*(1 - function2[j, i+3]) + (1 - 0.18*0.45 - 0.66*0.45 - 0.16*0.45)*(1 - function2[j, i])
  ];


function2[j_,i_]:=function2[j,i]=
  Which[
    j >= 40, 1.0,
    i >= 40, 0.0,
    True, 0.18*0.55*(1 - function1[i, j+1]) + 0.66*0.55*(1 - function1[i, j+2]) + 0.16*0.55*(1 - function1[i, j+3]) + (1 - 0.18*0.55 - 0.66*0.55 - 0.16*0.55)*(1 - function1[i, j+0.01])
  ];

And my timings:

function1[0, 0]//AbsoluteTiming
(* Out: {8.56012,0.218588} *)

Any idea how to improve speed? I am positive that my approach is just too amateur, as I am not experienced with Mathematica. Thank you in advance.

Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
Frank Muller
  • 93
  • 3
  • 1
    You might think you could speed it up by compiling, but you would be wrong - just goes and goes without answering – Jason B. Mar 01 '16 at 11:30
  • @JasonB Currently there seems to be no way to compile this type of function. Some related posts: http://mathematica.stackexchange.com/questions/108121/how-to-speed-up-calculation-of-this-equation-findroot/108123#comment293404_108121 – xzczd Mar 01 '16 at 11:44
  • I have the suspicion that this problem can be vectorized but I cannot seem to wrap my mind around it as they say. If that can be done it will likely be a large improvement. – Mr.Wizard Mar 01 '16 at 14:26

2 Answers2

9

I found that it is considerably faster, if you process the fractional part of j as a seperate integer variable:

Clear[function1, function2]

function1[i_, j_, jc_] := function1[i, j, jc] = 
  Which[i >= 40, 1.0, j >= 40, 0.0, jc >= 40, 0.0, True, 
    0.18*0.45*(1 - function2[j, i + 1, jc]) + 
    0.66*0.45*(1 - function2[j, i + 2, jc]) + 
    0.16*0.45*(1 - function2[j, i + 3, jc]) + (1 - 0.18*0.45 - 
    0.66*0.45 - 0.16*0.45)*(1 - function2[j, i, jc])];


function2[j_, i_, jc_] := function2[j, i, jc] = 
  Which[j >= 40, 1.0, i >= 40, 0.0, True, 
    0.18*0.55*(1 - function1[i, j + 1, jc]) + 
    0.66*0.55*(1 - function1[i, j + 2, jc]) + 
    0.16*0.55*(1 - function1[i, j + 3, jc]) + (1 - 0.18*0.55 - 
    0.66*0.55 - 0.16*0.55)*(1 - function1[i, j, jc + 1])];
Berg
  • 1,059
  • 5
  • 14
  • It can be sped up for about 25% more if you FullSimplify the expressions for the True condition. – shrx Mar 01 '16 at 11:41
  • It' interesting that the two methods disagree by a small amount, about 10^-6 on my machine, also there's some difference in the memoized results. If you run OP's example, then ask for function1[33, 38.39], you get 0.00225468. But if you run your code, starting with function1[0, 0, 0] then the memoized result for function1[33, 38, 39] is 0. – Jason B. Mar 01 '16 at 11:57
  • 1
    This is an inconsistency in the OPs code: If you run only function1[33,38.39] it returns 0. If you run function1[0,0] first it returns 0.00225468. In my code 0 is always returned. Presumably a rounding error in j. – Berg Mar 01 '16 at 12:14
  • Thanks for your reply, i confirm the speed improvement. in my case, I will always call for intenger values (i,j) where results seem to be the same. – Frank Muller Mar 01 '16 at 13:18
  • @FrankMuller - but in fact you don't always call with integer values. When you call function1[0,0], it recursively calls function1 and function2 a few hundred thousand times, giving decimal values as the arguments. So the difference is relevant. I agree with Berg that it appears to be the result of a cummulative rounding error – Jason B. Mar 01 '16 at 13:32
  • @JasonB I mean I ask for function1[0,0] not function1[0,0.1] for example. I understand what you mean, some decimals take place in order to stop endless recursions. Hence I used fractionalpart in original function. But I will always ask for function1[n,m] where n,m Intengers – Frank Muller Mar 01 '16 at 13:35
6

The primary inefficiency of your code comes from quite a lot of unnecessary evaluations and memoizations because of pattern matching on machine numbers. If we eliminate these the function is several times faster. Simply changing function1[i, j+0.01] in your last line of code to function1[i, j + 1/100] is sufficient.

With your original:

function1[0, 0] // AbsoluteTiming
function1 // DownValues // Length
function2 // DownValues // Length

{6.90192, 0.218588}

342540


323897

With that sole change:

function1[0, 0] // AbsoluteTiming
function1 // DownValues // Length
function2 // DownValues // Length

{1.48465, 0.218587}

69036


67436

Berg's code fundamentally works in the same way.

Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
  • This is really helpful, thanks. But I still cannot understand why it works. Why expressing 0.01 as 1/100 makes such a difference in skipping unnessesary calculations. I am trying to implement this logic in other calculations with similar recursive functions, but I gain nothing. Is there a documentation explaining that? Thank you in advance. – Frank Muller Mar 10 '16 at 13:46
  • @FrankMuller The issue is that the inexact arithmetic may produce slightly different values for what "should" be the same number, and pattern matching sees these as different, causing the calculation and memoization to needlessly be performed for each. It is possible that other operations will not produce slightly different values and you will not see the same improvement, or that certain operations are sufficiently slower on exact numbers that you lose the advantage. I still have the suspicion that there is a better way to handle this particular problem but failed to find it. – Mr.Wizard Mar 11 '16 at 16:02