13
$\begingroup$

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
$\endgroup$
  • 4
    $\begingroup$ You should be happy that Python is only 25x slower (and not 250x). $\endgroup$ – stali Oct 19 '14 at 2:10
  • 4
    $\begingroup$ Because you are calling a python-function in a loop over and over again (hid by np.vectorize). Try doing calculations on the entire array at once. It that is not possible, have a look at numba or also Cython, but I hope the latter is not necessary. $\endgroup$ – sebix Oct 19 '14 at 15:53
  • 2
    $\begingroup$ "global adaptive quadrature" indicates that it adapts until it reaches a certain precision. To make sure that you're comparing the same thing look for the parameter (there surely is one) that sets the precision and set it for both. $\endgroup$ – bgschaid Oct 19 '14 at 21:49
  • 2
    $\begingroup$ Regarding @bgschaid's comment, integral's default absolute and relative tolerances are 1e-10 and 1e-6, respectively. integrate.quad specifies these both as 1.49e-8. I don't see where integrate.quad is described as a "global adaptive" method and it is most certainly different from the (adaptive Gauss-Kronrod, I believe) method used by integral. I'm not sure what the "global" part means, myself. Also, it's never a good idea to be using cputime instead of tic/toc or time it. $\endgroup$ – horchler Oct 20 '14 at 3:13
  • 5
    $\begingroup$ Before anything I'd check whether the problem is the algorithm or the language: add a global counter variable that is incremented inside the functions. After the integration this should tell you how often each function is evaluated. If these counters differ significantly then at least part of the problem is that MATLAB uses the better algorithm $\endgroup$ – bgschaid Oct 21 '14 at 8:34
15
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.