Why Numerical Stability Matters
A product-integral computation applies a sequence of transition updates in time order. Even when every update is valid, a computer represents most real numbers approximately. Small arithmetic errors can therefore accumulate over many multiplications. The practical question is not whether finite-precision error exists, but whether it is small enough to leave the reported transition probabilities trustworthy.
For a finite-state model, use the row convention: row \(i\) describes transitions from origin state \(i\), and a row vector of state probabilities is multiplied on the right. If \(A_1,\ldots,A_K\) are the successive update matrices, the cumulative transition matrix is \(P=A_1A_2\cdots A_K\). The order is chronological and must not be changed. When the update matrices are valid stochastic matrices, their entries are nonnegative and every row sums to one. The exact product has those same properties.
This tutorial focuses on the numerical computation of that product, not on the statistical construction of the transition updates. The earlier material on discrete event-time computation gives the product-integral setting. Here, assume the updates have been correctly formed and are stochastic; numerical checks can help identify when that assumption is violated in the computation.
What Rounding Can Change
In floating-point arithmetic, an operation is typically modeled as the exact result times \(1+\delta\), where \(|\delta|\leq u\) and \(u\) is the unit roundoff, provided underflow and overflow are not involved. The exact value of \(u\) depends on the numeric format. A matrix product contains many additions and multiplications, and a long product repeats those operations. Thus, retaining full working precision at every update matters.
Rounding every intermediate probability to a few decimal places is not the same as ordinary floating-point computation. It deliberately discards substantially more information. The following simple model shows why this can create a large error even when each individual step seems mild.
Worked Example: Repeated Rounding in a 100-Step Follow-Up
Suppose there are two states: state 0 and an absorbing state 1. At each of 100 successive updates, a probability of 0.01 moves from state 0 to state 1. The update matrix at each step is
Starting in state 0, the probability of remaining in state 0 after \(K\) updates is \(q_K=0.99^K\). The full transition matrix is
Using unrounded intermediate values, \(q_{100}=0.99^{100}\approx0.366032\). Thus the probability of being in the absorbing state is \(1-0.366032\approx0.633968\).
Now instead round \(q_K\) to two decimal places after every step. Starting at \(q_0=1.00\), the update is \(q_{K+1}=\operatorname{round}(0.99q_K,2)\). When the rounded value reaches \(0.50\), the next unrounded value is \(0.50(0.99)=0.495\). Under common decimal half-up or half-to-even rounding, this rounds to \(0.50\) again. Some binary implementations may settle at \(0.49\) instead. Continued repeated rounding can therefore keep this probability near \(0.50\), with the exact rounded fixed value depending on tie-breaking and implementation, although its correct value after 100 updates is about \(0.366\). The discrepancy is not a small final-digit effect: the estimated absorbing-state probability would be near \(0.50\), rather than about \(0.634\).
A stable calculation retains the working-precision value at every step and rounds only when presenting the final result. For this scalar product, an alternative calculation is \(q_{100}=\exp\{100\log(0.99)\}\approx0.366032\). The logarithmic calculation is useful for this product of positive scalars; it does not replace ordered matrix multiplication for general transition updates.
Stochasticity as a Check on Matrix Products
A stochastic matrix has nonnegative entries and rows summing to one. These properties are preserved by exact multiplication. To see this, let \(A\) and \(B\) be stochastic matrices of the same dimension. Their product has entries \(C_{ij}=\sum_{\ell} A_{i\ell}B_{\ell j}\), so each entry is nonnegative. Also,
The argument extends repeatedly to a product of any number of stochastic updates. In a numerical calculation, computed row sums may differ slightly from one; that is a diagnostic signal, not by itself proof of a substantive error. Check the size and pattern of the discrepancy against the number of operations, the working precision, and the scale of the computation. Also inspect whether any entries are meaningfully negative or greater than one.
Worked Example: Checking the Order and Row Sums of Two Updates
Consider two hypothetical stochastic updates,
Each row sums to one. Applying \(A\) first and then \(B\) gives \(AB\):
The row sums are \(0.78+0.22=1\) and \(0.30+0.70=1\). Reversing the order gives
These rows also sum to one, but \(BA\ne AB\). A row-sum check can confirm a basic probability property; it cannot confirm that updates were applied in the right chronological order. Both checks are needed.
How Large Can One Multiplication’s Error Be?
A dot product in a matrix multiplication is a sum of products. Under the standard floating-point model, a computed dot product of length \(m\), with no underflow or overflow, has an absolute rounding error bounded by approximately \(\gamma_m\) times the sum of the absolute products, where
This is a bound for one dot product; it is not a promise that every calculation attains the bound. For a single product of exact stochastic matrices \(A\) and \(B\), all products \(A_{i\ell}B_{\ell j}\) are nonnegative. Applying the dot-product bound to each entry and summing across row \(i\) gives an approximate upper bound \(\gamma_m\) on the absolute row-sum discrepancy from that multiplication:
Here \(C=AB\) is the exact product and \(\widehat C\) is the computed product. This justification assumes the input matrices are exact and stochastic and that the stated arithmetic model applies. In a long calculation, the inputs to later multiplications are themselves computed values, so their errors also matter. The bound should therefore guide expectations, not be used as a universal tolerance for an entire product.
Worked Example: A Row-Sum Residual After a Transition Update
Suppose a hypothetical three-state update is
and the probability row vector before the update is \(\pi=(0.70,0.20,0.10)\). The update gives
The entries sum to \(0.595+0.276+0.129=1.000\). If the three entries are instead rounded independently to two decimal places, the displayed vector becomes \((0.60,0.28,0.13)\), with sum \(1.01\). Its row-sum residual is \(1.01-1=0.01\). That residual is explained by coarse display rounding here; it is not evidence that the exact update failed to preserve probability. The calculation should retain the unrounded values, and the report should state a suitable number of digits.
A Stable Workflow for Long Products
A reliable implementation separates the statistical updates from the arithmetic used to accumulate them. The following sequence makes common computational failures easier to locate.
Before multiplication, check that entries are finite, no entry is materially negative, and each row sums to one within a tolerance appropriate to the precision used to construct the update.
Starting from the identity matrix, update the cumulative product on the right: \(P_k=P_{k-1}A_k\). Do not reorder updates merely to change the computation.
Do not round probabilities or matrix entries after each event time. Use the available standard precision, or higher precision when the required accuracy or diagnostics make it necessary.
Track row sums and the smallest and largest entries. A steadily growing row-sum defect, substantial negative entry, or entry above one calls for investigation rather than silent correction.
Recompute a small case independently, use higher precision, or compare with a direct calculation when feasible. If the discrepancy persists, check the input updates and multiplication code, not only the displayed rounding.
For dot products, pairwise summation or compensated summation can reduce the loss of low-order digits when many terms are added. These approaches improve arithmetic; they do not repair incorrect updates or the wrong product order. Likewise, multiplying in higher precision cannot make an invalid input matrix valid.
Diagnostics, Alternatives, and Limitations
A row-sum check is useful because probability mass should be preserved by a stochastic update. Check each row of a transition matrix, not only the row vector for one initial distribution: an error in an unoccupied origin state may otherwise go undetected. In addition, inspect individual entries. Tiny negative values near the scale of roundoff can arise in some computations, but a substantial negative probability or an entry substantially greater than one is not explained by ordinary rounding.
Do not automatically renormalize each computed row to sum to one. Renormalization changes the numerical result and can conceal an omitted transition, an incorrectly formed diagonal entry, or an update applied in the wrong way. First locate the source of a material discrepancy and recompute if needed. A documented correction may be defensible when the deviations are demonstrably only tiny rounding residuals, but it should not replace diagnosis.
| Approach | Useful for | Important limitation |
|---|---|---|
| Sequential full-precision matrix updates | Computing an ordered product of event-time transition updates | Roundoff can accumulate; monitor stochasticity throughout the calculation. |
| Higher precision or improved summation | Reducing arithmetic error when ordinary-precision diagnostics are inadequate | Cannot correct invalid updates, incorrect ordering, or a mistaken model. |
| Matrix exponential for piecewise-constant intensities | Computing transitions under a constant-intensity segment | It is a different computation and does not generally replace the product of event-time updates. |
Diagnostics have limits. A matrix with rows summing to one may still be wrong because its updates were ordered incorrectly or its entries were computed from incorrect inputs. Conversely, a very small row-sum residual may be an ordinary rounding effect. Interpretation depends on the update construction, matrix dimension, number of operations, precision, and size of the residual. For a long follow-up calculation, save the input updates and intermediate diagnostics so the first point of divergence can be identified.
Common Mistake
A related mistake is rounding the running probabilities for a cleaner-looking table. As the 100-update example shows, repeated rounding can produce a stable value far from the correct probability. Round for presentation, not after every computational step.
Check Your Understanding
Try each question, then compare your reasoning with the solutions.
- Two stochastic matrices \(A\) and \(B\) are applied in that chronological order under the row convention. Which product gives the cumulative update, and does checking row sums establish that the order is correct?
- For a probability that remains in a state with probability \(0.99\) at each of 100 updates, give its unrounded value to four significant figures. Why can rounding to two decimals after every update be misleading?
- A computed transition-matrix row sums to \(1.0000002\), and its entries are all between zero and one. Is this automatically a substantive error? What would you examine?
- Why is automatic row renormalization not a sufficient response to a row sum of \(0.97\)?
Solutions. (1) The cumulative update is \(AB\), because row distributions multiply on the right and the first update is applied first. The row sums do not establish the order: \(BA\) can also be stochastic while differing from \(AB\). (2) The value is \(0.99^{100}\approx0.3660\). Rounding at every step discards information repeatedly and can keep the computed value near \(0.50\), far from the correct result. (3) Not automatically. A residual of \(2\times10^{-7}\) may be consistent with the precision and number of operations, but assess it in context, check earlier updates and intermediate residuals, and compare with a higher-precision or independent computation if needed. (4) A row sum of \(0.97\) may indicate a missing transition, a malformed update, or an implementation error. Renormalization hides the discrepancy and changes the entries without identifying its cause; inspect and correct the computation first.