Message Boards Message Boards

GROUPS:

Performance of recursion code with radial polynomials

Posted 2 months ago
433 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

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 2 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.

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 2 months ago

I mean [0,1] is 0<=r<=1.

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.

  1. 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 2 months ago

the code is impressive. Thanks.

Posted 2 months ago

May I ask what is the further step to optimize this code so that the runtime will be less than 1 sec?

I think the code above by Neil can be speeded up if you change Table to ParallelTable.

Posted 2 months ago

I am afraid I cannot agree with you. It will spend more time

@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 2 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 2 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-1ST

The 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]

enter image description here

Posted 2 months ago

thanks

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

Group Abstract Group Abstract