为什么Python的numpy.polyval()如此缓慢?

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.

我正致力于可视化Julia & Mandelbrot集合和牛顿分形——为此，我需要在复平面中计算许多值。我可以使用任何类型的数学函数，但是对于多项式来说已经足够了。

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.

我需要计算函数/多项式的导数和值，所以我查看了numpy模块，找到了numpy.polyder()和numpy.polyval()。这似乎正是我所需要的，但我的脚本突然变得非常慢。

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):

基本上这个脚本计算了很多多项式g(C) = C* 2的值。时间结果为(程序的实际输出):

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?

有什么方法可以让它更快?我知道调用另一个函数是很耗时的，但仍然如此。我是否应该重新考虑仅仅在一个地方改变多项式系数的优点而不是在时间上改变它的缺点?关于如何处理这类问题，还有其他建议吗?

2 个解决方案

#1

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.

使用Horner方法时，乘法的顺序出现了一个奇怪的现象，其速度大约是原来的5倍。它被固定在numpy 1。10中，或者你可以使用numpi .多项式。多项式。polyval函数，已经被固定了。可以进一步向量化，例如使用numpi .多项式中的一个二维系数数组。polyval版本，但我不太清楚您确切的要求是什么。还需要注意的是，可以直接将多项式迭代为p(p)，但在某种程度上，我认为这会导致数值问题。

#2

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:

这本身并不是一个答案，但是我编写了一个函数polyval_factory来生成一个mypolyval函数，它给出了一个系数数组，比如pol。它以字符串形式生成快速表达式，如C*C + 1.0 *C*C *C*C + 2.0，然后将它们装入lambda函数。基本上，它使用字符串来代替完整的符号代数库。这个函数在下面的示例中定义并测试，其工作速度几乎与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:

编辑:polyval_factory:现在运行mypolyval = polyval_factory([2.0,3.0,1.0])做了适当的事情并打印:

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