Why is Python numpy.polyval () so slow?

Question

Update An error has occurred in the script.

I am working on visualizing a set of Julia and Mandelbrot, as well as Newton fractals - for this I need to calculate many values ​​in the complex plane. I can use any type of mathematical function that I want, but that’s enough for polynomials.

I need to calculate the derivative and function / polynomial value, so I looked into the numpy module and found out about numpy.polyder() and numpy.polyval() . It was like what I needed, but suddenly my scripts became very slow.

I tried to come up with some simple test to show the time difference. For this, I wrote the following script:

 import numpy as np import cmath import time from itertools import product C = 0.37 + 0.45j pol = [1,0,0] start_time = time.time() for i in xrange(100000): C = np.polyval(pol, C) print "Polyval: {}".format( time.time() - start_time ) print C C = 0.37 + 0.45j # forgot to reassign the initial value of C start_time = time.time() for i in xrange(100000): C = C**2 print "Standard: {}".format( time.time() - start_time ) print C 

Basically this script calculates a lot of values ​​for the polynomial g (C) = C ** 2. Results in time (actual program output):

 Polyval: 2.34903216362 0j Standard: 0.0198249816895 0j 

I would not put this test in the best way, having done something similar for the first time. But even if there was some kind of error, running my other scripts shows a big time difference.

Question

Is there any way to do this faster? I understand that calling another function takes a lot of time, but still. Should I rethink the advantage of being able to change polynomial coefficients in only one place versus lack of time? Any other suggestions on how to deal with such a problem?

Answer 1

Oddness in multiplication order using the Horner method takes into account a factor of about 5x in speed. This is fixed in numpy 1.10, or you can use numpy.polynomial.polynomial.polyval , which is already fixed. Perhaps you can even further vectorize, for example, using a 2d array of coefficients in the version of numpy.polynomial.polynomial.polyval , but I do not quite understand what exactly your requirements are. Also note that you can iterate the polynomials directly as p(p) , but at some point I expect this to cause numerical problems.

Answer 2

Not an answer in itself, but I wrote a function polyval_factory to generate the function mypolyval , given an array of coefficients of type pol . It generates quick expressions, such as C*C + 1.0 * C *C*C + 2.0 , in string form, and then their bottles in lambda functions. It mainly uses strings instead of a full library of symbolic algebra, such as sympy . This function is defined and tested in the example below and works almost as fast as C*C :

 import numpy as np import cmath import time from itertools import product def polyval_factory(pol): """ Generate a lambda function for evaluating a given polynomial pol : a list of coefficients with highest degree first Note: this function basically uses strings in lieu of a fully symbolic algebra package like sympy """ poly_string_list = [] for i, coeff in enumerate(pol[::-1]): if np.abs(coeff) > 1e-10: if i > 1: poly_string_list.append( repr(coeff) + '*' + '*'.join(['x']*i)) elif i == 1: poly_string_list.append(repr(coeff)+'*x') elif i ==0: poly_string_list.append(repr(coeff)) lambda_str = 'lambda x :' + '+'.join(poly_string_list) print "The prepared lambda function is: \""+lambda_str + "\"" return eval(lambda_str) C = 0.37 + 0.45j pol = [1,0,0] numiter = 30000 start_time = time.time() for i in xrange(numiter): C = np.polyval(pol, C) print "Polyval: {}".format( time.time() - start_time ) print C C = 0.37 + 0.45j # forgot to reassign the initial value of C print "" print "Generating lambda function..." mypolyval = polyval_factory(pol) # generate the lambda function print "" start_time = time.time() for i in xrange(numiter): C = mypolyval(C) print "Polyval_factory: {}".format( time.time() - start_time ) print C C = 0.37 + 0.45j # forgot to reassign the initial value of C print "" start_time = time.time() for i in xrange(numiter): C = C**2 print "Standard: {}".format( time.time() - start_time ) print C 

Output:

 Polyval: 0.738290071487 0j Generating lambda function... The prepared lambda function is: "lambda x :1*x*x" Polyval_factory: 0.013610124588 0j Standard: 0.00678110122681 0j 

Edit : polyval_factory : now working mypolyval = polyval_factory([2.0,3.0,1.0]) does the corresponding thing and prints:

  The prepared lambda function is: "lambda x :1.0+3.0*x+2.0*x*x"