# Performance of recursion code with radial polynomials

Posted 8 months ago
1121 Views
|
14 Replies
|
3 Total Likes
|
 Consider the following code: Clear[r, re, p, pmax, delta, imagesize, delta] ClearSystemCache[] re[0, r_] := Sqrt[8/Pi]*((1 - r)/r)^(1/4)*1; re[1, r_] := Sqrt[8/Pi]*((1 - r)/r)^(1/4)*-1*2*(1 - 2*r); re[p_, r_] := re[p, r] = Sqrt[8/Pi]*((1 - r)/r)^(1/4)*(-1)^p*(re[1, r]*re[p - 1, r] - re[p - 2, r]); imagesize = 32; pmax = 10; delta = 2/imagesize; Table[r = Sqrt[x^2 + y^2]; re[pmax, r], {x, -1 + delta/2, 1 - delta/2, delta}, {y, 1 - delta/2, -1 + delta/2, -delta}]; when the imagesize and pmax increase, the time will become unacceptable. So, I would ask if we can use compile of other methods to speed up,like: for imagesize is 256 and pmax is 120, the time will be about 10 seconds. In my code, I also use the memoization to store the value during the evaluation which I will use in the future.
14 Replies
Sort By:
Posted 8 months ago
 Your first problem is that you are not computing the expressions numerically. They have exact integer values so MMA is keeping these huge expressions in memory. Add N[] around your expressions so they become numerical. This will give a huge speed up. Second, your results are coming out complex. Is this correct? I have some other suggestions but I wanted to verify if the expressions are correct before diving deeper.Regards,Neil
Posted 8 months ago
 Sorry, I forgot to define the r, r should be [0,1]. Or the result would be complex which I do not need.
Posted 8 months ago
 You have a bug. r is defined to be Sqrt[x^2+y^2] but also 0.1? I think you need to fix some notation. I also would not use r to be the pattern variable in your function definitions. (Three uses of the same variable name can cause strange results). Can you fix the code and repost it in complete form?Regards,Neil
Posted 8 months ago
 I mean [0,1] is 0<=r<=1.
Posted 8 months ago
 Your original code never completes (I killed it after 5 minutes). I had to protect the evaluation of re[] for values of r that make it complex by adding an if statement. If[r <= 1, re[pmax, r], 0] If you make the table limits floating point, (delta/2. vs delta/2) it takes about 24.5 seconds. You have several issues. You recompute Sqrt[8/Pi] many times -- do it once. 2.This code takes 4.9 seconds or a speedup of 5 times faster. Clear[r,re,p,pmax,delta,imagesize,delta] ClearSystemCache[]; sqrt8 = Sqrt[8/Pi]; cf1=Compile[{{rr,_Real}},sqrt8*((1-rr)/rr)^(1/4)*1]; re[0,r_]:=cf1[r]; cf2=Compile[{{rr,_Real}},sqrt8*((1-rr)/rr)^(1/4)*-1*2*(1-2*rr)]; re[1,r_]:=cf2[r]; cf3 = Compile[{{pp,_Integer},{rr,_Real}},sqrt8*((1-rr)/rr)^(1/4)*(-1)^pp]; re[p_,r_]:=re[p,r]=cf3[p,r]*(re[1,r]*re[p-1,r]-re[p-2,r]); imagesize=256; pmax=120; delta=2/imagesize; {foo,bar}=Timing[Table[r=Sqrt[x^2+y^2];If[r<=1,re[pmax,r],0],{x,-1+delta/2.,1-delta/2.,delta},{y,1-delta/2.,-1+delta/2.,-delta}]];  There are some further speedups like maybe not looping through values of r that make no sense but this is a start.Regards,Neil
Posted 8 months ago
 the code is impressive. Thanks.
Posted 8 months ago
 May I ask what is the further step to optimize this code so that the runtime will be less than 1 sec?
Posted 8 months ago
 I think the code above by Neil can be speeded up if you change Table to ParallelTable.
Posted 8 months ago
 I am afraid I cannot agree with you. It will spend more time
Posted 8 months ago
 @Sompob Saralamba That is a good suggestion, however the way the code is currently setup various kernels would need to calculate the same variables over again because of the recursion. Is there a way for the kernels to share variables? I have not looked into that topic recently.My latest thought is to restructure the calculation to see if that order can be changed to take advantage of the Parallelism or reduce computation.My suggestion:It would be faster if you can do the calculation in variable "r" as opposed to "x" and y". If you look at the result it only depends on r and it is symmetric. In my posted version of your calculation I am computing many zeros and many duplicate values. I would try to do the calculation for unique values of r and not, x,y. It also depends a bit on what form you need the answer in at the end but I would try that (you might not need any parallelization if this approach works because there are not many unique values. To prove this point, you compute 256*256 = 65,536 values Do cnt = Counts[Flatten[bar]] From this you see that it computes zero 14,068 times! If you do Keys[cnt] // Length you get only 4006 unique values other than 0. This means you are doing 65536 iterations to obtain 4006 unique values. With some cleverness you should be able to directly compute the 4000 values you actually need.Regards,Neil
Posted 8 months ago
 Also，I think my code is not right at the beginning, cause It dose calucate this part "Sqrt[8/Pi]((1 - r)/r)^(1/4)(-1)^p" in every recursion which is not right. it is my fault. But still good code I can learn from. Thanks.
Posted 8 months ago
 re = Compile[{{pmax, _Integer}, {r, _Real}}, Block[{C1, Cold, Cp, buffer}, If[0. <= r <= 1., If[pmax == 0, Cp = 1., C1 = 2. (1. - 2. r); Cold = 1.; Cp = C1; Do[buffer = Cp; Cp = C1 Cp - Cold; Cold = buffer, {p, 2, pmax}]; ]; Sqrt[8./Pi] ((1 - r)/r)^(1/4) (-1.)^pmax Cp, 0. ] ], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True] imagesize = 256; delta = 2./imagesize; x = Subdivide[-1. + delta/2, 1. - delta/2, imagesize - 1]; y = Subdivide[-1. + delta/2, 1. - delta/2, imagesize - 1]; result = re[120, Sqrt[0.5 Outer[Plus, x^2, y^2]]]; // AbsoluteTiming // First This is another answer from stackexchange.com
 @look xu - Please make sure you know the rules: https://wolfr.am/READ-1STThe rules explain how to format your code properly. If you do not format code, it may become corrupted and useless to other members. Please EDIT your post and make sure code blocks start on a new paragraph and look framed and colored like this. int = Integrate[1/(x^3 - 1), x]; Map[Framed, int, Infinity]