2

我想通过以下方式在 python 中执行数学积分：

[1] 在 scipy.optimize.fsolve 的帮助下求解一个隐式方程以找到被积函数的最大位置

[2] 在 scipy.integrate.quad 的帮助下，将被积函数的最大值移到零并从 -10 积分到 10

由于被积函数有一个自由参数 xi，我想在 numpy 的帮助下使用一系列 xi 值来执行此操作，因此我使用了 numpy.vectorize。

现在有两种方法可以向量化这个算法：

[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/sec

单独积分和 fsolve 矢量化：速度 = 2136.75134429 ints/sec

由于这只是我实际问题的简化版本，我需要知道它为什么会这样。有人可以解释一下吗？谢谢！

4

1 回答 1

0

总之，xc当您使用“一次”方法时，这是变量的问题。这是一个ndarray，当math.exp()用xor xi（两个浮点数）调用时，代码会变得更慢。如果您xc=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。Evenmath.exp在“单独”方法中的速度是原来的两倍。

这表明在每种方法中评估的类型是不同的。事实上，这就是重点。当“一次”运行时，您可以打印type(xc)以查看它是一个numpy.ndarray, float64，而在​​“单独”方法中它只是一个float64. 在 a 中对这些类型求和似乎不是一个好主意math.exp()，如下所示：

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和pifromnumpy可以减少差异，但它会使您的代码变慢，可能是因为这些函数也可能返回 a ndarray：

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

ndarray在这种情况下，不要转换为 a 似乎是个好主意。好主意是转换为floatthen ，如下所示（在“一次”方法中是必需的）：

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 回答