Tutorials › Stochastic Processes › Comparing Uniformization and Direct Matrix Exponentials

Kolmogorov Equations and Transient Behavior · Tutorial 512 of 994

Comparing Uniformization and Direct Matrix Exponentials

Learn to compare uniformization and direct matrix-exponential calculations by checking their accuracy, computational work, and behaviour on sparse generators.

Advanced 12 min read

What You'll Learn

  • Relate the uniformization series to the direct matrix exponential.
  • Distinguish truncation error from numerical error in evaluating the matrix exponential.
  • Compare the two methods on a two-state chain and a sequential-transition model.
  • Track sparse matrix-vector work using the nonzero entries of the embedded transition matrix.
  • Choose a method with attention to accuracy, sparsity, and the requested output.

Two Ways to Compute Transient Probabilities

For a finite-state continuous-time Markov chain, the transient transition matrix at time \(t\) is \(P(t)=e^{tQ}\), where \(Q\) is the generator. A direct matrix-exponential method evaluates this exponential, often using scaling and squaring with a rational approximation. Uniformization instead represents the same transition matrix as a Poisson-weighted sum of powers of a stochastic matrix. The two methods target the same answer, but their computational work and sources of error differ.

We compare these methods for a process \((X_t)_{t\ge0}\) on a finite state space \(S\), defined on a probability space \((\Omega,\mathcal F,\mathbb P)\). The generator \(Q=(q_{ij})\) has nonnegative off-diagonal entries and rows summing to zero. A distribution is a row vector \(\mu\), so its distribution at time \(t\) is \(\mu P(t)\). Choose a uniformization rate \(\nu\ge \max_i(-q_{ii})\), with \(\nu>0\), and define \(R=I+Q/\nu\). Then \(R\) is row-stochastic: its entries are nonnegative and each row sums to one. Uniformization expresses the transition matrix as

$$ P(t)=e^{tQ}=e^{-\nu t}\sum_{k=0}^{\infty}\frac{(\nu t)^k}{k!}R^k. $$

Thus, uniformization amounts to averaging discrete-time transitions \(R^k\) over a Poisson random number of steps with mean \(\nu t\). A direct matrix-exponential routine does not generally form these powers or the Poisson weights. It approximates \(e^{tQ}\) using matrix operations, typically with scaling and squaring. The identity above is justified below; the uniformization formula and Poisson-tail bounds are also standard results from the earlier material on uniformization.

Comparison principle: Both methods compute the same \(P(t)\). A finite uniformization sum has a controllable truncation error, while a direct matrix-exponential routine has approximation and rounding error determined by its implementation and tolerances. Work estimates depend on the matrix structure and algorithm; no method is uniformly faster for every input.

Why the Two Formulas Agree

The central algebraic fact is that \(Q=\nu(R-I)\). Since \(R\) commutes with \(I\), the exponential of this sum factors. This gives the uniformization series directly from the power series for the matrix exponential.

Proposition: If \(Q\) is a finite-state generator and \(\nu\ge\max_i(-q_{ii})\) with \(\nu>0\), then \(R=I+Q/\nu\) is row-stochastic and the uniformization identity holds for every \(t\ge0\).

Proof. For \(i\ne j\), \(R_{ij}=q_{ij}/\nu\ge0\), and \(R_{ii}=1+q_{ii}/\nu\ge0\) by the condition on \(\nu\). For each row \(i\), \(\sum_j R_{ij}=1+\nu^{-1}\sum_jq_{ij}=1\), so \(R\) is row-stochastic. Also, \(Q=\nu(R-I)\). The matrices \(R\) and \(I\) commute, so the exponential power series gives \(e^{tQ}=e^{\nu t(R-I)}=e^{-\nu t}e^{\nu tR}\). Expanding the remaining matrix exponential yields \(e^{-\nu t}\sum_{k\ge0}(\nu t)^kR^k/k!\). The series converges because it is the matrix exponential series in finite dimension. \(\square\)

For computation, stop the sum at \(K\), obtaining \(P_K(t)=e^{-\nu t}\sum_{k=0}^K(\nu t)^kR^k/k!\). This is generally substochastic: its row sums are the probability that a Poisson random variable with mean \(\nu t\) is at most \(K\). The omitted mass is the Poisson upper tail, not a defect in the generator or in the direct matrix exponential.

Proposition (truncation error): For any probability row vector \(\mu\), let \(N\) have the Poisson distribution with mean \(\nu t\). Then \(\|\mu P(t)-\mu P_K(t)\|_1=\mathbb P(N>K)\). In particular, the \(L^1\) error is at most the chosen Poisson tail.

