# Transition probability matrix rows not summing to 1

I'm implementing a method of sampling from a CTMC from here.

I'm trying to calculate the transition probability matrix but the rows are not summing to 1, except at 0.

I've diagonalised my rate matrix $Q$ to give me a matrix $D$ such that $Q=UDU^{-1}$ where $U$ is the matrix with it's columns the eigenvectors of $Q$ and $D$ the diagonal matrix with diagonal entries the eigenvalues of $Q$.

I should then be able to calculate the transition probability matrix of the CTMC using the equation $$P(t)=e^{Qt}=Ue^{tQ}U^{-1}.$$

I've implemented this using R, with the following code:

eig <- eigen(rateM)
U <- eig$vectors invU <- solve(U) Pt <- U%*%diag(exp(eig$values*t))%*%invU


where

> rateM
[,1]       [,2]       [,3]       [,4]
[1,] -1.1224490  0.6122449  0.3061224  0.2040816
[2,]  0.4081633 -0.9183673  0.3061224  0.2040816
[3,]  0.2040816  0.3061224 -0.9183673  0.2040816
[4,]  0.2040816  0.3061224  0.6122449 -1.1224490
> t
[1] 3


The rows of the rate matrix all sum to zero, $UU^{-1}=I$ as it should, but the larger t gets, the further from 1 Pt gets. E.g.

> sum((U%*%diag(exp(eig$values*0))%*%invU)[1,]) [1] 1 > sum((U%*%diag(exp(eig$values*1))%*%invU)[1,])
[1] 0.9773404
> sum((U%*%diag(exp(eig$values*2))%*%invU)[1,]) [1] 0.9316814 > sum((U%*%diag(exp(eig$values*3))%*%invU)[1,])
[1] 0.8799868
> sum((U%*%diag(exp(eig$values*10))%*%invU)[1,]) [1] 0.5701228  Is this owing to some sort of numerical error, or am I doing something incorrectly? • Your rateM does not have rows that sum to zero. Even though the sums are close to zero, the exponentiation is making those little differences blow up. – whuber Commented Jan 8, 2016 at 1:12 • Following on to whuber. The 3rd row of rateM sums to -0.2040817 . That's not good. Commented Jan 8, 2016 at 1:57 • D'oh. I made an error when building my rate matrix. I corrected it and this now works as it should! Thanks! Commented Jan 8, 2016 at 10:27 ## 1 Answer Converting comment to answer, since the comment was on the money, thereby allowing this to be closed. Credit to whuber for his initial comment. The 3rd row of rateM sums to -0.2040817, not to zero, and is therefore erroneous. Therefore the transition probability matrix will be erroneous. • +1. I would add a remark about numerical stability. When the rows sum to zero, the matrix is prima facie singular and must have a zero eigenvalue. The computation of eigen may render that zero as a small number. If, however, any such small number is replaced by a true zero (as in zapsmall in R), then its exponential will truly be$1$, thereby guaranteeing the sum-to-unity property of the matrix exponential (up to floating point roundoff in the components of$U$) no matter what value$t\$ might have.
– whuber
Commented Jan 8, 2016 at 14:28