Pourquoi la surperformance intégrale de Matlab integr.quad dans Scipy?

13

Je ressens une certaine frustration quant à la façon dont matlab gère l'intégration numérique par rapport à Scipy. J'observe les différences suivantes dans mon code de test ci-dessous:

  1. La version de Matlab tourne en moyenne 24 fois plus vite que mon équivalent python!
  2. La version de Matlab est capable de calculer l'intégrale sans avertissements, tandis que python retourne nan+nanj

Que puis-je faire pour m'assurer d'obtenir les mêmes performances en python par rapport aux deux points mentionnés? Selon la documentation, les deux méthodes devraient utiliser une "quadrature adaptative globale" pour approximer l'intégrale.

Vous trouverez ci-dessous le code dans les deux versions (assez similaire, bien que python nécessite la création d'une fonction intégrale pour pouvoir gérer des intégrales complexes.)

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
Dipôle
la source
4
Vous devriez être heureux que Python soit seulement 25x plus lent (et non 250x).
stali
4
Parce que vous appelez une fonction python en boucle encore et encore (caché par np.vectorize). Essayez de faire des calculs sur l'ensemble du tableau à la fois. Si ce n'est pas possible, jetez un œil à numba ou aussi à Cython, mais j'espère que ce dernier n'est pas nécessaire.
sebix
2
La "quadrature adaptative globale" indique qu'elle s'adapte jusqu'à atteindre une certaine précision. Pour vous assurer que vous comparez la même chose, recherchez le paramètre (il y en a sûrement un) qui définit la précision et définissez-le pour les deux.
bgschaid
2
En ce qui concerne le commentaire de @ bgschaid, integralles tolérances absolue et relative par défaut sont 1e-10et 1e-6, respectivement. integrate.quadspécifie ces deux comme 1.49e-8. Je ne vois pas où integrate.quadest décrite comme une méthode "globale adaptative" et elle est très certainement différente de la méthode (adaptative Gauss-Kronrod, je crois) utilisée par integral. Je ne sais pas trop ce que signifie la partie "globale". De plus, ce n'est jamais une bonne idée d'utiliser à la cputimeplace de tic/ tocou time it.
horchler
5
Avant tout, je vérifierais si le problème est l'algorithme ou le langage: ajoutez une variable de compteur globale qui est incrémentée à l'intérieur des fonctions. Après l'intégration, cela devrait vous dire à quelle fréquence chaque fonction est évaluée. Si ces compteurs diffèrent de manière significative, alors au moins une partie du problème est que MATLAB utilise le meilleur algorithme
bgschaid

Réponses:

15

La question comporte deux sous-questions très différentes. Je ne traiterai que du premier.

La version de Matlab tourne en moyenne 24 fois plus vite que mon équivalent python!

Le second est subjectif. Je dirais que faire savoir à l'utilisateur qu'il y a un problème avec l'intégrale est une bonne chose et ce comportement SciPy surpasse celui de Matlab pour le garder silencieux et essayer de le traiter en interne de la manière connue uniquement par les ingénieurs de Matlab qui a décidé que c'était le meilleur.

J'ai changé la durée d'intégration pour passer de 0 à 30 (au lieu de 0 à np.inf ) pour éviter les avertissements NaN et j'ai ajouté une compilation JIT. Pour comparer la solution, j'ai répété l'intégration 300 fois, les résultats proviennent de mon ordinateur portable.

Sans compilation JIT:

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

Avec la compilation JIT:

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

De cette façon, l'ajout de deux lignes de code entraîne un facteur d'accélération d'environ 40 fois le code Python par rapport à une version non JIT. Je n'ai pas de Matlab sur mon ordinateur portable pour fournir une meilleure comparaison, cependant, s'il évolue bien sur votre PC que 24/40 = 0,6, donc Python avec JIT devrait être presque deux fois plus rapide que Matlab pour cet algorithme utilisateur particulier. Code complet:

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

Mettez en commentaire la ligne @jit pour voir la différence pour votre PC.

kostyfisik
la source
1

Parfois, la fonction à intégrer ne peut pas être JIT. Dans ce cas, l'utilisation d'une autre méthode d'intégration serait la solution.

Je recommanderais scipy.integrate.romberg (ref) . rombergpeut intégrer des fonctions complexes et peut évaluer la fonction avec un tableau.

calsina yo
la source