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

Question

Я пытаюсь написать быстрый алгоритм для вычисления логарифмической функции log .В настоящее время моя реализация кажется наивной и просто повторяется 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 реализация и моя буквально в 100 000 раз медленнее.Так что я делаю что-то очень неправильное или очень наивное (вероятно, оба).Хотя мои ответы, по крайней мере, правильны с точностью до 4 десятичных знаков в худшем случае по сравнению со scipy.

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

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

Answer 1

Мне удалось увеличить производительность примерно в 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)

Times:

#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 я не знаю, как его разбить дальше

Answer 2

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

Пример

Обтекание функций Cython cdef довольно простои переносимый, если задействованы только простые типы данных (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

#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функция.

Answer 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 достаточно мало.Другими словами, число итераций должно быть не менее 1032 *, предпочтительно на порядок больше.

Однако я бы придерживался проверенной реализации scipy.если это вообще возможно.