# Improve speed in this tree - memoization won't suffice

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.

-
You might think you could speed it up by compiling, but you would be wrong - just goes and goes without answering – JasonB Mar 1 at 11:30
@JasonB Currently there seems to be no way to compile this type of function. Some related posts: mathematica.stackexchange.com/questions/108121/… – xzczd Mar 1 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 1 at 14:26

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.

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.

-
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 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 at 16:02

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])];

-
It can be sped up for about 25% more if you FullSimplify the expressions for the True condition. – shrx Mar 1 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. – JasonB Mar 1 at 11:57
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 1 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 1 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 – JasonB Mar 1 at 13:32