Why Scale Before Computing an Exponential?
For a finite-state continuous-time Markov chain, the transition matrix at time \(t\) is \(P(t)=e^{tQ}\), where \(Q\) is the generator. For a large or awkward matrix, directly summing the exponential series can be inefficient or inaccurate. Scaling and squaring addresses this by first computing an exponential of a smaller matrix, then recovering the desired result through repeated squaring.
Let \((X_t)_{t\ge0}\) be a continuous-time Markov chain on a finite state space \(S\), defined on a probability space \((\Omega,\mathcal{F},\mathbb{P})\), and let \(Q=(q_{ij})\) be its generator. The off-diagonal entries of \(Q\) are nonnegative and every row sums to zero. For \(t\ge0\), \(P(t)=e^{tQ}\) has entries \(p_{ij}(t)=\mathbb{P}(X_t=j\mid X_0=i)\). Thus \(P(t)\) is row-stochastic: its entries are nonnegative and each row sums to one. The computation below approximates this matrix; it does not change the definition of the transition probabilities.
A simple choice for \(R\), useful for learning the method, is the degree-\(m\) Taylor polynomial \(R_m(B)=I+B+B^2/2!+\cdots+B^m/m!\). Practical numerical libraries often use Padé approximants with scaling and squaring instead; the choice of approximation and thresholds is implementation-dependent. The worked calculations here use Taylor polynomials so each step can be checked by hand.
Why the Procedure Works, and What Controls Its Error
Proof. Put \(C=A/2^s\) and \(N=2^s\). The exponential series is absolutely convergent for every finite matrix. Multiplying the series for \(e^C\) by itself and collecting powers of \(C\) gives \(e^C e^C=e^{2C}\): the coefficient of \(C^k\) in the product is \(\sum_{r=0}^k 1/(r!(k-r)!)=2^k/k!\). Repeating this identity \(s\) times gives \((e^C)^N=e^{NC}=e^A\). For the bound, absolute convergence and the triangle inequality give \(\|e^B-R_m(B)\|_\infty\le\sum_{k=m+1}^\infty\|B\|_\infty^k/k!\), because the induced infinity norm is submultiplicative. Write \(k=m+1+r\). Since \((m+1+r)!\ge(m+1)!r!\), each term is at most \(\|B\|_\infty^{m+1}\|B\|_\infty^r/((m+1)!r!)\). Summing over \(r\ge0\) gives the stated bound. \(\square\)
The bound explains the role of scaling: reducing \(\|B\|_\infty\) can make a fixed-degree Taylor approximation more accurate. It is a bound on the truncation error before squaring, not a guarantee that all roundoff errors are negligible. Squaring also propagates the approximation error. For example, if \(R\) approximates \(E=e^{A/2^s}\), then the algebraic identity \(R^N-E^N=\sum_{k=0}^{N-1}R^k(R-E)E^{N-1-k}\) gives $$ \|R^N-E^N\|_\infty\le N\max(\|R\|_\infty,\|E\|_\infty)^{N-1}\|R-E\|_\infty. $$
To justify this last inequality, expand the sum: consecutive terms cancel, leaving \(R^N-E^N\). Apply the triangle inequality and submultiplicativity to each summand. Its norm is at most \(\max(\|R\|_\infty,\|E\|_\infty)^{N-1}\|R-E\|_\infty\); there are \(N\) summands. Thus a small initial error is useful, but it should not be treated as an exact answer.
Worked Example: The Zero Generator
Consider a chain with one state, \(S=\{1\}\), and generator \(Q=[0]\). Its sole row sums to zero, as required. For any \(t\ge0\), \(A=tQ=[0]\), so \(\|A\|_\infty=0\). There is no need to scale: take \(s=0\). The Taylor approximation of any degree \(m\ge0\) is \(R_m(A)=I=[1]\), since all positive powers of the zero matrix are zero. With \(2^0=1\), the result is \(P(t)=[1]\).
Check: The only row sum is \(1\), and the entry is nonnegative. Also \(e^{tQ}=e^{[0]}=[1]\), so the computed transition probability is exact. This case is a useful boundary check: if the matrix norm is zero, scaling cannot improve the calculation and is unnecessary.
A Two-State Calculation
Worked Example: A Symmetric Switching Chain
Let \(S=\{1,2\}\), with rate \(1\) from either state to the other. The generator is $$ Q=\begin{bmatrix}-1&1\\1&-1\end{bmatrix}. $$
Each generator row sums to zero: \(-1+1=0\) in the first row and \(1-1=0\) in the second. Compute the transition matrix at \(t=1\), so \(A=tQ=Q\). Since each row of \(A\) has absolute sum \(2\), \(\|A\|_\infty=2\). Choose \(s=2\), giving \(B=A/4\), with \(\|B\|_\infty=1/2\). Use a degree-two Taylor approximation. Direct multiplication gives \(Q^2=-2Q\), hence \(B^2=-Q/8\). Therefore $$ R_2(B)=I+B+\frac{B^2}{2}=I+\frac{Q}{4}-\frac{Q}{16} =\begin{bmatrix}\tfrac{13}{16}&\tfrac{3}{16}\\[2pt]\tfrac{3}{16}&\tfrac{13}{16}\end{bmatrix}. $$
Each row of this small-step approximation sums to \(13/16+3/16=1\). Square twice, as \(s=2\) calls for. The first square is $$ R_2(B)^2=\begin{bmatrix}\tfrac{89}{128}&\tfrac{39}{128}\\[2pt]\tfrac{39}{128}&\tfrac{89}{128}\end{bmatrix}, $$ because \(13^2+3^2=178\) and \(2(13)(3)=78\), each divided by \(16^2=256\). Squaring again gives $$ R_2(B)^4=\begin{bmatrix}\tfrac{4721}{8192}&\tfrac{3471}{8192}\\[2pt]\tfrac{3471}{8192}&\tfrac{4721}{8192}\end{bmatrix}. $$ Indeed, the diagonal numerator is \(89^2+39^2=7921+1521=9442\), and the off-diagonal numerator is \(2(89)(39)=6942\), each divided by \(128^2=16384\).
Check: Each final row sums to \((4721+3471)/8192=8192/8192=1\), and all entries are nonnegative. As an independent check of the approximation, the exact exponential has diagonal entries \((1+e^{-2})/2\) and off-diagonal entries \((1-e^{-2})/2\): the identity \(Q^2=-2Q\) reduces its power series to that form. The computed rational entries approximate those values; they are not claimed to equal them.
A Three-State Transient Calculation
Worked Example: Two Sequential Moves and an Absorbing State
Consider \(S=\{1,2,3\}\), where the chain moves from \(1\) to \(2\) at rate \(1\), from \(2\) to \(3\) at rate \(1\), and state \(3\) is absorbing. Thus $$ Q=\begin{bmatrix}-1&1&0\\0&-1&1\\0&0&0\end{bmatrix}. $$
The row sums are \(-1+1+0=0\), \(0-1+1=0\), and \(0+0+0=0\). Find an approximation at \(t=1/2\). The matrix \(A=tQ=Q/2\) has infinity norm \(1\): its first and second absolute row sums are both \(1\), and its third is \(0\). Choose \(s=1\), so \(B=A/2=Q/4\), and use degree two. Multiplication gives $$ Q^2=\begin{bmatrix}1&-2&1\\0&1&-1\\0&0&0\end{bmatrix}, \qquad R_2(B)=I+\frac{Q}{4}+\frac{Q^2}{32} =\begin{bmatrix}\tfrac{25}{32}&\tfrac{6}{32}&\tfrac{1}{32}\\[2pt]0&\tfrac{25}{32}&\tfrac{7}{32}\\[2pt]0&0&1\end{bmatrix}. $$
The first two row sums are \(25/32+6/32+1/32=1\) and \(25/32+7/32=1\); the third is \(1\). Since \(s=1\), square once. Multiplying \(R_2(B)\) by itself gives $$ R_2(B)^2= \begin{bmatrix} \tfrac{625}{1024}&\tfrac{300}{1024}&\tfrac{99}{1024}\\[2pt] 0&\tfrac{625}{1024}&\tfrac{399}{1024}\\[2pt] 0&0&1 \end{bmatrix}. $$
For example, the first diagonal entry is \((25/32)^2=625/1024\). The first-row, second-column entry is \(2(25/32)(6/32)=300/1024\). The first-row, third-column entry is \((25/32)(1/32)+(6/32)(7/32)+1/32=25/1024+42/1024+32/1024=99/1024\). In the second row, the last entry is \((25/32)(7/32)+7/32=175/1024+224/1024=399/1024\).
Check: The first row sums to \((625+300+99)/1024=1\); the second sums to \((625+399)/1024=1\); and the absorbing row sums to \(1\). Starting in state \(1\), the approximation assigns probabilities \(625/1024\), \(300/1024\), and \(99/1024\) to states \(1,2,3\), respectively. The entries are nonnegative and sum to one, as expected for a transition-probability approximation.
Choosing the Scaling and Avoiding a Common Error
In an actual computation, the scaling integer \(s\) is chosen so that \(A/2^s\) is small enough for the selected approximation method. The Taylor error bound gives one way to reason about the choice: for fixed degree \(m\), its right-hand side decreases as \(\|B\|_\infty\) is reduced. Increasing \(s\) also increases the number \(2^s\) of factors being combined, so the approximation and arithmetic costs must be balanced. Numerical libraries typically choose both a scaling amount and an approximation method using criteria designed for floating-point arithmetic.
Common mistake: scaling but not squaring enough. In the two-state example, \(s=2\) means \(B=A/4\), and recovering \(e^A\) requires four factors of \(e^B\). It is not enough to square once. For instance, using \(R_2(B)^2\) as the final answer would give diagonal entries \(89/128\) and off-diagonal entries \(39/128\), not the calculated four-factor result \(4721/8192\) and \(3471/8192\). The correction is to perform two successive squarings: first form \(R_2(B)^2\), then square that matrix to obtain \(R_2(B)^4\). A row-sum check alone would not catch this error: the one-squaring matrix also has row sums one, but it approximates the exponential at a different time scale.
Scaling and squaring is useful because it converts a potentially difficult exponential into a small-step approximation plus ordinary matrix multiplications. For a CTMC, check the generator’s row sums before computing, carry out exactly \(s\) successive squarings, and inspect the resulting rows and entries as numerical checks. Such checks help catch arithmetic or implementation errors, but an approximate matrix should not be mistaken for the exact exponential.
Check Your Understanding
For each question, distinguish the small matrix being approximated from the final matrix after squaring.
- If \(s=3\), what is the scaled matrix, and how many factors of its exponential are combined to recover \(e^A\)?
- For \(A=\begin{bmatrix}-2&2\\0&0\end{bmatrix}\), compute \(\|A\|_\infty\). What is \(\|A/4\|_\infty\)?
- Why does the Taylor error bound become smaller when \(\|B\|_\infty\) is reduced while the degree \(m\) is held fixed?
- In the sequential three-state example, verify the first-row sum of the final approximation using its three displayed numerators.
- Why does a row-sum check fail to detect using only one squaring when \(s=2\)?