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!