Быстрый алгоритм для лог-гамма-функции

Я пытаюсь написать быстрый алгоритм для вычисления лог-гамма-функции. В настоящее время моя реализация кажется наивной и просто повторяется 10 миллионов раз для вычисления журнала гамма-функции (я также использую numba для оптимизации кода).

import numpy as np
from numba import njit
EULER_MAS = 0.577215664901532 # euler mascheroni constant
HARMONC_10MIL = 16.695311365860007 # sum of 1/k from 1 to 10,000,000

@njit(fastmath=True)
def gammaln(z):
"""Compute log of gamma function for some real positive float z"""
    out = -EULER_MAS*z - np.log(z) + z*HARMONC_10MIL
    n = 10000000 # number of iters
    for k in range(1,n+1,4):
        # loop unrolling
        v1 = np.log(1 + z/k)
        v2 = np.log(1 + z/(k+1))
        v3 = np.log(1 + z/(k+2))
        v4 = np.log(1 + z/(k+3))
        out -= v1 + v2 + v3 + v4

    return out

Я сравнил свой код с реализацией scipy.special.gammaln, а мой буквально в 100000 раз медленнее. Так что я делаю что-то очень неправильное или очень наивное (вероятно, оба). Хотя мои ответы, по крайней мере, верны с точностью до 4 десятичных знаков в худшем случае по сравнению со скупым.

Я пытался прочитать код _ufunc, реализующий функцию gipmal от scipy, однако я не понимаю код cython, в котором написана функция _gammaln.

Есть ли более быстрый и оптимизированный способ расчета гамма-функции журнала? Как я могу понять реализацию scipy, чтобы я мог включить ее в свою?

Источник

3 ответа

Решение

Время выполнения вашей функции будет масштабироваться линейно (до некоторой постоянной нагрузки) с количеством итераций. Таким образом, сокращение количества итераций является ключом к ускорению алгоритма. В то время как вычисление HARMONIC_10MIL заранее умная идея, она на самом деле приводит к худшей точности, когда вы усекаете ряд; Вычисление только части серии дает более высокую точность.

Код ниже является модифицированной версией кода, опубликованного выше (хотя с использованием cython вместо numba).

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

Он может получить близкое приближение даже после 100 приближений, как показано на рисунке ниже.

На 100 итерациях его время выполнения того же порядка, что и scipy.special.gammaln:

%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)

Остающийся вопрос, конечно, сколько итераций использовать. Функция log1p(t) может быть расширен как серия Тейлор для малых t (что актуально в пределе большого k). Особенно,

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

такой, что для большого k, аргумент суммы становится

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

Следовательно, аргумент суммы от нуля до второго порядка в t что незначительно, если t достаточно маленький. Другими словами, число итераций должно быть не менее z предпочтительно по меньшей мере на порядок больше.

Тем не менее, я бы придерживался scipy Хорошо проверенная реализация, если это вообще возможно.

Источник

Мне удалось добиться увеличения производительности примерно в 3 раза, попробовав параллельный режим numba и используя в основном векторизованные функции (к сожалению, numba не может понять numpy.substract.reduce)

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)

Времена:

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

Тем не менее, вы будете гораздо лучше, используя функцию Сципи. Без кода на C я не знаю, как его сломать

Источник

Что касается ваших предыдущих вопросов, я думаю, что пример оборачивания scipy.special функции для Numba также полезны.

пример

Обертывание функций cdef в Cython довольно простое и переносимое, если в нем задействованы только простые типы данных (int, double, double*,...). Для документации о том, как вызывать функции scipy.special , посмотрите на это. Имена функций, которые вам действительно нужны для переноса функции, находятся в scipy.special.cython_special.__pyx_capi__, Названия функций, которые можно вызывать с разными типами данных, искажены, но определить правильный довольно легко (достаточно взглянуть на типы данных)

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

Использование в Нумбе

#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

Задержки

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

Конечно, вы можете легко распараллелить эту функцию и превзойти однопоточную реализацию gammaln в scipy, и вы можете эффективно вызывать эту функцию в любой скомпилированной функции Numba.

Источник
Другие вопросы по тегам