Proof. The full and truncated vectors differ by \(e^{-\nu t}\sum_{k>K}(\nu t)^k\mu R^k/k!\). Each \(\mu R^k\) is a probability row vector, since \(R\) is stochastic; hence every entry of the difference is nonnegative. Its \(L^1\) norm is therefore its sum of entries. Since \(\sum_j(\mu R^k)_j=1\), that sum is \(e^{-\nu t}\sum_{k>K}(\nu t)^k/k!\), which is exactly \(\mathbb P(N>K)\). \(\square\)

Worked Comparisons

Worked Example: A Two-State Chain

Consider \(S=\{1,2\}\) and \(Q=\begin{bmatrix}-2&2\\1&-1\end{bmatrix}\). The row sums are \(-2+2=0\) and \(1-1=0\), as required. Take \(\nu=2\), so \(R=I+Q/2=\begin{bmatrix}0&1\\ \tfrac12&\tfrac12\end{bmatrix}\). Its row sums are \(0+1=1\) and \(\tfrac12+\tfrac12=1\); it is a valid transition matrix.

Starting in state \(1\), the direct computation solves the forward equations for \(p(t)=(p_1(t),p_2(t))\). Since \(p_1(t)+p_2(t)=1\) and \(p_1'(t)=-2p_1(t)+p_2(t)=1-3p_1(t)\), with \(p_1(0)=1\), solving gives \(p_1(t)=\tfrac13+\tfrac23e^{-3t}\). Therefore \(p_2(t)=1-p_1(t)=\tfrac23(1-e^{-3t})\). At \(t=1\), the direct answer is \((\tfrac13+\tfrac23e^{-3},\tfrac23(1-e^{-3}))\).

Uniformization gives the same answer from \(e^{-2}\sum_{k\ge0}2^kR^k/k!\). To see why it agrees without evaluating every power, substitute \(Q=2(R-I)\) into the matrix exponential: \(e^Q=e^{-2}e^{2R}\). The direct solution above is \(e^Q\) applied to the row vector \((1,0)\), so the two calculations agree entry by entry. This example has only three nonzero entries in \(R\). Either method is small here; for a larger sparse matrix, forming a dense approximation may cost more than sparse matrix-vector products.

Worked Example: Truncation on a Sequential Chain

Take \(S=\{1,2,3\}\), with transitions \(1\to2\) and \(2\to3\) each at rate \(1\), and state \(3\) absorbing. Thus \(Q=\begin{bmatrix}-1&1&0\\0&-1&1\\0&0&0\end{bmatrix}\). The row sums are \(-1+1=0\), \(0-1+1=0\), and \(0\). Choose \(\nu=1\); then \(R=\begin{bmatrix}0&1&0\\0&0&1\\0&0&1\end{bmatrix}\), whose row sums are all \(1\).

For a direct calculation, let \(p(t)\) be the row of transition probabilities starting in state \(1\). The forward equations are \(p_1'=-p_1\), \(p_2'=p_1-p_2\), and \(p_3'=p_2\), with \(p(0)=(1,0,0)\). Solving in order gives \(p_1(t)=e^{-t}\), \(p_2(t)=te^{-t}\), and \(p_3(t)=1-(1+t)e^{-t}\); the last expression follows from \(p_1+p_2+p_3=1\). At \(t=1\), the exact row is \((e^{-1},e^{-1},1-2e^{-1})\).

For uniformization truncated at \(K=2\), the vectors \(v_k=(1,0,0)R^k\) are \(v_0=(1,0,0)\), \(v_1=(0,1,0)\), and \(v_2=(0,0,1)\). The Poisson weights at mean \(1\), for \(k=0,1,2\), are \(e^{-1}\), \(e^{-1}\), and \(e^{-1}/2\). Hence the truncated row is \((e^{-1},e^{-1},e^{-1}/2)\). Its difference from the exact row is \((0,0,1-\tfrac52e^{-1})\), which is exactly \(\mathbb P(N>2)=1-e^{-1}(1+1+\tfrac12)\) for \(N\) Poisson with mean \(1\). Thus the first two state probabilities are already exact in this truncation, while the third is short by the omitted Poisson tail. This illustrates that a global error bound can be conservative for individual entries.

Worked Example: Sparse Work and a Direct Check

Use \(S=\{1,2,3,4\}\), where \(1\to2\) and \(2\to3\) each have rate \(1\), and states \(3\) and \(4\) are absorbing. State \(4\) is an isolated state included to show the sparsity pattern. The generator is \(Q=\begin{bmatrix}-1&1&0&0\\0&-1&1&0\\0&0&0&0\\0&0&0&0\end{bmatrix}\); each row sum is \(0\). With \(\nu=2\), \(R=\begin{bmatrix}\tfrac12&\tfrac12&0&0\\0&\tfrac12&\tfrac12&0\\0&0&1&0\\0&0&0&1\end{bmatrix}\). Its row sums are \(\tfrac12+\tfrac12=1\), \(0+\tfrac12+\tfrac12=1\), \(1\), and \(1\).

There are exactly six nonzero entries in \(R\): two in row 1, two in row 2, one in row 3, and one in row 4. Starting from state 1, successive sparse row-vector multiplications give \(v_0=(1,0,0,0)\), \(v_1=(\tfrac12,\tfrac12,0,0)\), \(v_2=(\tfrac14,\tfrac12,\tfrac14,0)\), \(v_3=(\tfrac18,\tfrac38,\tfrac12,0)\), and \(v_4=(\tfrac1{16},\tfrac14,\tfrac{11}{16},0)\). Each vector sums to \(1\), checking the multiplication results. These four multiplications involve at most \(4\cdot6=24\) nonzero matrix-entry contributions. The count is an upper bound: some entries need not contribute for a particular input vector.

At \(t=1\), the Poisson mean is \(\nu t=2\). The weights for \(k=0,\ldots,4\) are \(e^{-2}\), \(2e^{-2}\), \(2e^{-2}\), \(\tfrac43e^{-2}\), and \(\tfrac23e^{-2}\). Multiplying and adding the vectors gives the truncated result \(e^{-2}(\tfrac{65}{24},\tfrac83,\tfrac{13}{8},0)\). Its entries sum to \(7e^{-2}\), equal to the probability that a Poisson random variable with mean \(2\) is at most \(4\). The exact direct result at \(t=1\), from the preceding sequential-chain equations, is \((e^{-1},e^{-1},1-2e^{-1},0)\). The \(L^1\) truncation error is therefore \(1-7e^{-2}\), the omitted Poisson tail. The direct result has row sum \(e^{-1}+e^{-1}+1-2e^{-1}=1\); the truncated result is not normalized because its tail has been omitted.

Accuracy, Runtime, and Method Choice

These examples separate three issues that are easy to conflate:

  • Model validity. The generator must have nonnegative off-diagonal entries and zero row sums. For uniformization, the chosen \(\nu\) must also make \(R\) stochastic.
  • Approximation error. Truncating uniformization at \(K\) omits exactly a Poisson tail in \(L^1\) for any starting distribution. A direct matrix-exponential routine has its own approximation and rounding errors, controlled through its algorithm and tolerance rather than a Poisson truncation index.
  • Computational work. Uniformization can exploit sparsity by repeatedly multiplying a row vector by \(R\). The cost per multiplication is proportional to the number of nonzero entries when sparse storage and multiplication are used. Direct methods such as scaling and squaring can be effective, especially for moderate dense matrices or when the full exponential is required, but intermediate matrix products can create additional nonzero entries.

These are structural comparisons, not universal runtime rankings. Actual elapsed time depends on matrix size, sparsity, hardware, implementation, requested tolerance, and whether the task needs one initial distribution or the entire transition matrix. Uniformization is particularly attractive when the generator is sparse, only a few starting distributions are needed, and a Poisson-tail cutoff is manageable. A direct exponential can be preferable when a robust matrix-exponential library is available, the matrix is small or dense, or the full \(P(t)\) is needed. A very large value of \(\nu t\) can require many uniformization terms, so sparsity alone does not guarantee a speed advantage.

Common mistake. A finite uniformization sum is not automatically a probability distribution. For example, in the three-state calculation truncated at \(K=2\), the entries sum to \(e^{-1}+e^{-1}+\tfrac12e^{-1}=\tfrac52e^{-1}\), not \(1\). Renormalizing these entries changes the approximation and invalidates the exact Poisson-tail error statement. The correction is to report the omitted tail \(1-\tfrac52e^{-1}\), increase \(K\) until it meets the required tolerance, or use a separately justified approximation method. Also, counting all matrix positions as computational work can badly misrepresent sparse calculations: count nonzero contributions instead, and do not infer a machine-time comparison from that count alone.

Check Your Understanding

Use the definitions and calculations above to check each conclusion.

  1. For a generator \(Q\), what condition on \(\nu\) ensures that \(R=I+Q/\nu\) is row-stochastic?
  2. For Poisson mean \(3\) and truncation at \(K=2\), write the exact \(L^1\) truncation error without evaluating it numerically.
  3. Why does a truncated uniformization sum generally have row sum less than \(1\)?
  4. In the four-state example, how many nonzero entries does \(R\) have, and what upper bound follows for four successive sparse row-vector multiplications?
  5. Give one feature of a problem that may favour uniformization and one that may favour a direct matrix exponential.