Why does Matlab's integral outperform integrate.quad in Scipy?

Question

I am experiencing some frustration over the way matlab handles numerical integration vs. Scipy. I observe the following differences in my test code below:

1. Matlab's version runs on average 24 times faster than my python equivalent!
2. Matlab's version is able to calculate the integral without warnings, while python returns nan+nanj

What can I do to ensure I get the same performance in python with respect to the two points mentioned? According to documentation both methods should be using a "global adaptive quadrature" to approximate the integral.

Below is the code in the two versions (fairly similar although python requires that an integral function is created so that it can handle complex integrands.)

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

Answer 1

The question has two very different subquestions. I will address the first one only.

Matlab's version runs on average 24 times faster than my python equivalent!

The second one is subjective. I would say that letting know the user that there is some problem with the integral is a good thing and this SciPy behavior outperforms the Matlab`s one to keep it silent and somehow try to deal with it internally the way only known by Matlab engineers who decided it to be the best.

I changed integration span to be from 0 to 30 (instead of from 0 to np.inf) to avoid NaN waring and added a JIT compilation. To benchmark the solution I repeated the integration for 300 times, results are from my laptop.

Without JIT compilation:

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

With JIT compilation:

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

This way adding two lines of code leads to the speedup factor of about 40 times of Python code compared to a non-JIT version. I have no Matlab on my laptop to provide a better comparison, however, if it scales well to your PC than 24/40 = 0.6, so Python with JIT should be almost twice as fast as Matlab for this particular user algorithm. Full code:

#!/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)

Comment out the @jit line to see the difference for your PC.

Answer 2

Sometimes the function to integrate cannot be JITed. In that case, using another integration method would be the solution.

I would recommend scipy.integrate.romberg (ref). romberg can integrate complex functions, and can evaluate the function with array.