Numerical methods for quasi-stationary distributions

TL;DR

This paper addresses the computation of quasi-stationary distributions (QSDs) for Markov processes with absorbing states, addressing both discrete and continuous state spaces and multiple absorbing states. It revisits two numerical approaches: (i) a generalized iterative scheme that solves the defining non-linear equation for the QSD, and (ii) a Monte Carlo resetting method, including a novel single-trajectory variant that uses the trajectory’s own history for resets. The authors provide detailed implementation guidance, compare accuracy and efficiency across a diverse set of models, and offer practical recommendations for parameter choices. Overall, the iterative method is typically faster and more accurate for simple boundaries, while the Monte Carlo resetting approach excels for problems with complex boundaries or large tail probabilities; together they offer a versatile toolkit for QSD computation and analysis.

Abstract

In stochastic processes with absorbing states, the quasi-stationary distribution provides valuable insights into the long-term behaviour prior to absorption. In this work, we revisit two well-established numerical methods for its computation. The first is an iterative algorithm for solving the non-linear equation that defines the quasi-stationary distribution. We generalise this technique to accommodate general Markov stochastic processes, either with discrete or continuous state space, and with multiple absorbing states. The second is a Monte Carlo method with resetting, for which we propose a novel single-trajectory approach that uses the trajectory's own history to perform resets after absorption. In addition to these methodological contributions, we provide a detailed analysis of implementation aspects for both methods. We also compare their accuracy and efficiency across a range of examples. The results indicate that the iterative algorithm is generally the preferred choice for problems with simple boundaries, while the Monte Carlo approach is more suitable for problems with complex boundaries, where the implementation of the iterative algorithm is a challenging task.

Figures (12)

