3

I am trying to optimize a Python code where the following mathematical function has to be evaluated many (many) times

from math import pow, pi
from scipy.special import gamma

def norm_spirals(a,k):
    return pi*pow(abs(gamma(0.25+0.5*k+0.5j*a)/gamma(0.75+0.5*k+0.5j*a)),2)

I have already optimized as much as possible my code using cProfile and timeit. I also got rid of the call to the procedure and direcly embedded this calculation in the code. The last step of optimization I could think of, is to tune the order of the mathematical operations in order to accelerate the evaluation. The formula above is the fastest form that I have been able to obtain using timeit.

  • Would you think of another order of calculation for this formula that could be faster ?
  • Would you have any ideas of other optimizations I could use ?

Thank you in advance for your help !

jibe
  • 198
  • 6
  • 1
    the next level would be to see how this function is called - i.e. in a for loop, in a list comprehension, in an inefficient while/next structure, etc. but sounds like you have already done that. – Corley Brigman Mar 07 '14 at 14:36
  • @CorleyBrigman This function is basically an integrand function called via `scipy.integrate.quad`. That's why it is evaluated so many times... – jibe Mar 07 '14 at 14:37
  • what kind of values do a and k take? might it be possible to store the results of this or gamma in some table (I note that the two calls to gamma are basically a difference of 1 in k), or are they different all the time? – RemcoGerlich Mar 07 '14 at 14:48
  • @RemcoGerlich `k` is an integer (positive or negative), whereas `a` is a float (on which I will integrate). The difference between the two arguments of the `gamma` function is exactly 1/2, so that there is no straigthforward mathematical simplification ([Wikipedia](http://en.wikipedia.org/wiki/Gamma_function)). I had thought of storing a table of the values of `norm_spirals`, but the problem is that I can't predict where the adaptative algorithm of `scipy.integrate.quad` will ask me to evaluate the function... – jibe Mar 07 '14 at 14:55
  • You may be sure that code is your only bottleneck, but if I were you, I would just spend a little effort to prove it with [*this technique*](http://stackoverflow.com/a/4299378/23771). Also, you are calling `pow(...,2)` and relying on it to be smart enough not to call `log` and `exp`. I wouldn't take that chance - I would just put `...` in a local and square it myself. – Mike Dunlavey Mar 07 '14 at 16:52
  • "there is no straigthforward mathematical simplification": on the opposite, the duplication formula will allow you to trade two calls to gamma for one. And instead of squaring abs, you could evaluate Re^2+Im^2 to avoid squaring a square root. –  Mar 07 '14 at 18:16
  • @YvesDaoust Unfortunately, the duplication formula look like `gamma(z) gamma(z+1/2) ~ gamma(2z)`. But in my case, I have to study a term of the form `gamma(z)/gamma(z+1/2)`, so that I always have to deal with two calls to gamma... However, if you find out any faster evaluation of this term, you're more than welcome ! – jibe Mar 10 '14 at 10:32
  • Ooops, yes, I didn't see the '/' :( –  Mar 10 '14 at 10:38

3 Answers3

3

Here are four simple transformations to speed-up your code:

1) Replace pow(x, 2) with x ** 2.0:

def norm_spirals(a,k):
    return pi* abs(gamma(0.25+0.5*k+0.5j*a)/gamma(0.75+0.5*k+0.5j*a)) ** 2.0

2) Factor-out the common subexpression: 0.25+0.5*k+0.5j*a 

def norm_spirals(a,k):
    se = 0.25 + 0.5*k + 0.5j*a
    return pi * abs(gamma(se)/gamma(se + 0.5)) ** 2.0

3) Localize the global lookups:

def norm_spirals(a,k, pi=pi, abs=abs, gamma=gamma):
    se = 0.25 +0.5*k + 0.5j*a
    return pi * abs(gamma(se)/gamma(se + 0.5)) ** 2.0

4) Replace the square operation with b * b. This is called strength reduction:

def norm_spirals(a,k, pi=pi, abs=abs, gamma=gamma):
    se = 0.25 + 0.5*k + 0.5j * a
    b = abs(gamma(se)/gamma(se + 0.5))
    return b * b * pi
Raymond Hettinger
  • 216,523
  • 63
  • 388
  • 485
  • 1
    Is not `value*value` also faster than `value**2` as well? –  Mar 07 '14 at 23:10
  • Yes, that is correct. I've modified the answer accordingly. Thank you. – Raymond Hettinger Mar 07 '14 at 23:23
  • 1
    @RaymondHettinger I used your rewritting and also directly embedded this calculation in the code, so that this procedure is no more called. It sped up the calculation by about 20%, which is definitely non-negligible. – jibe Mar 10 '14 at 10:55
0

I get it 15% faster (python 2.7) by introducing a temp var and remove the call to pow:

def norm_spirals(a, k):
    _ = 0.5 * k + 0.5j * a
    return pi * abs(gamma(0.25 + _) / gamma(0.75 + _)) ** 2

For python 3 the gain was just 11%.

Disclaimer:
I didn't have scipy installed so for gamma I used

def gamma(a):
    return a

So "15% faster" might be a bit misleading.

MyGGaN
  • 1,766
  • 3
  • 30
  • 41
0

If norm_spirals() is invoked many times with the same arguments, then you can try memoization to cache the results and see if it makes a difference.

Community
  • 1
  • 1
Christophe Vu-Brugier
  • 2,645
  • 26
  • 21