Computing the exponential of an intensity matrix (Q1801855)

The probability state vector of a discrete-state continuous-time Markov process can be expressed as \(\pi(t) = e^{Qt}\pi(0)\), where \(Q\) is the intensity matrix satisfying \(q_{ij} \geq 0\) for all \(i \neq j\) and \(\sum^ n_{j=1}q_{ij} = 0\) for all \(i\). The scaling and squaring method for computing the matrix exponential is based on the identity \(e^{Qt} = e^{-qt}[e^{2^{-m}Q't}]^{2^ m}\), where \(Q' = Q+qI\) and the number \(m\) is chosen so that \(2^{m-1} \leq \| Q't\| < 2^ m\). The number \(q\) is chosen so that \(\| Q'\|\) is minimized and can be computed exactly for an intensity matrix. With this preconditioning, the exponential \(e^{2^{-m}Q't}\) can be computed using either Padé or Taylor approximants, and the desired matrix \(e^{Qt}\) by successive squaring of the result. This procedure circumvents many of the reliability problems associated with the standard procedure.