Matrix Factorials

Specific Domains Numerics
3

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