Group Abstract Group Abstract

Message Boards Message Boards

Slow Summation of values produced by recursive function

Posted 12 years ago
POSTED BY: Tom Zinger
13 Replies
POSTED BY: Daniel Lichtblau
Posted 10 years ago

this is really impressive thank you. It its too advanced for me, so I have to study those functions and methodology of yours, nevertheless thank you, I appreciate it.

POSTED BY: Tom Zinger
Posted 10 years ago

I am still trying to understand you code above, and I really have difficulties. @Daniel Lichtblau

I would appreciate it, if you could spend 5 minutes to explain it in more detail. If it is much trouble though, then its ok, no need to devote more time. Anyhow, thanks again

POSTED BY: Tom Zinger

The short answer is that as this was a year ago I have no recollection of the problem, let alone why I saw fit to code it as I did. The next best thing I can do is show intermediate results from running a smaller example. This might make it easier to see what happens at each step.

It also might be helpful to look at what ListConvolve does. Here is a short example.

ListConvolve[{a, b, c, d}, {1, 2, 3, 4, 5, 6}]

(* Out[682]= {4 a + 3 b + 2 c + d, 5 a + 4 b + 3 c + 2 d, 
 6 a + 5 b + 4 c + 3 d} *)

In effect it lines up the first vector under the second one, multiplies elementwise and sums, then slides the first vector one place to the right, repeats, and continues in this way until they are aligned on their right ends (they start out aligned on the left ends).

Also have a look at documentation for Nest. It is a neat way to iterate a function application without forming an explicit looping construct.

Here is the code with some Print statements in place for showing intermediate results.

