为什么 Matlab 的积分在 Scipy 中的表现优于Integrated.quad?

计算科学 matlab Python 正交 麻木
2021-11-30 01:31:02

我对 matlab 处理数值积分与 Scipy 的方式感到有些沮丧。我在下面的测试代码中观察到以下差异:

  1. Matlab 版本的运行速度平均比我的 python版本快 24 倍!
  2. Matlab 的版本能够在没有警告的情况下计算积分,而 python 返回nan+nanj

对于提到的两点,我能做些什么来确保我在 python 中获得相同的性能?根据文档,这两种方法都应该使用“全局自适应正交”来近似积分。

下面是两个版本的代码(虽然python要求创建一个积分函数以便它可以处理复杂的被积函数,但非常相似。)

Python

import numpy as np
from scipy import integrate
import time

def integral(integrand, a, b,  arg):
    def real_func(x,arg):
        return np.real(integrand(x,arg))
    def imag_func(x,arg):
        return np.imag(integrand(x,arg))
    real_integral = integrate.quad(real_func, a, b, args=(arg))
    imag_integral = integrate.quad(imag_func, a, b, args=(arg))   
    return real_integral[0] + 1j*imag_integral[0]

vintegral = np.vectorize(integral)


def f_integrand(s, omega):
    sigma = np.pi/(np.pi+2)
    xs = np.exp(-np.pi*s/(2*sigma))
    x1 = -2*sigma/np.pi*(np.log(xs/(1+np.sqrt(1-xs**2)))+np.sqrt(1-xs**2))
    x2 = 1-2*sigma/np.pi*(1-xs)
    zeta = x2+x1*1j
    Vc = 1/(2*sigma)
    theta =  -1*np.arcsin(np.exp(-np.pi/(2.0*sigma)*s))
    t1 = 1/np.sqrt(1+np.tan(theta)**2)
    t2 = -1/np.sqrt(1+1/np.tan(theta)**2)
    return np.real((t1-1j*t2)/np.sqrt(zeta**2-1))*np.exp(1j*omega*s/Vc);

t0 = time.time()
omega = 10
result = integral(f_integrand, 0, np.inf, omega)
print time.time()-t0
print result

MATLAB

function [ out ] = f_integrand( s, omega )
    sigma = pi/(pi+2); 
    xs = exp(-pi.*s./(2*sigma));
    x1 = -2*sigma./pi.*(log(xs./(1+sqrt(1-xs.^2)))+sqrt(1-xs.^2));
    x2 = 1-2*sigma./pi.*(1-xs);
    zeta = x2+x1*1j;
    Vc = 1/(2*sigma);
    theta =  -1*asin(exp(-pi./(2.0.*sigma).*s));
    t1 = 1./sqrt(1+tan(theta).^2);
    t2 = -1./sqrt(1+1./tan(theta).^2);
    out = real((t1-1j.*t2)./sqrt(zeta.^2-1)).*exp(1j.*omega.*s./Vc);
end

t=cputime;
omega = 10;
result = integral(@(s) f_integrand(s,omega),0,Inf)
time_taken = cputime-t
2个回答

这个问题有两个非常不同的子问题。我将只讨论第一个。

Matlab 版本的运行速度平均比我的 python版本快 24 倍!

第二个是主观的。我想说让用户知道积分存在问题是一件好事,而且这种 SciPy 行为优于 Matlab 的行为,可以保持沉默,并以某种方式尝试在内部以只有 Matlab 工程师知道的方式处理它决定它是最好的。

我将集成跨度更改为从030(而不是从0np.inf)以避免 NaN 警告并添加了 JIT 编译。为了对解决方案进行基准测试,我重复了 300 次集成,结果来自我的笔记本电脑。

没有 JIT 编译:

$ ./test_integrate.py
34.20992112159729
(0.2618828053067563+0.24474506983644717j)

使用 JIT 编译:

$ ./test_integrate.py
0.8560323715209961
(0.261882805306756+0.24474506983644712j)

与非 JIT 版本相比,这种方式添加两行代码导致 Python 代码的加速因子约为 40 倍。我的笔记本电脑上没有 Matlab 来提供更好的比较,但是,如果它比 24/40 = 0.6 可以很好地扩展到您的 PC,那么对于这个特定的用户算法,带有 JIT 的 Python 应该几乎是 Matlab 的两倍。完整代码:

#!/usr/bin/env python3
import numpy as np
from scipy import integrate
from numba import complex128,float64,jit
import time

def integral(integrand, a, b,  arg):
    def real_func(x,arg):
        return np.real(integrand(x,arg))
    def imag_func(x,arg):
        return np.imag(integrand(x,arg))
    real_integral = integrate.quad(real_func, a, b, args=(arg))
    imag_integral = integrate.quad(imag_func, a, b, args=(arg))   
    return real_integral[0] + 1j*imag_integral[0]

vintegral = np.vectorize(integral)


@jit(complex128(float64, float64), nopython=True, cache=True)
def f_integrand(s, omega):
    sigma = np.pi/(np.pi+2)
    xs = np.exp(-np.pi*s/(2*sigma))
    x1 = -2*sigma/np.pi*(np.log(xs/(1+np.sqrt(1-xs**2)))+np.sqrt(1-xs**2))
    x2 = 1-2*sigma/np.pi*(1-xs)
    zeta = x2+x1*1j
    Vc = 1/(2*sigma)
    theta =  -1*np.arcsin(np.exp(-np.pi/(2.0*sigma)*s))
    t1 = 1/np.sqrt(1+np.tan(theta)**2)
    t2 = -1/np.sqrt(1+1/np.tan(theta)**2)
    return np.real((t1-1j*t2)/np.sqrt(zeta**2-1))*np.exp(1j*omega*s/Vc);

t0 = time.time()
omega = 10
for i in range(300): 
    #result = integral(f_integrand, 0, np.inf, omega)
    result = integral(f_integrand, 0, 30, omega)
print (time.time()-t0)
print (result)

注释掉 @jit 行以查看您的 PC 的差异。

有时,要集成的功能无法 JIT。在这种情况下,使用另一种集成方法将是解决方案。

我会推荐scipy.integrate.romberg (ref)romberg可以集成复杂的函数,可以用数组对函数求值。