Algorithme rapide pour la fonction log gamma

Le temps d'exécution de votre fonction sera mis à l'échelle linéairement (jusqu'à une surcharge constante) avec le nombre d'itérations. Réduire le nombre d'itérations est donc essentiel pour accélérer l'algorithme. Bien que calculer à l'avance HARMONIC_10MIL soit une idée intelligente, cela conduit en fait à une moins bonne précision lorsque vous tronquez la série; le calcul d'une partie seulement de la série s'avère donner une plus grande précision.

Le code ci-dessous est une version modifiée du code posté ci-dessus (bien qu'en utilisant cython au lieu de numba ) .

t - log1p(t) = t ** 2 / 2 + ...

Il est capable d'obtenir une approximation proche même après 100 approximations comme le montre la figure ci-dessous.

À 100 itérations, son exécution est du même ordre de grandeur que scipy.special.gammaln:

log1p(t) = t - t ** 2 / 2 + ...

La question restante est bien sûr le nombre d'itérations à utiliser. La fonction log1p (t) peut être étendue comme une série de Taylor pour les petits t (ce qui est pertinent dans la limite des grands k ). En particulier,

%timeit special.gammaln(5)
# 932 ns ± 19 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit gammaln(5, 100)
# 1.25 µs ± 20.3 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

de telle sorte que, pour un grand k , l'argument de la somme devienne

from libc.math cimport log, log1p
cimport cython
cdef:
    float EULER_MAS = 0.577215664901532 # euler mascheroni constant

@cython.cdivision(True)
def gammaln(float z, int n=1000):
    """Compute log of gamma function for some real positive float z"""
    cdef:
        float out = -EULER_MAS*z - log(z)
        int k
        float t
    for k in range(1, n):
        t = z / k
        out += t - log1p(t)

    return out

Par conséquent, l'argument de la somme est nul jusqu'au second ordre dans t ce qui est négligeable si t est suffisamment petit. En d'autres termes, le nombre d'itérations doit être au moins aussi grand que z , de préférence au moins un ordre de grandeur plus grand.

Cependant, je m'en tiens à l'implémentation bien testée de scipy si possible.

J'ai réussi à obtenir une augmentation des performances d'environ 3x en essayant le mode parallèle de numba et en utilisant principalement des fonctions vectorisées (malheureusement, numba ne comprend pas numpy.substract.reduce )

#Your function:
%timeit gammaln(1.5)
48.6 ms ± 1.23 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

#My function:
%timeit gammaln_vec(1.5)
15 ms ± 340 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

#scpiy's function
%timeit gammaln_sp(1.5)
1.07 µs ± 18.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

Times:

from functools import reduce
import numpy as np
from numba import njit

@njit(fastmath=True, parallel=True)
def gammaln_vec(z):
    out = -EULER_MAS*z - np.log(z) + z*HARMONC_10MIL
    n = 10000000

    v = np.log(1 + z/np.arange(1, n+1))

    return out-reduce(lambda x1, x2: x1-x2, v, 0)

Alors encore, vous serez beaucoup mieux en utilisant la fonction de scipy. Sans code C, je ne sais pas comment le décomposer davantage

En ce qui concerne vos questions précédentes, je suppose qu'un exemple d'encapsulation des fonctions scipy.special dans Numba est également utile.

Exemple

Emballage Les fonctions cdef de Cython sont assez simples et portables tant qu'il n'y a que des types de données simples impliqués (int, double, double *, ...). Pour la documentation sur la façon d'appeler les fonctions scipy.special consultez ceci . Les noms de fonction dont vous avez réellement besoin pour envelopper la fonction se trouvent dans scipy.special.cython_special .__ pyx_capi__ . Les noms de fonction, qui peuvent être appelés avec différents types de données sont mutilés, mais déterminer le bon est assez facile (il suffit de regarder les types de données)

data=np.random.rand(1_000_000)
Test_func(A): 39.1ms
gammaln(A):   39.1ms

Utilisation dans Numba

#Numba example with loops
import numba as nb
import numpy as np
@nb.njit()
def Test_func(A):
  out=np.empty(A.shape[0])
  for i in range(A.shape[0]):
    out[i]=numba_gammaln(A[i])
  return out

Timings

#slightly modified version of https://github.com/numba/numba/issues/3086
from numba.extending import get_cython_function_address
from numba import vectorize, njit
import ctypes
import numpy as np

_PTR = ctypes.POINTER
_dble = ctypes.c_double
_ptr_dble = _PTR(_dble)

addr = get_cython_function_address("scipy.special.cython_special", "gammaln")
functype = ctypes.CFUNCTYPE(_dble, _dble)
gammaln_float64 = functype(addr)

@njit
def numba_gammaln(x):
  return gammaln_float64(x)

Bien sûr, vous pouvez facilement paralléliser cette fonction et surpassent l'implémentation de gammaln à thread unique dans scipy et vous pouvez appeler cette fonction efficacement dans n'importe quelle fonction compilée par Numba.