Recently I needed to build a Monte Carlo simulator of a continuous-time Markov chain. This is a pretty straightforward exercise; the only catch was that I wanted it to perform well, so I had to use a fast algorithm for matrix exponentiation.
One way to calculate a matrix exponent efficiently is to calculate the eigendecomposition of the matrix: , where is a diagonal matrix. Then (details here). The process of finding such eigendecomposition is also known as diagonalisation.
Not every matrix can be diagonalised, though. One particularly nice case is when is Hermitian: then there exists a decomposition where is real, and is unitary. Numerical diagonalisation of a Hermitian matrix is implemented in many software packages.
Even if the matrix is question is not Hermitian, not all is lost: there is a big class of non-Hermitian matrices that can be easily diagonalised. This is what this post is about.
If we can find matrix such that
- is Hermitian,
- can be decomposed as ( denotes conjugate transpose of ), where is invertible,
then there exists a decomposition where is real. Note: in this post, I discuss a general case where elements of are complex numbers. All of that is also applicable, as a special case, to real-valued matrices. For them, “Hermitian” means “symmetrical”, and “unitary” means “orthogonal”. The elements of an intensity matrix of a Markov chain are, of course, real.
To find the eigendecomposition of a non-Hermitian matrix , we start with the eigendecomposition of the matrix . We use abbreviation . Matrix is Hermitian, which is easy to see by checking that using the fact that .
Suppose we have found the eigendecomposition
(using our favourite matrix manipulation package). Then, by definition of , we have
Multiplying both sides of this equation by on the left and on the right, we get
Let , then . Thus,
which is the eigendecomposition we were looking for.