我想通过以下方式在 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 整数/秒
分别对矢量化和 fsolve 进行积分:速度 = 2136.75134429 整数/秒
因为这只是我实际问题的简化版本，所以我需要知道为什么它会这样。有人可以解释一下吗？谢谢!
请您参考如下方法:
总之，当您使用“立即”方法时，这是 xc 变量的问题。它是一个 ndarray，当在 math.exp() 中使用 x 或 xi (均为 float )调用时使代码变得更慢。如果您在“一次”方法中生成 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，它在 integral_with_fsolve 中的运行时间几乎是原来的 7 倍。数值积分 quad 也是如此。即使是 math.exp 在“单独”方法中也快了一倍。
这表明在每种方法中评估的类型是不同的。事实上，这就是重点。当“一次”运行时，您可以打印 type(xc) 以查看它是一个 numpy.ndarray, float64，而在“分别”方法中它只是一个 float64。在 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。使用 numpy 中的 exp、sqrt 和 pi 减少了差异，但它使您的代码慢得多，可能是因为这些函数也可能返回一个 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 似乎是个好主意。好主意是转换为 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
