Data-driven computation methods for quasi-stationary distribution and sensitivity analysis

We first give definition of the QSD and the exponential killing rate $\lambda$ of a Markov process with an absorbing state. Let $X=(X_{t}:t\geq 0)$ be a continuous-time Markov process taking values in a measurable space $(\mathcal{X},\mathcal{B}(\mathcal{X}))$. Let $P^{t}(x,\cdot)$ be the transition kernel of $X$ such that $P^{t}(x,A)=\mathbb{P}[X_{t}\in A\,|\,X_{0}=x]$ for all $A\in\mathcal{B}$. Now assume that there exists an absorbing set $\partial\mathcal{X}\subset\mathcal{X}$. The complement $\mathcal{X}^{a}:=\mathcal{X}\backslash\partial\mathcal{X}$ is the set of allowed states.

The process $X_{t}$ is killed when it hits the absorbing set, implying that $X_{t}\in\partial\mathcal{X}$ for all $t>\tau$, where $\tau=\inf\{t>0:X_{t}\in\partial\mathcal{X}\}$ is the hitting time of set $\partial\mathcal{X}$. Throughout this paper, we assume that the process is almost surely killed in finite time, i.e. $\mathbb{P}[\tau<\infty]=1$.

For the sake of simplicity let $\mathbb{P}_{x}$ (resp. $\mathbb{P}_{\mu}$) be the probability conditioning on the initial condition $x\in\mathcal{X}$ (resp. the initial distribution $\mu$).

In the analysis of QSD, we are particularly interested in a parameter $\lambda$, called the killing rate of the Markov process. If the distribution of the killing time $\mathbb{P}_{x}(\tau>t)$ has an exponential tail, then $\lambda$ is the rate of this exponential tail. The following theorem shows that the killing time is exponentially distributed when the process starts from a QSD[5].

Throughout this paper, we assume that $X$ admits a QSD denoted by $\mu$ with a strictly positive killing rate $\lambda$.

Consider a stochastic differential equation

where $X_{t}\in\mathbb{R}^{d}$ and $X_{t}$ is killed when it hits the absorbing set $\partial\mathcal{X}\subset\mathbb{R}^{n}$; $f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is a continuous vector field; $\sigma$ is an $d\times d$ matrix-valued function; and $\text{d}W_{t}$ is the white noise in $\mathbb{R}^{d}$. The following well known theorem shows the existence and the uniqueness of the solution of equation (2.3).

There are quite a few recent results about the existence and convergence of QSD. Since the theme of this paper is numerical computations, in this paper we directly assume that $X_{t}$ admits a unique QSD $\mu$ on set $\mathcal{X}^{a}$ that is also the quasi-limit distribution. The detailed conditions are referred in [18, 13, 20, 9, 21].

Let $u$ be the probability density function of $\mu$. We refer [2] for the fact that $u$ satisfies

where $D=\sigma^{T}\sigma$, and $\lambda$ is the killing rate.

It remains to review the simulation algorithm for QSDs. In order to compute the QSD, one needs to numerically simulate a long trajectory of $X$. Once $X_{t}$ hits the absorbing state, a new initial value is sampled from the empirical QSD. The re-sampling step can be done in two different ways. We can either use many independent trajectories that form an empirical distribution [3] or re-sample from the history of a long trajectory [19]. In this paper we use the latter approach.

Let $\hat{X}^{\delta}=\{\hat{X}^{\delta}_{n},\ n\in\mathbb{Z}_{+}\}$ be a long numerical trajectory of the time-$\delta$ sample chain of $X_{t}$, where $\hat{X}^{\delta}_{n}$ is an numerical approximation of $X_{n\delta}$, then the Euler-Maruyama numerical scheme is given by

where $\hat{X}^{\delta}_{n}=X_{0}$, $W_{(n+1)\delta}-W_{n\delta}\sim\mathcal{N}(0,\delta\mathrm{Id}_{d}),\ n\in\mathbb{Z}_{+}$ is a $d$-dimensional normal random variable.

Another widely used numerical scheme is called the Milstein scheme, which reads

where $I$ is a $d\times d$ matrix with its $(i,j)$-th component being the double It$\hat{o}$ integral

and $L\in\mathbb{R}^{d}$ is a vector of left operators with $i$-th component

