Fast Upper Triangular Matrix Exponentiation

Question

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)$.

Answer 1

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.

Answer 2

When I teach the exponential of matrices, I tell the students that the converging series is not a practical tool for calculation. It is way better to solve the differential equation. This turns out to be true here. Say $T$ is upper triangular and $M$ denotes $\exp T$. For the sake of simplicity, set $t_i=t_{ii}$Then one finds $$m_{ii}=e^{t_{i}},\qquad m_{i-1,i}=t_{i-1,i}\frac{e^{t_{i}}-e^{t_{i-1}}}{t_{i}-t_{i-1}}$$ and $$m_{i-2,i}=\frac{t_{i-2,i}t_{i-1,i}}{t_{i}-t_{i-2}}\left(\frac{e^{t_{i}}-e^{t_{i-2}}}{t_{i}-t_{i-2}}-\frac{e^{t_{i-1}}-e^{t_{i-2}}}{t_{i-1}-t_{i-2}}\right)$$ and so on. All the entries $m_{ij}$ are obtained by successive finite differences. As usual, when an argument is repeated, the finite difference is a derivative.

Answer 3

The exponential $e^{Q}$ of any $n\times n$ upper triangular matrix $Q$ can be computed efficiently by solving a set of $n$ first-order differential equations, $u_{i}'(t)=\sum_{j}Q_{ij}u_j(t)$; these $n$ equations can be solved one by one: solve first for $i=n$ and then work back to $i=n-1,n-2,\ldots 1$; at each step you have a first-order linear ODE of the form $y'(t)+ay(t)+b(t)=0$, that can be readily solved with the method of integrating factors. The $j$-th column of $e^{Q}$ is the solution $u(1)$ with $u_i(0)=\delta_{ij}$.