Why is Python's numpy.polyval() so slow?

Question

Update There was a mistake in the script.

I am working on visualizing Julia & Mandelbrot set and also Newton fractals - for this I need calculating a lot of values in the complex plane. I can use any type of mathematical functions I want, but it is enough for polynomials.

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

I tried to think of some easy test to show the difference in time. For that purpose 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. The results in time are (the actual output of the program):

Polyval: 2.34903216362
0j
Standard: 0.0198249816895
0j

I might have not set this test in the best way, doing something like this for the first time. But even if there would be any mistake, running my other scripts show great difference in the time.

Question

Is there a way how to make it faster? I understand calling another function is time consuming, but still. Should I rethink the advantage of being able to change the polynomial coefficients only at one place against the disadvantage in time? Any other suggestions on how to deal with such issue?

Answer 1

Not an answer per se, but I wrote a function polyval_factory to generate a mypolyval function given an array of coefficients like pol. It generates fast expressions like C*C + 1.0 * C *C*C + 2.0, in string form, and then bottles them in a lambda function. Basically it uses strings in lieu of a full symbolic algebra library like sympy. This function 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

The Output is :

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 running mypolyval = polyval_factory([2.0,3.0,1.0]) does the appropriate thing and prints:

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

Answer 2

An oddity in the order of multiplication when using the Horner method accounts for 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. It may be possible to vectorize further, for instance by using a 2d array of coefficients in the numpy.polynomial.polynomial.polyval version, but I am not clear on exactly your requirements are. Also note that you can iterate polynomials directly as p(p), but at some point I expect that will cause numerical problems.