gSum[t_, d_] := 
 Module[{iter = Ceiling[101 - t], vec, row}, 
  row = ConstantArray[0, 2*d + 1 + 10*iter];
  row[[d + 5*iter + 1]] = 1;
  vec = Reverse@{pa1, 0, pa2, 0, pa3, p0, pb3, 0, pb2, 0, pb1};
  Print[row];
  Print[vec];
  Range[-d, d].Nest[(Print[#]; ListConvolve[vec, #]) &, row, 
    100 - t + 1]]

Now we try it out on a "small" example, where we decrease the number of iterations and also d.

gSum[90, 2]

(*During evaluation of In[677]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}

During evaluation of In[677]:= {0.0187268926056,0,0.0112361355634,0,0.00216035798122,0.926427083333,0.00278755868545,0,0.0144982394366,0,0.0241637323944}

During evaluation of In[677]:= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}

During evaluation of In[677]:= {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.0187268926056,0.,0.0112361355634,0.,0.00216035798122,0.926427083333,0.00278755868545,0.,0.0144982394366,0.,0.0241637323944,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}

During evaluation of In[677]:= {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.000350696506663,0.,0.000420835807996,0.,0.000207164326207,0.0346982009931,0.000152952774553,0.0208189205958,0.000610323866969,0.004002828287,0.859510016593,0.00516493972516,0.000613427203787,0.0268631233495,0.000185234010799,0.0447718722491,0.000344914590979,0.,0.000700663155872,0.,0.000583885963227,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}

During evaluation of In[677]:= {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,6.56745581745*10^-6,0.,0.0000118214204714,0.,9.36574227431*10^-6,0.000974684225409,7.07879972564*10^-6,0.00116962107049,0.0000198532031615,0.000575767927495,0.0482636460291,0.000425098778452,0.0289711879004,0.00169626167989,0.00558446399082,0.798576225484,0.00720328490394,0.00170488672572,0.0373819616221,0.000514817413076,0.0622758072744,0.000958614655559,0.0000265440325772,0.00194733997168,0.000012691141759,0.00162278330973,0.0000201204224617,0.,0.0000253959554956,0.,0.0000141088641642,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}

During evaluation of In[677]:= {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.22988039786*10^-7,0.,2.95171295485*10^-7,0.,3.22406388303*10^-7,0.0000243370757515,2.81644340655*10^-7,0.0000438067363527,5.99729963139*10^-7,0.0000347067091938,0.00180739653767,0.0000262319671333,0.00216888988998,0.000073570180399,0.00106816189014,0.0597296922676,0.000789050728104,0.0358859906172,0.00314548877856,0.00695241079309,0.743026294279,0.00896168069601,0.00316154466301,0.046303702426,0.000955595604231,0.0770710708598,0.00177829166283,0.0000983644427217,0.00361104884579,0.0000470296697758,0.00300921960408,0.0000745604171864,1.09488243574*10^-6,0.0000941100039131,6.99648969278*10^-7,0.0000522833355073,8.93710434366*10^-7,0.,8.18214763327*10^-7,0.,3.40922818053*10^-7,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}

During evaluation of In[677]:= {0.,0.,0.,0.,0.,0.,0.,2.30318381284*10^-9,0.,6.90955143853*10^-9,0.,9.61995269296*10^-9,5.69697254917*10^-7,9.87743724844*10^-9,1.3672734118*10^-6,1.76981031756*10^-8,1.49343004982*10^-6,0.0000564089441028,1.30461472525*10^-6,0.000101522355044,2.77803040269*10^-6,0.000080444698405,0.00279517589237,0.0000608197879561,0.00335427612817,0.000170500598533,0.00165345256966,0.0693649220775,0.00122266773579,0.0417305935458,0.00486469222763,0.00814554217797,0.692325575669,0.0104891127162,0.0048897146358,0.0538440661273,0.00148077578157,0.0895042858466,0.00275230472928,0.000227967410561,0.00558461451971,0.000109034821235,0.00465390209609,0.000172814143821,5.07164370766*10^-6,0.000218100788761,3.24086876983*10^-6,0.000121186008162,4.13978775527*10^-6,4.46896865368*10^-8,3.79008158365*10^-6,3.28806953291*10^-8,1.57920065985*10^-6,3.44083956859*10^-8,0.,2.47139032277*10^-8,0.,8.23796774257*10^-9,0.,0.,0.,0.,0.,0.,0.}

During evaluation of In[677]:= {0.,2.6276417913*10^-10,1.28023911727*10^-8,3.1441216378*10^-10,3.8407173518*10^-8,5.15850133524*10^-10,5.34731082909*10^-8,1.58455004758*10^-6,5.49043522853*10^-8,3.80206810302*10^-6,9.8376012633*10^-8,4.15304072344*10^-6,0.000104675517164,3.62854672792*10^-6,0.000188339532114,7.72501528272*10^-6,0.000149280603238,0.00389362244643,0.000112930238318,0.0046725276082,0.000316309549488,0.0023063885088,0.0774047329677,0.00170811881128,0.0466499269642,0.00677668471315,0.00919583452073,0.645996802126,0.0118261138072,0.00681193987006,0.0601899893402,0.00206878930601,0.0998792329999,0.00383836226103,0.000422941409518,0.00777936345274,0.000202436666432,0.00648296101122,0.00032067065861,0.0000141032691988,0.000404611625896,9.01350947964*10^-6,0.000224889311753,0.0000115120471474,2.4841041572*10^-7,0.0000105393921256,1.8276940003*10^-7,4.39256640586*10^-6,1.91261217945*10^-7,1.79845104264*10^-9,1.3737377571*10^-7,1.38006989711*10^-9,4.579125857*10^-8,1.21270717041*10^-9,0.}

During evaluation of In[677]:= {3.34528274247*10^-9,1.73468577083*10^-7,3.43045647599*10^-6,1.78126343116*10^-7,8.22757025872*10^-6,3.19130770453*10^-7,8.98775868828*10^-6,0.000170089310524,7.85517643323*10^-6,0.000305912923289,0.0000167165526912,0.000242576350162,0.00506618154608,0.000183671994521,0.00607980780256,0.000513781540176,0.00300642321839,0.0840554447009,0.00223110506019,0.0507695757096,0.00881793681269,0.0101295224527,0.603613575418,0.0130061526652,0.00886450023817,0.065503472249,0.00270234721934,0.108462156308,0.00500198737982,0.000687033480414,0.0101222973509,0.000329200186709,0.00843559244134,0.000521033179746,0.0000305200018317,0.000657198653216,0.0000195110786687,0.00036545123372,0.0000249128926341,8.05850296226*10^-7,0.0000228072133743,5.92941711153*10^-7,9.51048289645*10^-6,6.20452868253*10^-7,1.16629362777*10^-8}

During evaluation of In[677]:= {7.89130328169*10^-7,0.000016681600964,0.000252881212855,0.0000145864210549,0.000454573353685,0.000031022471229,0.000360666238107,0.0062829361827,0.000273409964517,0.00754024424973,0.000763480107174,0.00373692672535,0.0894978915393,0.0027802147369,0.0541983663809,0.0109364682763,0.010968068657,0.564794626503,0.0140568028427,0.0109952836409,0.0699249096885,0.00336765739379,0.115486402647,0.00621522393068,0.00102103123626,0.0125536771221,0.000489947480649,0.0104620440197,0.000774587837456,0.0000566422181736,0.000976574614961,0.0000362260046162,0.000543384380967,0.0000462370474296,1.99271236802*10^-6}

During evaluation of In[677]:= {0.0000518414577311,0.000503109638511,0.00751957532076,0.000381953707819,0.0090246886675,0.00106428669011,0.00448452289808,0.0938905651949,0.00334638939414,0.0570306069733,0.0130900472949,0.0117289770542,0.529198751593,0.0150007259605,0.0131619696786,0.0735759755613,0.00405377678648,0.121156480681,0.00745554634111,0.00142348260056,0.0150249590323,0.000684296304377,0.012521894581,0.00108034722126,0.0000946620041171}

During evaluation of In[677]:= {0.00141383871961,0.00523853304996,0.0973723655727,0.00392248282427,0.0593480258782,0.015244653906,0.0124264739085,0.496520342016,0.015856497841,0.0153304809202,0.0765621207482,0.00475207649088,0.12565161287,0.0087049421142,0.00189128026671}

Out[677]= 0.00376821139042*)
POSTED BY: Daniel Lichtblau
Posted 10 years ago

