6

numpy を使用して多項式 (3 次元) を評価しようとしています。より単純な 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 レベルのループはかなりのことを台無しにします。

これは、代替の派手な実装です。

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(もちろんできますが、おそらくに頼らなければそれほどではありません)低レベルのコードを書く)。ただし、次数が高い場合 (明示的な式をある種のループに置き換える必要がある場合)、numpy アプローチ (たとえばg、x

于 2014-06-05T17:44:50.263 に答える