For simplicity, we introduce the algorithm for $n=2$, specifically, we solve $u$ in equation (2.4) numerically on a 2D domain $D=[a_{0},b_{0}]\times[a_{1},b_{1}]$. Firstly, we construct an $N\times M$ grid on $D$ with grid size $h=\frac{b_{0}-a_{0}}{N}=\frac{b_{1}-a_{1}}{M}$. Each small box in the mesh is denoted by $O_{i,j}=[a_{0}+(i-1)h,a_{0}+ih]\times[a_{1}+(j-1)h,a_{1}+jh]$. Let $\mathbf{u}=\{u_{i,j}\}^{i=N,j=M}_{i=1,j=1}$ be the numerical solution on $D$ that we are interested in, then $\mathbf{u}$ can be considered as a vector in $\mathbb{R}^{N\times M}$. Each element $u_{i,j}$ approximates the density function $u$ at the center of each $O_{i,j}$, with coordinate $(ih+a_{0}-h/2,jh+a_{1}-h/2)$. Generally speaking, we count the number of $\hat{X}^{\delta}$ falling into each box and set the normalized value as the approximate probability density at $O_{i,j}$. As we are interested in the QSD, the main difference from the general Monte Carlo is that the Markov process will be regenerated as the way of its empirical distribution uniformly once it is killed. The details of simulation is shown in Algorithm 1 as following.

Sometimes the Euler-Maruyama method underestimates the probability that $X$ moves to the absorbing set, especially when $X_{t}$ is close to $\partial\mathcal{X}$. This problem can be fixed by introducing the Brownian bridge correction. We refer to [4] for details. For a sample falling into a small box which are closed to $\partial\mathcal{X}$, the probability of them falling into the trap $\partial\mathcal{X}$ is relatively high. In fact, this probability is exponentially distributed and the rate is related to the distance from $\partial\mathcal{X}$. Let $B^{T}_{B}=W_{t}-\frac{t}{T}W_{T}$ be the Brownian Bridge on the interval $[0,T]$. In the 1D case, the law of the infimum and the supremum of the Brownian Bridge can be computed as follows: for every $z\geq\max(x,y)$

where $x=\hat{X}^{\delta}_{n}\in\mathcal{X}^{a},y=\hat{X}^{\delta}_{n+1}\in\mathcal{X}^{a}$, and $\phi=\sigma(\hat{X}^{\delta}_{n})$ is the strength coefficient of Brownian Bridge. This means that at each step $n$, if $\hat{X}^{\delta}_{n+1}\in\mathcal{X}^{a}$, one can compute, with the help of the above properties, a Bernoulli random variable $G$ with the parameter

The coupling method is used for the sensitivity analysis of QSDs.

The definition of coupling can be extended to any two random variables that take value in the same state space. Now consider two Markov processes $X=(X_{t}:t\geq 0)$ and $Y=(Y_{t}:t\geq 0)$ with the same transition kernel $P$. A coupling of $X$ and $Y$ is a stochastic process $(\widetilde{X},\widetilde{Y})$ on the product state space $\mathcal{X}\times\mathcal{X}$ such that

(i) The marginal processes $\widetilde{X}$ and $\widetilde{Y}$ are Markov processes with the transition kernel $P$;

(ii) If $\widetilde{X_{s}}=\widetilde{Y_{s}}$, we have $\widetilde{X_{t}}=\widetilde{Y_{t}}$ for all $t>s$.

The first meeting time of $X_{t}$ and $Y_{t}$ is denoted as $\tau^{C}:=\inf_{t\geq 0}\{X_{t}=Y_{t}\}$, which is called the coupling time. The coupling $(\widetilde{X},\widetilde{Y})$ is said to be successful if the coupling time is almost surely finite, i.e. $\mathbb{P}[\tau^{C}<\infty]=1$.

In order to give an estimate of the sensitivity of the QSD, we need the following two metrics.

In this paper, without further specification, we assume that the 1-Wasserstein distance is induced by $d(x,y)=\max\{1,\|x-y\|\}.$

Consider a Markov coupling $(\widetilde{X},\widetilde{Y})$ where $\widetilde{X}_{t}$ and $\widetilde{Y}_{t}$ are two numerical trajectories of the stochastic differential equation described in (2.5). Theoretically, there are many ways to make stochastic differential equations couple. But since numerical computation always has errors, two numerical trajectories may miss each other when the true trajectories couple. Hence we need to apply a mixture of the following coupling methods in practice.