Thank you, you are right. i really appreciate your help , i will go through what you suggest. Your clarifications are really helpfull. @Daniel Lichtblau

POSTED BY: Tom Zinger
Posted 12 years ago

I'm sorry. That was not attitute.

You said "unfortunately I cannot provide community with actual code" even after several different people have specifically asked you for that again and again in several different ways.

You then provided one fragment, but said "I have more functions like f[t,d], which means that my overall speed is significantly lower.", It was my mistake when I wrote "show the function f", when I should have written, as others had, "show all the code." It seems it still isn't possible to understand where and why the code is slow or have a chance of providing your "superfast" solution in a single try.

I sincerely hope it works out for you and someone else can provide what you want.

POSTED BY: Bill Simpson
Posted 12 years ago

Thanks a lot for your assistance!

@Daniel Lichtblau : I will try to implement what you said, thanks. With respect to excel, maybe I was not clear. My original algorithm didn't involve recursive functions. I did it in excel, where I created a huge spreadsheet where each state t (0 until 100) was a column of cells whose values where calculated based on previous step. So it was really big and painful. But also really fast since it didn't involve any function calling it self. But it was really hard to maintain it, while with the recursive function f (and mathematica), it was easier for me. The slow speed apparently has to do with me not writing "correct" code. Nevertheless I appreciate your support.

@Bill Simpson I really do not understand your attitude. f is provided already, so why the irony. Nevertheless thank you too for trying to help in your own way. I appreciate it.

POSTED BY: Tom Zinger
Posted 12 years ago

Since you know a great deal about your function f, but this function must remain a secret, take a few minutes, hours or days and construct a different f, one that is different enough that your secret f will not be compromised, but which has enough similarities to the secret f that others are able to test this, see exactly the time behavior you are seeing and perhaps be able to to give you a some concrete advice on what to do.

Until then the problem is in the class "I typed some stuff in, some stuff came out, what do I do?"

POSTED BY: Bill Simpson

One way to possibly improve speed would be to recast as a double loop, iterating for t from 100 to 0 and d from something like -600 to 600 when t is 100, decreasing this range as t decreases since fewer such values will be needed. It could be done with Compile.

That said, I have to wonder what Excel might be doing to get a result instantaneously. These loops imply a complexity of O(n^2) at a minimum, where n is 100 in this case. Not huge by any means, but also not negligeable. Do you have the Excel code available?

POSTED BY: Daniel Lichtblau
Posted 12 years ago

