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)$.
Best Answer
You probably want the Schur-Parlett method for computing matrix functions. It is a method to compute a generic function of a triangular matrix. Essentially, you apply the function to its diagonal elements and then use a recursion (derived from the identity $f(A)A=Af(A)$) to reconstruct the elements in positions $(i,i+1)$, then $(i,i+2)$ and so on, one superdiagonal at a time.
Trouble may ensue when two diagonal elements are close, as you correctly identify, but there are techniques to work around this issue.
See Chapters 9 and 10 of Higham's Functions of Matrices for a thorough treatment.