Independent coupling. Independent coupling means the noise term in the two marginal processes $X_{t}$ and $Y_{t}$ are independent when running the coupling process $(\widetilde{X},\widetilde{Y})$. That is

where $(W^{(1)}_{(n+1)\delta}-W^{(1)}_{n\delta})$ and $(W^{(2)}_{(n+1)\delta}-W^{(2)}_{n\delta})$ are independent random variables for each $n$.

Reflection coupling Two Wiener processes meet less often than the 1D case when the state space has higher dimensions. This fact makes the independent coupling less effective. The reflection coupling is introduced to avoid this case. Take the Euler-Maruyama scheme of the SDE

as an example, where $\sigma$ is a constant matrix. The Euler-Maruyama scheme of $X_{t}$ reads as

where $W$ is a standard Wiener process. The reflection coupling means that we run $\hat{X}^{\delta}_{n}$ as

while run $\hat{Y}^{\delta}_{n}$ as

where $P=I-2e_{n}e^{T}_{n}$ is a projection matrix with

Nontechnically, reflecting coupling means that the noise term is reflected against the hyperplane that orthogonally passes the midpoint of the line segment connecting $\hat{X}^{\delta}_{n}$ and $\hat{Y}^{\delta}_{n}$. In particular, $e_{n}=-1$ when the state space is 1D.

Maximal coupling Above coupling schemes can bring $\hat{X}^{\delta}_{n}$ moves close to $\hat{Y}^{\delta}_{n}$ when running numerical simulations. However, a mechanism is required to make $\hat{X}^{\delta}_{n+1}=\hat{Y}^{\delta}_{n+1}$ with certain positive probability. That’s why the maximal coupling is involved. One can couple two trajectories whenever the probability distributions of their next step have enough overlap. Denote $p^{(x)}(x)$ and $p^{(y)}(x)$ as the probability density functions of $\hat{X}_{n+1}$ and $\hat{Y}_{n+1}$ respectively. The implementation of the maximal coupling is described in the following algorithm.

For discrete-time numerical schemes of SDEs, we use reflection coupling when $\hat{X}^{\delta}_{n}$ and $\hat{Y}^{\delta}_{n}$ are far away from each other, and maximal coupling when they are sufficiently close. The threshold of changing coupling method is $2\sqrt{\delta}\|\sigma\|$ in our simulation, that is, the maximal coupling is applied when the distance between $\hat{X}^{\delta}_{n}$ and $\hat{Y}^{\delta}_{n}$ is smaller than the threshold.

Let $\hat{X}^{\delta}=\{\hat{X}^{\delta}_{n},\ n\in\mathbb{Z}_{+}\}$ be a long numerical trajectory of $X_{t}$ as described in Algorithm 1. Let ${\bm{\tau}}=\{\tau_{m}\}^{M}_{m=0}$ be recordings of killing times of the numerical trajectory such that $X_{t}$ hits $\partial\mathcal{X}$ at $\tau_{0},\tau_{0}+\tau_{1},\tau_{0}+\tau_{1}+\tau_{2},\cdots$ when running Algorithm 1. Note that ${\bm{\tau}}$ is an 1D vector and each element in ${\bm{\tau}}$ is a sample of the killing time. It is well known that if the QSD $\mu$ exists for a Markov process, then there exists a constant $\lambda>0$ such that

Therefore, we can simply assume that the killing times ${\bm{\tau}}$ be exponentially distributed and the rate can be approximated by

One pitfall of the previous approach is that Algorithm 1 only gives a QSD when the time approaches to infinity. It is possible that ${\bm{\tau}}$ has not converged close enough to the desired exponential distribution. So it remains to check whether the limit is achieved. Our approach is to check the exponential tail in a log-linear plot. After having ${\bm{\tau}}$, it is easy to choose a sequence of times $t_{0},t_{1},\cdots,t_{n}$ and calculate $n_{i}=|\{\tau_{m}>t_{i}\,|\,0\leq m\leq M\}|$ for each $i=0,\cdots n$. Then $p_{i}=n_{i}/M$ is an estimator of $\mathbb{P}_{\mu}[\tau>t_{i}]$. Now let $p_{i}^{u}$ (resp. $p_{i}^{l}$) be the upper (resp. lower) bound of the confidence interval of $p_{i}$ such that

