NumPy高效计算多项式

Question

NumPy高效计算多项式

6

我正在尝试使用numpy来评估三次多项式，但发现使用更简单的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中还有其他方法可以求解多项式吗？

- mclafee

2个回答

1

好的，看一下polyval的实现（当你eval一个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

一方面，避免使用幂运算应该有助于提高速度，另一方面，Python级别的循环几乎会使事情变得混乱。

以下是一种替代的numpy实现：

POW = np.arange(100)[::-1]
def g(m, x):
    return np.dot(m, x ** POW[m.size : ])

为了提高速度，我避免在每次调用时重新创建幂数组。此外，在与numpy进行基准测试时，应该从numpy数组开始，而不是列表，以避免在每次调用时将列表转换为numpy的惩罚。m = np.array(m)添加后，我的g大约只比你的f慢50%。

尽管在您发布的示例上速度较慢，但对于在标量x上评估低阶多项式，您确实无法比显式实现（如您的f）更快（当然您可以，但可能不会太多，而不必诉诸编写更低级别的代码）。但是，对于更高的阶数（需要用某种循环替换您的显式表达式），随着阶数的增加，numpy方法（例如g）将证明更快，并且对于矢量化评估，即x为向量时也是如此。

- shx2

“实现者决定包含一个显式循环似乎有些奇怪”，如果您尝试在一个块中完成所有操作，可能会导致内存错误。 - endolith

2

numpy源代码是在霍纳形式下进行评估，这可能不容易矢量化，但要比您编写的方法优越得多。（请参见此问题的其他答案。） - Mike

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Jaime · Accepted Answer

大约23年前，我从大学图书馆借阅了Press等人所著的 Numerical Recipes in C。这本书里有很多很酷的内容，但其中一段话至今仍深深印在我的脑海中，详情请见此处第173页：

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倍，但仍然比以前的循环慢得多：

In [27]: np_fast_poly = np.polynomial.polynomial.Polynomial(m[::-1])

In [28]: %timeit np_poly(12)
100000 loops, best of 3: 15.4 us per loop

In [29]: %timeit np_fast_poly(12)
100000 loops, best of 3: 8.01 us per loop