• Figure 1: Quasi-stationary distribution and mean time to extinction for the discrete processes in Sec. \ref{['sec:ExamplesDiscrete']}: (a) the linear branching process, (b) the immigration-death process, (c) the branching-annihilation-decay (B-A-D) process with $V = 250$, and (d) the biased voter model with $N = 100$. Main panels: Quasi-stationary distribution for different values of the parameter $R$ (see text). Black crosses come from the iterative algorithm, while coloured symbols come from the single-trajectory Monte Carlo method. Additionally, in panels (a)-(b), solid lines represent the analytical solution. For better visualisation, in cases (c) for $R = 1.0, 1.1, 1.2$ and (d), the numerical distribution is plotted only for odd values of $n$. In both numerical methods, the initial condition has been set to: (a) $n_0 = 10$, (b) $n_0 = R$ (metastable state) (c) $n_0 = 10$ if $R \leq 1$ and $n_0 = 1 + V(1 - 1/R)$ (metastable state) if $R > 1$, and (d) $n_0 = 10, 50, 90$ for $h = 0.9, 1.0, 1.1$, respectively. Additionally, in cases (a)-(c) we have set an upper bound for the number of particles at $N = 150$ in the iterative algorithm. Insets: Mean absorption time, $\tau$, for different choices of $R$. Grey bars come from the iterative algorithm, filled coloured bars come from the Monte Carlo method and, in cases (a) and (b), bars with diagonal lines represent the analytical solution.
• Figure 2: Quasi-stationary density function and mean absorption time for the continuous processes in Sec. \ref{['sec:ExamplesContinuous']}: (a) the biased random walk with $D = 0.4$ and $L = 1$, (b) the Ornstein--Uhlenbeck process with $D = 0.3$, (c) the linear branching process with $f = 5$, and (d) the continuous biased voter model with $N = 100$. Main panels: Quasi-stationary density function for different values of a characteristic parameter of the model (see text). Dotted coloured lines correspond to the iterative algorithm, while dashed coloured lines correspond to the Monte Carlo method. Moreover, in panels (a)-(c), solid grey lines represent the analytical solution. In most cases, the different lines are visually indistinguishable. In the numerical methods, we have used as initial condition: (a) $x_0 = 0.1$ for $f = -1$, and $x_0 = 0.9$ for $f = 0.5, 2.0$, (b)-(c) $x_0 = 0.1$, and (d) $x_0 = 0.1, 0.5, 0.9$ for $h = 0.9, 1.0, 1.1$, respectively. Additionally, in cases (b)-(c) we have introduced a reflecting upper boundary in the iterative algorithm at $L = 3$ and $L = 2$, respectively. Insets: Mean absorption time, $\tau$, for different choices of the characteristic parameter. Grey bars come from the iterative algorithm, filled coloured bars come from the Monte Carlo method and, in cases (a)-(c), bars with diagonal lines represent the analytical solution.
• Figure 3: Results for the two-dimensional processes in Sec. \ref{['sec:Examples2D']}. (a),(c) Quasi-stationary distribution obtained with the Monte Carlo method for: (a) the SIRS model (with $N = 20$ and $R = 1.4$) for $K = 10^9$ steps and initial condition $(m_0, n_0) = \left(\frac{N}{R}, N \frac{R - 1}{2R} \right)$ (metastable state), and (c) the two-dimensional diffusion process with circular symmetry (with $D = 0.5$, $r_0 = 0.3$ and $r_1 = 1.2$) for $K = 10^{11}$ steps, discretisation parameters $\Delta x = \Delta y = 5 \cdot 10^{-3}$ and $\Delta t = 10^{-4}$, and initial condition $(x_0, y_0)$ randomly sampled within the annular domain. (b) Marginal quasi-stationary distribution of the number of infected individuals, $Q(\cdot,n) = \sum_{m=1}^{N-1} Q(m,n)$, for the SIRS model (with $N = 20$ and varying $R$) computed using the iterative algorithm with a relaxation factor $s = 0.1$ (black crosses) and the Monte Carlo method (coloured symbols). The initial condition has been set to $(m_0, n_0) = (N-10, 10)$ for $R = 0.6, 1.0$ and identically as in (a) for $R = 1.4$. Inset: Mean time to extinction, defined as the average time the system takes to reach any absorbing state of the form $(m, 0)$ (i.e., with no infected individuals), as a function of $R$. (d) Radial distribution for the two-dimensional diffusion process with circular symmetry (with parameters as in (c)), obtained in the Monte Carlo simulation with $\Delta r = 10^{-3}$ (pink) and compared to the analytical solution (black).
• Figure 4: Comparison of efficiency and accuracy of the numerical methods. Computation time (measured in seconds) as a function of the total error in the numerical estimate of the quasi-stationary distribution. Blue and pink symbols come from the iterative algorithm and the Monte Carlo method, respectively. The initial conditions have been set as in Sec. \ref{['sec:Examples']}. (a) Data for the discrete linear branching process (empty symbols) and the immigration-death process (filled symbols) described in Sec. \ref{['sec:ExamplesDiscrete']}. In both cases we have set $R = 0.6$, used Eq. \ref{['eq:TotalErrorDiscrete']} to compute the errors, and set an upper bound at $N = 150$ in the iterative algorithm. (b) Results for the continuous linear branching process (empty symbols) for $D = 0.3$ and $f = 5$, and the Ornstein--Uhlenbeck process (filled symbols) for $D = 0.3$ and $f = 3.5$ (see Sec. \ref{['sec:ExamplesContinuous']}). Errors have been calculated using Eq. \ref{['eq:TotalErrorContinuous']}. In the iterative algorithm, we have introduced a reflecting boundary condition at $L = 2$ for the linear branching process, and at $L = 3$ for the Ornstein--Uhlenbeck process. All simulations in this paper have been run on an AMD Epyc Rome processor.
• Figure 5: Sampling of small probabilities. Quasi-stationary distribution and mean time to extinction for the immigration-death process introduced in Eq. \ref{['eq:Rates_Mix']} for large values of the reproduction ratio $R$. Main panel: Quasi-stationary distribution, obtained with the iterative algorithm (black crosses), the Monte Carlo method (coloured symbols), and the analytical solution (solid lines). For better visualisation, the numerical distribution is plotted only for odd values of $n$. In the numerical methods, we have set the initial condition at $n_0 = R$. Moreover, we have set $N = 500$ as the upper bound for the population size in the iterative algorithm. Inset: Mean time to extinction, $\tau$, as a function of $R$. Results from the iterative algorithm are shown as grey bars, results from the Monte Carlo method as filled coloured bars, and the analytical solution as bars with diagonal lines. As noted in the text, results for the mean time to extinction can only be obtained from the Monte Carlo method within manageable computing time for $R=20$.