fastest polynomial evaluation in python

Question

I'm working on an old project of building Newton Basins and I'm trying to make it as fast as possible. The first thing I'm trying to speed up is how to evaluate a polynomial function at a given complex point x0. I thought of 4 different ways of doing this and tested them with timeit. The code I used is the following:

import timeit
import numpy as np
import random
from numpy import polyval

class Test(object):
    re = random.randint(-40000, 40000)/10000
    im = random.randint(-40000, 40000)/10000
    x0 = complex(re, im)
    coefs = np.array([48,8,4,-10,2,-3,1])
    flip_coefs = np.flip(coefs)

    def solve0():
        y = np.array([Test.coefs[i]*(Test.x0**i) for i in range(len(Test.coefs))]).sum()
        return y

    def solve1():
        y = 0
        for i in range(len(Test.coefs)):
            y += Test.coefs[i]*(Test.x0**i)
        return y

    def solve2():
        y = np.dot(Test.coefs,Test.x0**np.arange(len(Test.coefs)))
        return y

    def solve3():
        y = polyval(Test.flip_coefs, Test.x0)
        return y


Test.solve0()


if __name__ == '__main__':
    print(timeit.timeit('Test.solve0()', setup="from __main__ import Test", number=10000))
    print(timeit.timeit('Test.solve1()', setup="from __main__ import Test", number=10000))
    print(timeit.timeit('Test.solve2()', setup="from __main__ import Test", number=10000))
    print(timeit.timeit('Test.solve3()', setup="from __main__ import Test", number=10000))

the thing is that I was pretty sure that numpy.polyval() would be the fastest, but it seems that, on Linux, np.dot(coefs,x**np.arange(len(coefs))) is more than twice as fast, regardless of the value of x0 (I don't know if it is the same in Windows and MacOS). This is an output example I've got:

0.1735437790002834
0.12607222800033924
0.0313361469998199
0.0796813930001008

This seems quite strange since numpy.polyval() was specifically built for solving polynomials. So, my questions are: Is there something I'm missing here (maybe related to the coefficients I chose)? Are there faster ways of evaluating polynomials?

user_14604598 · Answer

Examining the source code of the polyval() function of numpy you'll observe that this is a purely pythonic function. Numpy uses Horner's method for polynomial evaluation (and facilitates the evaluation of multiple points concurrently, though this case doesn't apply here).

To answer your questions about evaluation cost, I believe Horner's method is optimal from the algorithmic perspective: It allows the evaluation of a polynomial of degree n with only n multiplications and n additions, in contrast to your fastest method solve2() where multiplication operations scale polynomially with input size. In terms of speed of computation, I suspect that the increase in performance of your approach vs the polyval() function comes from the fact that the dot product in numpy is implemented in C and receives a substantial speedup on that front.

Addendum: For a comprehensive overview of the topic, I cannot recommend enough Knuth's The Art of Computer Programming. Vol. 2: Seminumerical Algorithms (3rd ed.), of which section 4.6.4 on page 485 refers to this situation. I believe Horner's method is optimal for a single polynomial evaluation, however if this polynomial is to be evaluated many times, then potentially further speedups can be gained if one allows preconditioning.

fastest polynomial evaluation in python

1 Answers