6

我正在尝试使用 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 中是否有另一种方法来评估多项式？

4

2 回答 2

8

大约在 23 年前，我从大学图书馆查阅了 Press et al Numerical Recipes in C的副本。那本书里有很多很酷的东西，但有一段话多年来一直困扰着我，第 173 页：

我们假设您知道的足够多，永远不会以这种方式评估多项式：

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

或者（更糟！），

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

（计算机）革命来了，所有犯有此类犯罪行为的人都将被即决处决，他们的程序不会被处决！然而，是否写是一个品味问题

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

或者

    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
于 2014-06-05T18:16:52.130 回答
1

好吧，看看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

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

这是另一种 numpy-ish 实现：

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上面的运行速度只比你的慢 50% f。

尽管在您发布的示例中速度较慢，但​​对于在标量上评估低次多项式x，您确实不能比显式实现（如您的）快得多f（当然您可以，但如果不求助于编写低级代码）。但是，对于更高的度数（您必须用某种循环替换显式表达式），g随着度数的增加，numpy 方法（例如）将证明更快，并且对于向量化评估，即何时x是向量。

于 2014-06-05T17:44:50.263 回答