2

次の方法でPythonで数学的積分を実行したい：

[1] scipy.optimize.fsolve を使用して暗黙の方程式を解き、被積分関数の最大位置を見つけます。

[2] 被積分関数の最大値をゼロにシフトし、scipy.integrate.quad を使用して -10 から 10 まで積分します。

被積分関数には 1 つの自由パラメーター xi があるため、numpy を使用して xi 値の範囲でこれを実行したいので、numpy.vectorize を使用します。

このアルゴリズムをベクトル化するには、次の 2 つの方法があります。

[A] [1] と [2] を別々にベクトル化し、vec_[1] の結果を vec_[2] に入力として与える

[B] [1] と [2] を同時に実行する関数をベクトル化する

[A] は [B] よりもはるかに高速であることに気付きました。コードは次のとおりです。

from scipy import optimize, integrate
import numpy as np
from math import exp, sqrt, pi
import time

def integral_with_fsolve(xi):
    xc = optimize.fsolve(lambda x: x+1./(1+exp(-x))-xi,0.)
    def integrand(x,xi):
        return exp(-(x-xi+xc)**2)/(2.*sqrt(2.*pi))/(1.+exp(x+xc))
    integral = integrate.quad(integrand,-10.,10.,args=(xi,),epsabs=0.)
    return integral[0]

def integral(xi,xc):
    def integrand(x,xi):
        return exp(-(x-xi+xc)**2)/(2.*sqrt(2.*pi))/(1.+exp(x+xc))
    integral = integrate.quad(integrand,-10.,10.,args=(xi,),epsabs=0.)
    return integral[0]

def fsolve(xi):
    return optimize.fsolve(lambda x: x+1./(1+exp(-x))-xi,0.)

vec_integral_with_fsolve = np.vectorize(integral_with_fsolve)
vec_integral = np.vectorize(integral)
vec_fsolve = np.vectorize(fsolve)

xi = np.linspace(0.,2.,1000)

t0=time.time()
something = vec_integral_with_fsolve(xi)
dur=(time.time()-t0)
speed = xi.size/dur
print('Integrate and fsolve vectorized in one: speed = {} ints/sec'.format(speed))

t0=time.time()
xc = vec_fsolve(xi)
something = vec_integral(xi,xc)
dur=(time.time()-t0)
speed = xi.size/dur
print('Integrate and fsolve vectorized seperately: speed = {} ints/sec'.format(speed))

出力は常に次のようになります

ベクトル化された積分と fsolve: 速度 = 298.151473998 ints/秒

個別にベクトル化された積分と fsolve: 速度 = 2136.75134429 ints/秒

これは実際の問題の単純化されたバージョンにすぎないため、なぜこのように動作するのかを知る必要があります。誰か説明してくれませんか？ありがとう！

4

1 に答える 1

0

要約すると、これはxc「一度に」アプローチを使用する場合の変数の問題です。or (両方のフロート)ndarrayで呼び出されると、math.exp()コードが大幅に遅くなります。「一度に」アプローチを行うと、「別々に」アプローチと実質的に同じパフォーマンスが得られます。xxixc=float(xc)

その下に、それを見つける方法に関する詳細な説明が続きます。

これを使用cProfileすると、ボトルネックがどこにあるかを簡単に確認できます。

AT ONCE:
ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 1001    0.002    0.000    2.040    0.002 quadpack.py:135(quad)
 1001    0.001    0.000    2.038    0.002 quadpack.py:295(_quad)
 1001    0.002    0.000    2.143    0.002 tmp.py:15(integral_with_fsolve)
231231   1.776    0.000    1.925    0.000 tmp.py:17(integrand)
470780   0.118    0.000    0.118    0.000 {math.exp}
 1001    0.112    0.000    2.037    0.002 {scipy.integrate._quadpack._qagse}

SEPARATELY:
ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 1001    0.001    0.000    0.340    0.000 quadpack.py:135(quad)
 1001    0.001    0.000    0.339    0.000 quadpack.py:295(_quad)
 1001    0.001    0.000    0.341    0.000 tmp.py:9(integral)
231231   0.200    0.000    0.278    0.000 tmp.py:10(integrand)
470780   0.054    0.000    0.054    0.000 {math.exp}
 1001    0.060    0.000    0.338    0.000 {scipy.integrate._quadpack._qagse}

全体のタイミングは次のとおりです。

AT ONCE:
1 loops, best of 3: 1.91 s per loop
SEPARATELY:
1 loops, best of 3: 312 ms per loop

最大の違いは でintegrand、実行に約 7 倍の時間がかかりますintegral_with_fsolve。同じことが数値積分にも起こりquadます。math.exp「個別」アプローチでは、2 倍の速さでもあります。

各アプローチで評価される型が異なることを示唆しています。実際、それがポイントです。「一度に」実行するtype(xc)と、印刷しnumpy.ndarray, float64てfloat64. math.exp()以下に示すように、これらの型を a 内で合計することはお勧めできません。

xa = -0.389760856858
xc = np.array([[-0.389760856858]],dtype='float64')

timeit for i in range(1000000): exp(xc+xa)
#1 loops, best of 3: 1.96 s per loop

timeit for i in range(1000000): exp(xa+xa)
#10 loops, best of 3: 173 ms per loop

どちらの場合もmath.exp()を返しますfloat。exp、sqrtおよびpifromを使用numpyすると差は縮まりますが、おそらくこれらの関数も を返す可能性があるため、コードが大幅に遅くなりますndarray。

AT ONCE:
1 loops, best of 3: 4.46 s per loop
SEPARATELY:
1 loops, best of 3: 2.14 s per loop

この場合、 a に変換しないことをお勧めしndarrayます。以下に示すように、 thenに変換することをお勧めしfloatます (「一度に」アプローチで必要なだけです)。

def integral_with_fsolve(xi):
    xc = optimize.fsolve(lambda x: x+1./(1+exp(-x))-xi,0.)
    xc = float(xc) # <-- SEE HERE
    def integrand(x,xi):
        return exp(-(x-xi+xc)**2)/(2.*sqrt(2.*pi))/(1.+exp(x+xc))
    integral = integrate.quad(integrand,-10.,10.,args=(xi,),epsabs=0.)
    return integral[0]

新しいタイミングで:

AT ONCE:
1 loops, best of 3: 321 ms per loop
SEPARATELY:
1 loops, best of 3: 315 ms per loop
于 2013-05-14T19:29:20.633 に答える