# How could we define the factorial of a matrix?

Suppose I have a square matrix $\mathsf{A}$ with $\det \mathsf{A}\neq 0$.

How could we define the following operation? $$\mathsf{A}!$$

Maybe we could make some simple example, admitted it makes any sense, with

$$\mathsf{A} = \left( \begin{matrix} 1 & 3 \\ 2 & 1 \end{matrix} \right)$$

-
Made me wonder: do you ask just out of curiosity (i.e. recreational mathematics), or do you have some application/usage in mind? – Oliphaunt Jan 31 at 21:19
I asked this as a separate question: math.stackexchange.com/questions/1637318/… – Oliphaunt Feb 2 at 12:21

For any holomorphic function $G$, we can define a corresponding matrix function $\bar{G}$ via (a formal version of) the Cauchy Integral Formula: We set $$\bar{G}(B) := \frac{1}{2 \pi i} \oint_C G(z) (z I - B)^{-1} dz ,$$ where $C$ is an (arbitrary) anticlockwise curve that encloses the eigenvalues of the (square) matrix $B$. Note that the condition on $C$ means that restrictions on the domain of $G$ determine restrictions on the domain of $\bar{G}$.

Now, for all nonnegative integers $n$, $n!$ coincides with the value $F(n)$ of the holomorphic function $$F: z \mapsto \Gamma(z + 1)$$ at $n$, where $\Gamma$ denotes the Gamma function. So, we can make sense of the factorial of a (square) matrix $B$ by forming $\bar{F}$ as above and declaring $B! := \bar{F}(B)$.

The Cauchy Integral Formula definition of $\bar G$ coincides, by the way, with the formal evaluation $\sum_{i = 0}^{\infty} a_k (B - z_0 I)^k$ of the power series $\sum_{i = 0}^{\infty} a_k (z - z_0)^k$ for $G$, at least on a set where the series has the appropriate convergence behavior, but NB this set may be empty! Indeed, I don't believe that power series for $F$ can be used directly to evaluate $A!$ for our particular $A$: The function $F$ has a pole on the line segment in $\Bbb C$ with endpoints the eigenvalues of $A$, so there is no basepoint $z_0$ for which the corresponding series converges at $A$.

We can define $B!$ in another way that coincides appropriately with Cauchy Integral Formula definition and is more amenable to explicit computation: If $B$ is diagonalizable, so that we can decompose $$B = P \pmatrix{\lambda_1 & & \\ & \ddots & \\ & & \lambda_n} P^{-1} ,$$ for eigenvalues $\lambda_a$ of $B$ and some matrix $P$, we define $$\bar{G}(B) := P \pmatrix{G(\lambda_1) & & \\ & \ddots & \\ & & G(\lambda_n)} P^{-1} .$$ Indeed, by substituting and rearranging, we can see that this coincides, at least formally, with the power series characterization. (There is a similar but more complicated form for nondiagonalizable $B$ that I won't write out here, but which is given in the Wikipedia article Matrix function; we won't need it here anyway.)

In our example, $A$ has distinct eigenvalues $\lambda_{\pm} = 1 \pm \sqrt{6}$, and can so can be diagonalized as $$P \pmatrix{1 - \sqrt{6} & 0 \\ 0 & 1 + \sqrt{6}} P^{-1} ;$$ indeed, we can take $$P = \pmatrix{\tfrac{1}{2} & \tfrac{1}{2} \\ -\frac{1}{\sqrt{6}} & \frac{1}{\sqrt{6}}}.$$

