Group Abstract Group Abstract

Message Boards Message Boards

0
|
8.6K Views
|
14 Replies
|
3 Total Likes
View groups...
Share
Share this post:

Performance of recursion code with radial polynomials

Posted 7 years ago

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.

POSTED BY: look xu
14 Replies
Posted 7 years 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

POSTED BY: look xu
POSTED BY: EDITORIAL BOARD
Posted 7 years ago

thanks

POSTED BY: look xu

@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 BY: Neil Singer
Posted 7 years ago
POSTED BY: look xu
POSTED BY: Neil Singer
Posted 7 years ago
POSTED BY: look xu
Posted 7 years ago

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

POSTED BY: look xu

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

POSTED BY: Sompob Saralamba
Posted 7 years ago
POSTED BY: look xu

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 BY: Neil Singer
Posted 7 years ago

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

POSTED BY: look xu

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 BY: Neil Singer
Posted 7 years 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 BY: look xu
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard