9
$\begingroup$

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

$\endgroup$
3
  • 1
    $\begingroup$ You already have a trivial diagonal + nilpotent decomposition, by separating the diagonal and the superdiagonal part. If you want one with all ones in the superdiagonal, just multiply by the matrix $\operatorname{diag}(\lambda_1,\lambda_2,\dots,\lambda_n)$ and its inverse on both sides. $\endgroup$ Dec 12, 2015 at 7:14
  • 1
    $\begingroup$ @Federico Poloni: unfortunately you also need the resulting pair of matrices to commute which isn't the case with that decomposition $\endgroup$
    – Alex R.
    Dec 12, 2015 at 9:38
  • $\begingroup$ I agree, this isn't going to be useful to compute an exponential. I was just pointing out that a decomposition of the kind that you mention is available, if you don't impose additional conditions. $\endgroup$ Dec 12, 2015 at 11:50

3 Answers 3

5
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ It looks like this is exactly the (non-blocked version of) the Schur-Parlett method that I was also mentioning. Interesting connection. $\endgroup$ Dec 12, 2015 at 11:55
3
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.