有效地计算多项式

Question

我试着用numpy来计算多项式（3'd次）。 我发现用更简单的python代码来做会更有效。在

import numpy as np
import timeit

m = [3,7,1,2]

f = lambda m,x: m[0]*x**3 + m[1]*x**2 + m[2]*x + m[3]
np_poly = np.poly1d(m)
np_polyval = lambda m,x: np.polyval(m,x)
np_pow = lambda m,x: np.power(x,[3,2,1,0]).dot(m)

print 'result={}, timeit={}'.format(f(m,12),timeit.Timer('f(m,12)', 'from __main__   import f,m').timeit(10000))
result=6206, timeit=0.0036780834198

print 'result={}, timeit={}'.format(np_poly(12),timeit.Timer('np_poly(12)', 'from __main__ import np_poly').timeit(10000))
result=6206, timeit=0.180546045303

print 'result={}, timeit={}'.format(np_polyval(m,12),timeit.Timer('np_polyval(m,12)', 'from __main__ import np_polyval,m').timeit(10000))
result=6206, timeit=0.227771043777

print 'result={}, timeit={}'.format(np_pow(m,12),timeit.Timer('np_pow(m,12)', 'from __main__ import np_pow,m').timeit(10000))
result=6206, timeit=0.168987989426

我错过什么了吗？在

在numpy中有没有其他方法来计算多项式？在

Answer 1

大约23年前，我从大学图书馆里查到了一本Press等人用C写的数字食谱。那本书里有很多很酷的东西，但有一段话多年来一直萦绕在我心头，page 173 here：

We assume that you know enough never to evaluate a polynomial this way:

    p=c[0]+c[1]*x+c[2]*x*x+c[3]*x*x*x+c[4]*x*x*x*x;

or (even worse!),

    p=c[0]+c[1]*x+c[2]*pow(x,2.0)+c[3]*pow(x,3.0)+c[4]*pow(x,4.0);

Come the (computer) revolution, all persons found guilty of such criminal behavior will be summarily executed, and their programs won't be! It is a matter of taste, however, whether to write

    p = c[0]+x*(c[1]+x*(c[2]+x*(c[3]+x*c[4])));

or

    p = (((c[4]*x+c[3])*x+c[2])*x+c[1])*x+c[0];

因此，如果您真的担心性能，您想尝试一下，对于高次多项式，差异将是巨大的：

In [24]: fast_f = lambda m, x: m[3] + x*(m[1] + x*(m[2] + x*m[3]))

In [25]: %timeit f(m, 12)
1000000 loops, best of 3: 478 ns per loop

In [26]: %timeit fast_f(m, 12)
1000000 loops, best of 3: 374 ns per loop

如果您想继续使用numpy，在我的系统上有一个更新的多项式类，它比poly1d快2倍，但仍然比前面的循环慢得多：

^{pr2}$

Answer 2

好吧，看看polyval的实现（这是在计算poly1d时最终调用的函数），实现者决定包含一个显式循环似乎很奇怪。。。从numpy 1.6.2的来源：

def polyval(p, x):
    p = NX.asarray(p)
    if isinstance(x, poly1d):
        y = 0
    else:
        x = NX.asarray(x)
        y = NX.zeros_like(x)
    for i in range(len(p)):
        y = x * y + p[i]
    return y

一方面，避免power操作在速度上应该是有利的，另一方面，python级别的循环几乎把事情搞砸了。在

以下是另一种新的实现方式：

^{pr2}$

为了提高速度，我避免在每次调用时重新创建电源阵列。另外，为了公平起见，在对numpy进行基准测试时，应该从numpy数组开始，而不是从列表开始，以避免在每次调用时将list转换为numpy的代价。在

所以，当添加m = np.array(m)时，我上面的g只比你的f慢大约50%。在

尽管在您发布的示例中速度较慢，但是对于计算标量x的低次多项式，您确实不能比显式实现快得多（比如您的f）（当然，您可以可以，但如果不编写较低级别的代码，则可能不会太快）。然而，对于更高的学位（你必须用某种循环来代替你的显式表达式），numpy方法（例如g）会随着度数的增加而证明更快，而且对于向量化的评估，即当x是一个向量时。在