¿Por qué el rendimiento superior integral de Matlab integra.quad en Scipy?

13

Estoy experimentando cierta frustración por la forma en que matlab maneja la integración numérica frente a Scipy. Observo las siguientes diferencias en mi código de prueba a continuación:

  1. ¡La versión de Matlab funciona en promedio 24 veces más rápido que mi equivalente en Python!
  2. La versión de Matlab es capaz de calcular la integral sin advertencias, mientras que Python regresa nan+nanj

¿Qué puedo hacer para asegurarme de obtener el mismo rendimiento en Python con respecto a los dos puntos mencionados? Según la documentación, ambos métodos deberían usar una "cuadratura adaptativa global" para aproximar la integral.

A continuación se muestra el código en las dos versiones (bastante similar, aunque Python requiere que se cree una función integral para que pueda manejar integrandos complejos).

Pitón

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
Dipolo
fuente
44
Debería estar contento de que Python sea solo 25 veces más lento (y no 250 veces).
stali
44
Porque estás llamando a una función python en un ciclo una y otra vez (oculto por np.vectorize). Intente hacer cálculos en toda la matriz a la vez. Si eso no es posible, eche un vistazo a Numba o también a Cython, pero espero que esto último no sea necesario.
sebix
2
La "cuadratura adaptativa global" indica que se adapta hasta que alcanza una cierta precisión. Para asegurarse de que está comparando lo mismo, busque el parámetro (seguramente hay uno) que establece la precisión y lo establece para ambos.
bgschaid
2
Con respecto al comentario de @ bgschaid, integrallas tolerancias absolutas y relativas predeterminadas son 1e-10y 1e-6, respectivamente. integrate.quadespecifica estos dos como 1.49e-8. No veo dónde integrate.quadse describe como un método "global adaptativo" y es ciertamente diferente del método (adaptable Gauss-Kronrod, creo) utilizado por integral. No estoy seguro de lo que significa la parte "global", yo mismo. Además, nunca es una buena idea usarlo en cputimelugar de tic/ toco time it.
horchler
55
Antes de nada, comprobaría si el problema es el algoritmo o el idioma: agregue una variable de contador global que se incremente dentro de las funciones. Después de la integración, esto debería decirle con qué frecuencia se evalúa cada función. Si estos contadores difieren significativamente, entonces al menos parte del problema es que MATLAB usa el mejor algoritmo
bgschaid

Respuestas:

15

La pregunta tiene dos preguntas muy diferentes. Me referiré solo al primero.

¡La versión de Matlab funciona en promedio 24 veces más rápido que mi equivalente en Python!

El segundo es subjetivo. Diría que dejarle saber al usuario que hay algún problema con la integral es algo bueno y que este comportamiento de SciPy supera al de Matlab para mantenerlo en silencio y de alguna manera tratar de tratarlo internamente de la manera que solo conocen los ingenieros de Matlab que Decidí que era el mejor.

Cambié el intervalo de integración para que sea de 0 a 30 (en lugar de 0 a np.inf ) para evitar la advertencia de NaN y agregué una compilación JIT. Para comparar la solución, repetí la integración 300 veces, los resultados son de mi computadora portátil.

Sin compilación JIT:

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

Con la compilación JIT:

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

De esta forma, agregar dos líneas de código conduce al factor de aceleración de aproximadamente 40 veces el código de Python en comparación con una versión que no sea JIT. Sin embargo, no tengo Matlab en mi computadora portátil para proporcionar una mejor comparación, si se adapta bien a su PC que 24/40 = 0.6, por lo que Python con JIT debería ser casi el doble de rápido que Matlab para este algoritmo de usuario en particular. Código completo:

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

Comente la línea @jit para ver la diferencia para su PC.

kostyfisik
fuente
1

Algunas veces la función a integrar no se puede JITAR. En ese caso, usar otro método de integración sería la solución.

Yo recomendaría scipy.integrate.romberg (ref) . rombergpuede integrar funciones complejas y puede evaluar la función con una matriz.

calsina yo
fuente