# Matrix Factorials

Hi,

I’ve been looking for a function to calculate the factorial of a matrix, but haven’t had any luck.

I was wondering whether someone could point me in the right direction?

Apologies if i have missed something obvious!

Are you sure you want the factorial and not the exponent? or perhaps the factorization (i.e. representation as a product of some structures matrices)?

I guess you mean to extend factorial to matrices by means of the gamma function (stackoverflow)?

Is your matrix symmetric positive definite? If so, you could simple apply the gamma function to the eigenvalues, e.g.

``````julia> using LinearAlgebra, SpecialFunctions

julia> A = randn(5,5); A = Symmetric(A'A) # random PSD matrix
5×5 Symmetric{Float64, Matrix{Float64}}:
2.91012   -0.227897  -1.31028   -0.198762   0.375903
-0.227897   3.00476   -3.01805   -1.45998   -1.30752
-1.31028   -3.01805    4.35062    1.72612    0.932463
-0.198762  -1.45998    1.72612    1.79473   -0.824075
0.375903  -1.30752    0.932463  -0.824075   2.65364

julia> F = eigen(A);

julia> Afact = F.vectors * Diagonal(gamma.(F.values .+ 1)) * F.vectors' # "factorial" of A
5×5 Matrix{Float64}:
1047.56    3786.34   -4796.89  -2069.53  -1358.66
3786.34   13823.4   -17490.0   -7541.97  -4976.91
-4796.89  -17490.0    22135.4    9545.66   6294.27
-2069.53   -7541.97    9545.66   4119.14   2711.5
-1358.66   -4976.91    6294.27   2711.5    1798.54
``````

(More generally, you could use diagonalization for other sorts of matrices, but for non-Hermitian matrices you can run into numerical problems for matrices that are nearly defective. For non-PSD Hermitian matrices you may also run into poles of the gamma function for negative eigenvalues, but of course those will show up no matter what computational method you use. The gamma function is nicest for PSD matrices.)

Or do you mean the elementwise factorial? That’s just `gamma.(A .+ 1)` in Julia, or `factorial.(A)` for integer matrices.

1 Like

Hi @stevengj , yes that is exactly what I am looking for, a much simpler implementation than I expected and the StackExchange link you sent is a great explanation. Thank you very much!