where $z=1.96$, $\tilde{n}_{i}=n_{i}+z^{2}$ and $\tilde{p}=\frac{1}{\tilde{n}}(n_{i}+\frac{z^{2}}{2})$ [1]. If $p_{i}^{l}\leq e^{-\lambda t_{i}}\leq p_{i}^{u}$ for each $0\leq i\leq n$, we accept the estimate $\lambda$. Otherwise we need to run Algorithm 1 for longer time to eliminate the initial bias in ${\bm{\tau}}$.

The data driven solver for the Fokker-Planck equation introduced in [16] can be modified to solve the QSD for the stochastic differential equation (2.3). We use the same 2D setting in Section 2.3 to introduce the algorithm. Let the domain $D$ and the boxes $\{O_{i,j}\}_{i=1,j=1}^{i=N,j=M}$ be the same as defined in Section 2.3. Let u be a vector in $\mathbb{R}^{N\times M}$ such that $u_{ij}$ approximates the probability density function at the center of the box $O_{i,j}$. As introduced in [7], we consider u as the solution to the following linear system given by the spatial discretization of the Fokker-Planck equation (2.4) with respect to each center point:

where $\mathbf{A_{0}}$ is an $(N-2)(M-2)\times(NM)$ matrix, which is called the discretized Fokker-Planck operator, and $\lambda$ is the killing rate, which can be obtained by the way we mentioned in previous subsection. More precisely, each row in $\mathbf{A_{0}}$ describes the finite difference scheme of equation (2.4) with respect to a non-boundary point in the domain $D$.

Similar to [16], we need the Monte Carlo simulation to produce a reference solution v, which can be obtained via Algorithm 1 in Section 2. Let $\hat{X}^{\delta}=\{\hat{X}^{\delta}_{n}\}^{N}_{n=1}$ be a long numerical trajectory of time -$\delta$ sample chain of process $X_{t}$ produced by Algorithm 1, and let $\textbf{v}=\{v_{i,j}\}^{i=N,j=M}_{i=1,j=1}$ such that

It follows from the ergodicity of (2.3) that v is an approximate solution of equation (2.4) when the trajectory is sufficiently long. However, as discussed in [16], the trajectory needs to be extremely long to make v accurate enough. Noting that the error term of v has little spatial correlation, we use the following optimization problem to improve the accuracy of the solution.

The solution to the optimization problem (3.2) is called the least norm solution, which satisfies $\mathbf{u}=\mathbf{v}-\mathbf{A^{T}(AA^{T})^{-1}(Av)}$, within $\mathbf{A=A_{0}}-\lambda\mathbf{I}$. [16]

An important method called the block data-driven solver is introduced in [7], in order to reduce the scale of numerical linear algebra problem in high dimensional problems. By dividing domain $D$ into $K\times L$ blocks $\{D_{k,l}\}^{k=1,l=1}_{k=K,l=L}$ and discretizing the Fokker-Planck equation, the linear constraint on $D_{k,l}$ is

where $\mathbf{A_{k,l}}$ is an $(N/K-2)(M/L-2)\times(NM/KL)$ matrix. The optimization problem on $D_{k,l}$ is

where $\mathbf{v^{k,l}}$ is a reference solution obtained from the Monte-Carlo simulation. Then the numerical solution to Fokker-Planck equation (3.1) is collage of all $\{u_{k,l}\}^{k=1,l=1}_{k=K,l=L}$ on all blocks. However, the optimization problem ”pushes” most error terms to the boundary of domain, which makes the solution is less accurate near the boundary of each block. Paper
citedobson2019efficient introduced the overlapping block method and the shifting blocks method to reduce the interface error. The overlapping block method enlarges the blocks and set the interior solution restricted to the original block as new solution, while the shifting block method moves the interface to the interior by shifting all blocks and recalculate the solution.

Note that in Section 3.1, we assume that $\lambda$ is a pre-determined value given by the Monte Carlo simulation. Theoretically one can also search for the minimum of $\|\mathbf{u-v}\|^{2}$ with respect to both $\lambda$ and v. But empirically we find that the result is not as accurate as using the killing rate $\lambda$ from the Monte Carlo simulation, possibly because $\mathbf{v}$ has too much error.

One natural question is that how the simulation error in $\lambda$ would affect the solution $\mathbf{u}$ to the optimization problem (3.2). Some linear algebraic calculation shows that the optimization problem (3.2) is fairly robust against small change of $\lambda$.