python - Factorial of a matrix elementwise with Numpy

Question

I'd like to know how to calculate the factorial of a matrix elementwise. For example,

import numpy as np
mat = np.array([[1,2,3],[2,3,4]])

np.the_function_i_want(mat)

would give a matrix mat2 such that mat2[i,j] = mat[i,j]!. I've tried something like

np.fromfunction(lambda i,j: np.math.factorial(mat[i,j]))

but it passes the entire matrix as argument for np.math.factorial. I've also tried to use scipy.vectorize but for matrices larger than 10x10 I get an error. This is the code I wrote:

import scipy as sp
javi = sp.fromfunction(lambda i,j: i+j, (15,15))
fact = sp.vectorize(sp.math.factorial)
fact(javi)

OverflowError: Python int too large to convert to C long

Such an integer number would be greater than 2e9, so I don't understand what this means.

Answer 1

There's a factorial function in scipy.misc which allows element-wise computations on arrays:

>>> from scipy.misc import factorial
>>> factorial(mat)
array([[  1.,   2.,   6.],
       [  2.,   6.,  24.]])

The function returns an array of float values and so can compute "larger" factorials up to the accuracy floating point numbers allow:

>>> factorial(15)
array(1307674368000.0)

You may need to adjust the print precision of NumPy arrays if you want to avoid the number being displayed in scientific notation.

---

Regarding scipy.vectorize: the OverflowError implies that the result of some of the calculations are too big to be stored as integers (normally int32 or int64).

If you want to vectorize sp.math.factorial and want arbitrarily large integers, you'll need to specify that the function return an output array with the 'object' datatype. For instance:

fact = sp.vectorize(sp.math.factorial, otypes='O')

Specifying the 'object' type allows Python integers to be returned by fact. These are not limited in size and so you can calculate factorials as large as your computer's memory will permit. Be aware that arrays of this type lose some of the speed and efficiency benefits which regular NumPy arrays have.