Let $Q_n$ be a $n\times n$ matrix with $Q_n=\begin{pmatrix} -\lambda_1-\mu_1 & \lambda_1 & 0 & \cdots\\ 0 & -\lambda_2-\mu_2 & \lambda_2 & \cdots\\ \vdots & \vdots & \vdots& \vdots\end{pmatrix}.$ In other words the diagonal of $Q$ is $\{-\lambda_i-\mu_i\}_{i=1}^n$ and superdiagonal is $\{\lambda_i\}_{i=1}^{n-1}$. Also, $\lambda_n=0$.

Letting $\mu=(\mu_1,\cdots,\mu_n)^T$, I need a fast way of calculating:

$$f(t):=p\exp(Qt)\mu,$$

where $t$ is a scalar and $p=(1,0,\cdots,0)$. I know there are a lot of matrix exponentiation algorithms but I'm hoping for one that could exploit the nice structure here. I would be happy with approximations, so long as there is good control of the maximum error. I was thinking of maybe trying to find a Jordan decomposition $Q=A+N$, where $A$ is diagonalizable and $N$ is nilpotent but this seems difficult.

Generally speaking, $Q_n$ will be diagonalizable as long as the diagonal entries are all distinct. However, I'm worried that when two values are close, there's a chance of high error.

Motivation: $f(t)$ is the density of a Coxian PhaseType distribution and for doing MLE it would be really nice to have a fast way of calculating $f(t)$.

isavailable, if you don't impose additional conditions. $\endgroup$ – Federico Poloni Dec 12 '15 at 11:50