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.