how to compute log factorial of an array of numbers

user_12627448 0
Viewed 548

The quantity to be computed is log(k!), where k could be 4000 or even higher, but of course the log will compensate. I tried computing sum(log(k)) which is the same.

So, I am given an large array with integers and I want to efficiently compute sum(log(k)). This was my attempt:

integers = np.asarray([435, 535, 242,])

score = np.sum(np.log(np.arange(1,integers+1)))

This would work, except that np.arange would generate an array of different size for each integer, so when I run that, it gives me an error (as it should).

The problem could be easily solved with a for loop as follows:

scores = []
for i in range(integers.shape[0]):
    score = np.sum(np.log(np.arange(1,integer[i]+1)))
    scores.append(score)

but that's too slow. My actual integers has millions of value to be computed.

Is there an efficient implementation for this that basically that doesn't need a for loop? I was thinking of a lambda function or something like that, but I am not really sure how to apply it. Any help is appreciated!

numpy arrays python performance
2 Answers

How about math.lgamma? Gamma function is factorial, and lgamma is log of gamma.

You don't need to compute factorial and then log.

There is also gammaln in the SciPy

Code, Python 3.9 x64 Win 10

import numpy as np
from scipy.special import gammaln

startf = 1 # start of factorial sequence
stopf  = 400 # end of of factorial sequence

q = gammaln(range(startf+1, stopf+1)) # n! = G(n+1)
print(q)

looks reasonable to me

user_40446960

You can vectorize with something like this:

mi = integers.max()
ls = np.log(np.arange(2, mi + 1))

Two optimizations so far: you only need the range up to the maximum, since the other numbers are covered by that, and you don't need log(1).

Now you take the cumulative sum:

cs = np.cumsum(ls)

The desired elements can be indexed directly:

result = cs[integers - 2]

If this is something you need to do many times, and you know the upper bound, this solution will be much faster than using math.lgmamma or scipy.special.gammaln once you precompute cs to the upper bound.

If this is a one-time call, here is the obligatory one-liner:

np.cumsum(np.log(np.arange(2, np.max(integers))))[integers - 2]

You can do most of the operations in-place if memory is a concern (I think it also makes them faster):

mi = integers.max()
cs = np.arange(2, mi + 1)
np.cumsum(np.log(cs, out=cs), out=cs)
user_29887300
Related