Now, $F(\lambda_{\pm}) = \Gamma(\lambda_{\pm} + 1) = \Gamma (2 {\pm} \sqrt{6}),$ and putting this all together gives that \begin{align*}\pmatrix{1 & 3 \\ 2 & 1} ! = \bar{F}(A) &= P \pmatrix{F(\lambda_-) & 0 \\ 0 & F(\lambda_+)} P^{-1} \\ &= \pmatrix{\tfrac{1}{2} & \tfrac{1}{2} \\ -\frac{1}{\sqrt{6}} & \frac{1}{\sqrt{6}}} \pmatrix{\Gamma (2 - \sqrt{6}) & 0 \\ 0 & \Gamma (2 + \sqrt{6})} \pmatrix{1 & -\frac{\sqrt{3}}{\sqrt{2}} \\ 1 & \frac{\sqrt{3}}{\sqrt{2}}} .\end{align*} Multiplying this out gives $$\color{#bf0000}{\boxed{\pmatrix{1 & 3 \\ 2 & 1} ! = \pmatrix{\frac{1}{2} \alpha_+ & \frac{\sqrt{3}}{2 \sqrt{2}} \alpha_- \\ \frac{1}{\sqrt{6}} \alpha_- & \frac{1}{2} \alpha_+}}} ,$$ where $$\color{#bf0000}{\alpha_{\pm} = \Gamma(2 + \sqrt{6}) \pm \Gamma(2 - \sqrt{6})}.$$ (Optionally, we can use the factorial-like identity $\Gamma(z + 1) = z \Gamma(z)$ to write $$\Gamma(2 \pm \sqrt{6}) = (1 \pm \sqrt{6}) \Gamma(1 \pm \sqrt{6}) = (6 \pm \sqrt{6}) \Gamma(\pm \sqrt{6})$$ and the reflection formula $$-z \Gamma(z) \Gamma(-z) = \frac{\pi}{\sin \pi z}$$ to write the entries as expressions algebraic in $\pi$, $\sin(\pi \sqrt{6})$, and $\Gamma(\sqrt{6})$ alone.) It's perhaps not very illuminating, but $A!$ has numerical value $$\pmatrix{1 & 3 \\ 2 & 1}! \approx \pmatrix{3.62744 & 8.84231 \\ 5.89488 & 3.62744} .$$

To carry out these computations, I wrote the following simple Maple procedure (in particular it implements the above formula and hence only works in the diagonalizable case):

with(LinearAlgebra):

MatrixFactorial := proc(A) local P, lambda;
P      := JordanForm(A, output = 'Q');
lambda := Diagonal(JordanForm(A, output = 'J'));
DiagonalMatrix(map(u -> GAMMA(u + 1), lambda));
P.%.MatrixInverse(P);
end proc;


After executing the above, one replicate the computation of $A!$ by running

A := Matrix([[1, 3], [2, 1]]);
MatrixFactorial(A);
evalf(%);


By the way, we need not have that $\det B \neq 0$ in order to define $B!$. In fact, proceeding as above we find that the factorial of the zero matrix is $$0! = I .$$ Likewise using the formula for nondiagonalizable matrices referenced above together with a special identity gives that the factorial of the $2 \times 2$ Jordan block of eigenvalue zero is, somewhat amusingly, $$\pmatrix{0 & 1\\0 & 0} ! = \pmatrix{1 & -\gamma \\ 0 & 1} ,$$ where $\gamma$ is the Euler-Mascheroni constant.

-
+1 because thou took time to answer in a very useful way! Thank you! – Time Master Jan 31 at 15:53
A truly excellent answer. – goblin Jan 31 at 15:57
@YoTengoUnLCD One can use, e.g., {\Large !} to increase the size of the factorial symbol, but the way MathJax aligns elements vertically makes this look strange for font sizes as large as you might like them. A kludge for forcing the vertical alignment is embedding the factorial symbol in a (bracketless) matrix, with something like \pmatrix{a&b\\c&d}\!\!\matrix{\Huge !}, which produces $$\pmatrix{a&b\\c&d}\!\!\matrix{\Huge !}$$ The commands \! are used to improve the kerning. – Travis Feb 1 at 11:06
@KimPeek *thou took'st, if I'm not mistaken :-D – The Vee Feb 1 at 15:44
@TobiasKienzler The function $\Gamma$ is holomorphic on its domain, and indeed, this much is necessary to guarantee path-independence of the integral in the definition of $\bar{F}$. This is also why we need to enclose all of the eigenvalues: The integrand $F(z) (z I - B)^{-1}$ has poles at the eigenvalues of $B$, and in general these poles contribute to the resulting integrand, so an integral over some loop not enclosing all the eigenvalues of $B$ will simply give a value other than $\bar{F}(B)$. – Travis Feb 2 at 11:59

The gamma function is analytic. Use the power series of it.

EDIT: already done: Some properties of Gamma and Beta matrix functions (maybe paywalled).

-
(+1) for mentioning this technique, but I don't believe it's possible to use power series to compute $A!$ for the particular example matrix $A$ given: The line segment connecting the eigenvalues of $A$ contains a pole of $z \mapsto \Gamma(z + 1)$, so no power series for that function (i.e., for any basepoint) converges at $A$. – Travis Jan 31 at 14:23
This would not have occurred to me! But the issue of convergence is a complicated one, I think. The gamma function has an infinite number of poles, after all, so it doesn't have a power series valid everywhere. – TonyK Jan 31 at 14:23
@Travis, obviously the convergence will depend of the matrix. – Martín-Blas Pérez Pinilla Jan 31 at 14:29
@VašekPotoček I don't think that's true; since the function $x \mapsto \Gamma(x + 1)$ is well-behaved at the eigenvalues of the given matrix $A$, I believe we can use the Jordan Canonical Form to make sense of $A!$. See my answer for more---comments and corrections are most welcome! – Travis Jan 31 at 14:51
@Martín-BlasPérezPinilla Yes, my above comment was restricted to the example in the question. But the same reasoning shows that the relevant power series will not converge (again, for any base point) for a large open set of matrices: I think this is the case, for example, if a matrix has an eigenvalue $\lambda$ with $\Re \lambda < -1$ and $|\Im \lambda| > \frac{1}{2}$. – Travis Jan 31 at 15:56

I don't have enough reputation points to comment on Travis' answer, but his numerical result is incorrect.

Using Julia I get

A = [1 3;2 1]
EVD = eigfact(A)
V = EVD[:vectors]
g = gamma(EVD[:values])
gammaA = V * diagm(g) * inv(V)
factA = A * gammaA

3.62744  8.84231
5.89488  3.62744



As long as cond(V) isn't too terrible, I've found the above procedure to be a practical way to evaluate arbitrary matrix functions.

-
Thanks for pointing this out. There was a bug in my Maple code, and it affected the exact value, too; I've since corrected both values in my point. – Travis Feb 1 at 10:37

I would start from the logical definition of the matrix factorial, without assuming that we want to cover all properties that we know from factorial in set of reals.

We define standard factorial as $1 \cdot (1+1) \cdot (1+1+1) \cdot ... \cdot (1+1+...+1+1)$

So first let us define $[n]!$ using the same logic replacing 1 with identity matrix. The obvious way to define it is

$$[n]!=\prod\limits_{k=1}^{n}\begin{bmatrix} k & 0\\ 0 & k \end{bmatrix}=\begin{bmatrix} n! & 0\\ 0 & n! \end{bmatrix}$$

All properties of the standard factorial are there. Now, we were defining Gamma function by simple extension $\Gamma (x+1)=x\Gamma (x)$ where $n!=\Gamma (n+1)$. That is all that is required. So we want to find matrix Gamma $\Gamma ([x]+I)=[x]\Gamma ([x])$

If we define

$$\Gamma (\begin{bmatrix} x & 0\\ 0 & x \end{bmatrix})=\begin{bmatrix} \Gamma (x) & 0\\ 0 & \Gamma (x) \end{bmatrix}$$

we are totally fine because

$$\begin{bmatrix} x & 0\\ 0 & x \end{bmatrix}\begin{bmatrix} \Gamma (x) & 0\\ 0 & \Gamma (x) \end{bmatrix}=\begin{bmatrix} x\Gamma (x) & 0\\ 0 & x\Gamma (x) \end{bmatrix}=\begin{bmatrix} \Gamma (x+1) & 0\\ 0 & \Gamma (x+1) \end{bmatrix}$$

There is nothing amiss if we start from $\begin{bmatrix} x & 0\\ 0 & y \end{bmatrix}$ because

$$\begin{bmatrix} x & 0\\ 0 & y \end{bmatrix}\begin{bmatrix} \Gamma (x) & 0\\ 0 & \Gamma (y) \end{bmatrix}=\begin{bmatrix} x\Gamma (x) & 0\\ 0 & y\Gamma (y) \end{bmatrix}=\begin{bmatrix} \Gamma (x+1) & 0\\ 0 & \Gamma (y+1) \end{bmatrix}$$

The remaining part is the other diagonal. What to do with $A=\begin{bmatrix} x_{0} & x_{1}\\ x_{2} & x_{3} \end{bmatrix}$?

So we start from what we would like to have $\Gamma([A]+I)=[A]\Gamma([A])$.

If we are able to diagonalize $A=P^{-1}\overline{A}P$ and to express in the same manner $\Gamma([A]) = P^{-1}\Gamma(\overline{A})P$ then we have

$$\Gamma([A]+I) = P^{-1} \overline{A} P P^{-1} \Gamma(\overline{A}) P = P^{-1} \overline{A} \Gamma(\overline{A}) P = P^{-1} \Gamma(\overline{A+I}) P=\Gamma(A+I)$$

so all should be fine.

Since $\overline{A}$ is diagonal with eigenvalues on the main diagonal $\lambda_{1},\lambda_{2}$ and we know how to deal with that type of matrix, we have the full definition of $\Gamma(A)$ even for matrices.

$$\Gamma(A)=P^{-1}\begin{bmatrix} \Gamma (\lambda_{1}) & 0\\ 0 & \Gamma (\lambda_{2}) \end{bmatrix}P$$

and now $A!=\Gamma(A+I)$ making it all

$$A!=P^{-1}\begin{bmatrix} \Gamma (\lambda_{1}+1) & 0\\ 0 & \Gamma (\lambda_{2}+1) \end{bmatrix}P$$

Instead of giving the solution just to the example I will give a general form for 2x2 matrix $\begin{bmatrix} a & b \\ c & d \end{bmatrix}$. Take discriminant $D=\sqrt{(a-d)^2+4bc} \neq 0, c \neq 0$. Then

$$\begin{bmatrix} a & b \\ c & d \end{bmatrix} ! = \begin{bmatrix} \frac{a-d-D}{2c} & \frac{a-d+D}{2c} \\ 1 & 1 \end{bmatrix} \begin{bmatrix} \Gamma (\frac{a+d-D}{2}+1 ) & 0 \\ 0 & \Gamma ( \frac{a+d+D}{2} +1)\end{bmatrix} \begin{bmatrix} -\frac{c}{D} & \frac{a-d+D}{2D} \\ \frac{c}{D} & -\frac{a-d-D}{2D} \end{bmatrix}$$

From here you can nicely conclude that the factorial matrix can be expressed using classical integer factorial if $a+d \pm D$ are even positive integers (including $0$).

For other values we use the extension of $\Gamma(x)$ itself.

-
Very nice! Could you fix the two little "bugs" in the LaTeX code? ^^ Then I'll read it with pleasure! – Time Master Feb 3 at 19:54
@KimPeek: still editing, I have to look for myself how it looks. I think it looks fine now – user195934 Feb 3 at 19:55
@AlexPeter: If you can add the factorial of the matrix which is given in the question that would be awesome. Any way very nice answer. ++1 – Nil Mal Feb 9 at 7:09

It would be good to mention that (almost) any matrix function can be made into a power-series expansion, which eventually involves the values of the function on the eigenvalues of the matrix multiplied by the eigenvectors.

In other words the matrix function is completely characterised by the values it takes on the eigenvalues of the matrix (even if a power-series expansion may be needed).

The above hold for matrices which are diagonalisable (i.e. the number of linearly independent eigenvectors is equal to the matrix dimension). There are ways to expand an arbitrary matrix into what is referred to as generalised eigenvectors, but this will not be pursued further here.

Furthermore, since any square, finite-dimensional, matrix satisfies its characteristic polynomial equation (if seen as a matrix function), aka Cayley-Hamilton theorem, the powers of $A^k$ for $k \ge n$ ($n$ is the dimension) can be expressed as a function of the powers of $A$ up to $n$. So eventually the matrix function power-series expansion collapses to polynomial expansion (for square matrices). Finally, this polynomial expansion, for a given function, can be found more easily by methods such as variation of parameters or polynomial modeling.

-
"any matrix function can be made into a power-series expansion, which eventualy involves the values of the function on the eigen-values of the matrix multiplied by the eigen-vectors." "every matrix function is completely characterised by the values it takes on the eigen-values of the matrix" Both statements are wrong for nondiagonalizable matrices. – Did Feb 4 at 6:42
Your remark about full-rank or non full-rank matrices is completely offtopic. Please refer to some simple examples of nondiagonalizable matrices to check if the statements in your answer are valid for such matrices (they are not). – Did Feb 4 at 12:30
Full-rank or not full-rank $\ne$ Diagonalisable or not. Please refer to some textbook on matrices. (How come that "eigenvectors" in your first comment mutate to "generalised eigenvectors"? Is this some kind of rhetorical trick? Are "full-rank" and "diagonalisable" supposed to become "generalized full-rank" and "generalized diagonalizable"?) – Did Feb 4 at 18:39
@Did, i'm not sure who is trying to use rhetorical tricks here. Anyway "The diagonalization theorem states that an n×n matrix A is diagonalizable if and only if A has n linearly independent eigenvectors, i.e., if the matrix rank of the matrix formed by the eigenvectors is n", i hope you do like Wolfram, as for the rest they are not touched upon. If you think the phrasing can be made better, no problem, else it is beating around the bush and poor use of my time – Nikos M. Feb 4 at 21:18
@NikosM. "Full-rank" does not mean that the eigenvectors have full span. Please revise what the rank of a matrix is. – Did Feb 5 at 6:46