Dlaczego całka Matlaba osiąga lepsze wyniki niż integracja.quad w Scipy?

13

Mam frustrację związaną ze sposobem, w jaki Matlab obsługuje integrację numeryczną vs. Scipy. Poniżej obserwuję następujące różnice w moim kodzie testowym:

  1. Wersja Matlaba działa średnio 24 razy szybciej niż mój odpowiednik w Pythonie!
  2. Wersja Matlaba jest w stanie obliczyć całkę bez ostrzeżeń, podczas gdy Python powraca nan+nanj

Co mogę zrobić, aby uzyskać taką samą wydajność w Pythonie w odniesieniu do dwóch wymienionych punktów? Zgodnie z dokumentacją obie metody powinny wykorzystywać „globalną kwadraturę adaptacyjną” w celu przybliżenia całki.

Poniżej znajduje się kod w dwóch wersjach (dość podobny, chociaż python wymaga utworzenia funkcji całkowej, aby mogła obsługiwać złożone całki).

Pyton

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
Dipol
źródło
4
Powinieneś być szczęśliwy, że Python jest tylko 25x wolniejszy (a nie 250x).
stali
4
Ponieważ wywołujesz funkcję python w pętli w kółko (ukryty przez np.vectorize). Spróbuj wykonać obliczenia dla całej tablicy jednocześnie. Nie jest to możliwe, spójrz na Numba lub Cython, ale mam nadzieję, że to drugie nie jest konieczne.
sebix,
2
„globalna kwadratura adaptacyjna” wskazuje, że dostosowuje się do osiągnięcia określonej precyzji. Aby upewnić się, że porównujesz to samo, poszukaj parametru (na pewno jest taki), który ustawia precyzję i ustaw ją dla obu.
bgschaid
2
Odnośnie @ komentarzu bgschaid, w integral„s domyślne bezwzględne i względne tolerancje 1e-10i 1e-6odpowiednio. integrate.quadokreśla oba jako 1.49e-8. Nie widzę, gdzie integrate.quadjest opisana jako metoda „globalnej adaptacji” i na pewno różni się ona od metody (adaptacyjnej metody Gaussa-Kronroda, jak sądzę) integral. Sam nie jestem pewien, co oznacza ta „globalna” część. Ponadto nigdy nie warto używać cputimezamiast tic/toc lub time it.
horchler
5
Przede wszystkim sprawdzę, czy problemem jest algorytm czy język: dodaj globalną zmienną licznika, która jest zwiększana wewnątrz funkcji. Po integracji powinno to powiedzieć, jak często każda funkcja jest oceniana. Jeśli te liczniki różnią się znacznie, to przynajmniej część problemu polega na tym, że MATLAB używa lepszego algorytmu
bgschaid

Odpowiedzi:

15

Pytanie ma dwa bardzo różne pytania. Zajmę się tylko pierwszym.

Wersja Matlaba działa średnio 24 razy szybciej niż mój odpowiednik w Pythonie!

Drugi jest subiektywny. Powiedziałbym, że poinformowanie użytkownika, że ​​jest jakiś problem z całką, jest dobrą rzeczą i to zachowanie SciPy przewyższa Matlaba, jeśli chodzi o milczenie i jakoś spróbować poradzić sobie z tym wewnętrznie w sposób znany tylko inżynierom Matlaba, którzy zdecydowałem, że będzie najlepszy.

Zmieniłem zakres integracji na z 0 na 30 (zamiast z 0 na np.inf ), aby uniknąć ostrzeżenia NaN i dodałem kompilację JIT. Aby przetestować rozwiązanie, powtórzyłem integrację 300 razy, wyniki pochodzą z mojego laptopa.

Bez kompilacji JIT:

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

Z kompilacją JIT:

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

W ten sposób dodanie dwóch wierszy kodu prowadzi do przyspieszenia około 40 razy kodu Pythona w porównaniu z wersją inną niż JIT. Nie mam Matlaba na moim laptopie, aby zapewnić lepsze porównanie, jednak jeśli skaluje się dobrze na twoim komputerze niż 24/40 = 0,6, więc Python z JIT powinien być prawie dwa razy szybszy niż Matlab dla tego konkretnego algorytmu użytkownika. Pełny kod:

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

Skomentuj linię @jit, aby zobaczyć różnicę na swoim komputerze.

kostyfisik
źródło
1

Czasami funkcja integracji nie może być JIT. W takim przypadku rozwiązaniem może być zastosowanie innej metody integracji.

Polecam scipy.integrate.romberg (ref) . rombergpotrafi integrować złożone funkcje i oceniać funkcję z tablicą.

Calsina Yo
źródło