# Slow Summation of values produced by recursive function

Posted 9 years ago
16228 Views
|
13 Replies
|
6 Total Likes
|
 Hi, I have defined a recursive function using memoization f[t_,x_]:=f[t,x]= ...... then i need to execute some sums like this marg = -Sum[j*f[0, j], {j, -100, 100}]; this sum takes up to 0.85 sec in order to be executed, depending on the function (actually I have more recursive functions, and i use one sum like the previous one for each. Hence, my algorithm speed overall is too slow, resulting in 4,50 seconds about).I tried two alternative methods to improve speed but unsuccesfully1) using tables, and then function total, like this pinakas = Table[i*f[0, i], {i, 0, 100}] Total[pinakas] 2) using for function instead of Sum s2=0 For[ i=0,i<=100,i++, s2=s2+f[0,?] ] All methods (i guess thats logical) yield the same times more or less. I have setup my algorithm in an excel file, and while serial, my algorithm is superfast. I cannot understand why, but if I print all values needed from my function, print is instant (which means that each f[0,i] is really fast calculated) but this sum is so slow. I though that It would be wise to limit precision somehow in my f[0,i] values, since I only need 3-4 digits, but I do not know how to exactly perform that. I even used decimals instead of integers (1.0 instead of 1 for example) as suggested in various guides, but this resulted in a decrease of speed.So here I am, asking for your assistance.. unfortunately I cannot provide community with actual code, but I would appreciate any tip or assistance.Regards
13 Replies
Sort By:
Posted 9 years ago
 With your definitions I note that difference is 0 and taf is 100. Running the summation gives a result of 3.75977112671 in around one second.We can instead define this as an iteration from t=100 down to 0, and construct an array of appropriate values at each step. This in turn can be done explicitly or, as I'll indicate using vector convolution.(Below is a revised version of the original. Slightly cleaner, using Nest rather than Do). 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}; Range[-d, d].Nest[ ListConvolve[vec, #] &, row, 100 - t + 1] ] In[292]:= gSum[0, 100] // Timing Out[292]={0.004000, 3.75977112671} So that's a speed gain of two+ orders of magnitude.
Posted 8 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 8 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 8 years ago
 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 8 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 9 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 9 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 9 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 9 years ago
 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 9 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 9 years ago
 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 9 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 9 years ago
 Hard to assess without knowing the function f[t,x]
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.