Thanks a lot, unfortunately as you point out code runs faster only after the initial run.. unfortunately the implementation that is going to take place (involving java and mathlnk) demands that every time code is executed kernel is quit and all memory is freed. Hence, I have to rum it every time from scratch... imagine that the final product will be dynamic, and variables a1,a2...,t will change maybe 2-3 times per minute. This is why i need to improve speed. Total 4,5 seconds in my pc (i7 3720qm) is not that good... :(

No, no compile..i do not know where should I implement it? Do you mean compile function f?

POSTED BY: Tom Zinger

Hi,

using parallel made it a bit faster for me. On my machine this is faster than the original (at least after the Kernels are launched).

CloseKernels[]; 
ClearAll["Global`*"]
    a1 = 0;
    a2 = 0;
    a3 = 0;
    a4 = 0;
    b1 = 0;
    b2 = 0;
    b3 = 0;
    b4 = 0;

    a = 20.50;
    b = 15.;
    pa1 = 0.0241637323943662`;
    pa2 = 0.014498239436619723`;
    pa3 = 0.0027875586854460006`;
    pb1 = 0.018726892605633805`;
    pb2 = 0.011236135563380286`;
    pb3 = 0.00216035798122065`;
    p0 = 0.9264270833333333`;

    aTotal = 5*a1 + 3*a2 + 1*(a3 + a4);
    bTotal = 5*b1 + 3*b2 + 1*(b3 + b4);
    difference = aTotal - bTotal;
    t = 0;
    taf = 100 - t;

    pa1 = 0.0241637323943662`;
    pa2 = 0.014498239436619723`;
    pa3 = 0.0027875586854460006`;
    pb1 = 0.018726892605633805`;
    pb2 = 0.011236135563380286`;
    pb3 = 0.00216035798122065`;
    p0 = 0.9264270833333333`;

    f[t_, d_] := 
      f[t, d] = 
       If[t > 100, 0, 
        If[d == difference && t == taf, 1, 
         pa1*f[t + 1, d - 5] + pa2*f[t + 1, d - 3] + pa3*f[t + 1, d - 1] +
           pb1*f[t + 1, d + 5] + pb2*f[t + 1, d + 3] + 
          pb3*f[t + 1, d + 1] + p0*f[t + 1, d]]];
    Quiet[LaunchKernels[]];
    InTime = AbsoluteTime[];
    (*expected=Sum[j*f[0,j],{j,-100,100}]*)
    expected = 
     Total[ParallelTable[j*f[0, j], {j, -100, 100}]]
    OutTime = AbsoluteTime[];
    OutTime - InTime

But there is a problem because the f function is not properly cleaned. If you rerun the program with different parameters it gives the same results unless you properly clean the function f. Mathematica obviously uses the assignments made in the first run and does not clear them even if you use ClearAll. If you clear the function and relaunch the Kernels etc. this method does not seem to help speed everything up. Perhaps it is possible to use Parallel somewhere outside of the larger structure that you describe in your code.

Have you tried compiling it?

Best wishes, M.

POSTED BY: Marco Thiel
Posted 12 years ago

hello, thanks for replying. I will add a part of my code to illustrate the problem:

ClearAll["Global`*"]

a1 = 0;
a2 = 0;
a3 = 0;
a4 = 0;
b1 = 0;
b2 = 0;
b3 = 0;
b4 = 0;

a = 20.50;
b = 15.;
pa1 = 0.0241637323943662`;
pa2 = 0.014498239436619723`;
pa3 = 0.0027875586854460006`;
pb1 = 0.018726892605633805`;
pb2 = 0.011236135563380286`;
pb3 = 0.00216035798122065`;
p0 = 0.9264270833333333`;

aTotal = 5*a1 + 3*a2 + 1*(a3 + a4);
bTotal = 5*b1 + 3*b2 + 1*(b3 + b4);
difference = aTotal - bTotal;
t = 0;
taf = 100 - t;

pa1 = 0.0241637323943662`;
pa2 = 0.014498239436619723`;
pa3 = 0.0027875586854460006`;
pb1 = 0.018726892605633805`;
pb2 = 0.011236135563380286`;
pb3 = 0.00216035798122065`;
p0 = 0.9264270833333333`;

f[t_, d_] := 
  f[t, d] = 
   If[t > 100, 0, 
    If[d == difference && t == taf, 1, 
     pa1*f[t + 1, d - 5] + pa2*f[t + 1, d - 3] + pa3*f[t + 1, d - 1] +
       pb1*f[t + 1, d + 5] + pb2*f[t + 1, d + 3] + 
      pb3*f[t + 1, d + 1] + p0*f[t + 1, d]]];

InTime = AbsoluteTime[];
expected = Sum[j*f[0, j], {j, -100, 100}]
OutTime = AbsoluteTime[];
OutTime - InTime

Now, I have more functions like f[t,d], which means that my overall speed is significantly lower. I understand that working with the very same algorithm but converted to serial (done it in excel) runs instantly. Is there any way I can speed up in mathematica while forking on such recursive functions?

Thank you in advance.

POSTED BY: Tom Zinger

Hard to assess without knowing the function f[t,x]

POSTED BY: Daniel Lichtblau
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard