A convergent discretisation method for transition path theory for diffusion processes

TL;DR

This paper introduces a Monte Carlo discretisation method for transition path theory in diffusion processes that avoids Markovian assumptions, providing rigorous error bounds and convergence guarantees based on Voronoi tessellations.

Contribution

The authors develop a novel Monte Carlo approach for TPT that does not rely on Markovian approximations and establish theoretical error bounds and convergence results.

Findings

Method accurately approximates TPT objects from trajectory data.

Error bounds are rigorously proven and demonstrate convergence.

Application to a triple-well potential illustrates the method's effectiveness.

Abstract

Transition path theory (TPT) for diffusion processes is a framework for analysing the transitions of multiscale ergodic diffusion processes between disjoint metastable subsets of state space. Most methods for applying TPT involve the construction of a Markov state model on a discretisation of state space that approximates the underlying diffusion process. However, the assumption of Markovianity is difficult to verify in practice, and there are to date no known error bounds or convergence results for these methods. We propose a Monte Carlo method for approximating the forward committor, probability current, and streamlines from TPT for diffusion processes. Our method uses only sample trajectory data and partitions of state space based on Voronoi tessellations. It does not require the construction of a Markovian approximating process. We rigorously prove error bounds for the approximate TPT objects and use these bounds to show convergence to their exact counterparts in the limit of arbitrarily fine discretisation. We illustrate some features of our method by application to a process that solves the Smoluchowski equation on a triple-well potential.

Equations

$$S_i:=\{x\in S:{\left|x-g_i\right|}_2\le {\left|x-g_j\right|}_2,j\ne i\}.$$

$$\text{cl}(S)=\bigcup _{i=1}^n{S}_i,S_i\cap S_j=\partial S_i\cap \partial S_j,i\ne j,$$

$$\rho (\{S_i{\}}_{i\in I}):=\underset{i\in I}{\sup }\text{diam}(S_i).$$

$$Y_t=\{\begin{array}{ll}i& X_t\in \text{int}(S_i),\\ j& \exists ϵ>0\text{ s.t. }\forall s\in (t-ϵ,t),X_s\in \text{int}(S_j).\end{array}$$

$${\tau }_{A\cup B}(X):=\inf \{t\ge 0:X_t\in A\cup B\}.$$

$${\tau }_{J\cup K}(Y):=\inf \{t\ge 0:Y_t\in J\cup K\}.$$

$$q:S\to [0,1],x\mapsto q(x):=\mathbb{P}(X_{{\tau }_{A\cup B}(X)}\in B|X_0=x).$$

$${\tilde{q}}_i:=\mathbb{P}\left(Y_{{\tau }_{J\cup K}(Y)}\in K|Y_0=i\right).$$

$$\tilde{q}(x):=\{\begin{array}{ll}{\tilde{q}}_i& x\in \text{int}(S_i)\\ a& x\notin {\cup }_{i\in I}\text{int}(S_i),\end{array}$$

$${\hat{q}}_i:=\frac{1}{\mu (S_i)}⟨q,1_{S_i}{⟩}_{\mu }.$$

$$\hat{q}(x):=\{\begin{array}{ll}{\hat{q}}_i& x\in \text{int}(S_i)\\ b& x\notin {\cup }_{i\in I}\text{int}(S_i),\end{array}$$

$$\Vert q-\hat{q}{\Vert }_{L^p(\mu )}\le C\rho ,$$

$$\left|q(x)-q(y)-⟨\nabla q(x),x-y⟩\right|\le H{\left|x-y\right|}_2^2\le H{\left|x-y\right|}_2.$$

$$\Vert \left(q-\hat{q}\right){|}_{S_i}{\Vert }_{L^p(\mu )}^p$$

$$={\int }_{S_i}{\left|q(x)-q(x_i)-⟨\nabla q(x_i),x-x_i⟩+⟨\nabla q(x_i),x-x_i⟩\right|}^p\mu (\mathrm{d}x)$$

$$\le 2^{p-1}\left({\int }_{S_i}{\left(H{\left|x-x_i\right|}_2\right)}^p\mu (\mathrm{d}x)+{\int }_{S_i}{\left|⟨\nabla q(x_i),x-x_i⟩\right|}^p\mu (\mathrm{d}x)\right)$$

$$\le 2^{p-1}\left({\int }_{S_i}{\left(H{\left|x-x_i\right|}_2\right)}^p\mu (\mathrm{d}x)+{\int }_{S_i}{\left({\Vert q\Vert }_{\text{Lip}}{\left|x-x_i\right|}_2\right)}^p\mu (\mathrm{d}x)\right)$$

$$\le 2^{p-1}{\rho }^p\mu (S_i)\left(H^p+{\Vert q\Vert }_{\text{Lip}}^p\right).$$

$$\Vert q-\hat{q}{\Vert }_{L^p(\mu )}^p$$

$$\Vert q-\tilde{q}{\Vert }_{L^p(\mu )}\le C\rho$$

$$\begin{array}{r}\underset{s\to 0^{+}}{\lim }\frac{1}{s}\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_{\mathcal{R}\cap \left[0,T\right]}{1}_C\left(X_t\right)1_{C^{\complement }}\left(X_{t+s}\right)-1_{C^{\complement }}\left(X_t\right)1_C\left(X_{t+s}\right)\mathrm{d}t\\ ={\int }_{\partial C}{J}_{AB}(y)\cdot n_C(y){\sigma }_C(\mathrm{d}y)\end{array}$$

$${\int }_{\partial S_i}{J}_{AB}(y)\cdot n_{S_i}(y){\sigma }_{S_i}(\mathrm{d}y)=\sum _{k\in \mathcal{N}(S_i)}{\int }_{\partial S_i\cap \partial S_k}{J}_{AB}(y)\cdot n_{ik}(y){\sigma }_{\partial S_i\cap \partial S_k}(\mathrm{d}y),$$

$$\mathcal{N}(S_i):=\{j\in I:\dim S_i\cap S_j=d-1\},$$

$$\begin{array}{rl}{\alpha }_{ik}:=& \underset{s\to 0^{+}}{\lim }\frac{1}{s}\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_{\mathcal{R}\cap \left[0,T\right]}{1}_{S_i}\left(X_t\right)1_{S_k}\left(X_{t+s}\right)-1_{S_k}\left(X_t\right)1_{S_i}\left(X_{t+s}\right)\mathrm{d}t\\ =& {\int }_{\partial S_i\cap \partial S_k}{J}_{AB}(y)\cdot n_{ik}{\sigma }_{\partial S_i\cap \partial S_k}(\mathrm{d}y).\end{array}$$

$${\hat{\alpha }}_{ik}:=\frac{{\alpha }_{ik}}{\sigma (\partial S_i\cap \partial S_k)}={\tilde{J}}_{AB,i}\cdot n_{ik},k\in \mathcal{N}(S_i).$$

$${\hat{\alpha }}_i=N_i{\tilde{J}}_{AB,i}\in {\mathbb{R}}^{\text{card}(\mathcal{N}(S_i))}.$$

$$M_i{\tilde{J}}_{AB,i}={\beta }_i,M_i:=N_i^{\top }{N}_i\in {\mathbb{R}}^{d\times d},{\beta }_i:=N_i^{\top }{\hat{\alpha }}_i\in {\mathbb{R}}^d$$

$${\tilde{J}}_{AB}(x)={\tilde{J}}_{AB,i(x)},i(x):=\min \left\{\text{argmax}\left\{{\left|{\tilde{J}}_{AB,i}\right|}_2|x\in S_i\right\}\right\}.$$

$${\sigma }_{\min }(S_i):={\sigma }_{\min }(N_i)=\sqrt{{\sigma }_{\min }(M_i)},{\sigma }_{\max }(S_i):={\sigma }_{\max }(N_i)=\sqrt{{\sigma }_{\max }(M_i)}.$$

$${\Vert {\tilde{J}}_{AB}-J_{AB}\Vert }_{L^2(S\setminus (A\cup B),\mu ;{\mathbb{R}}^d)}\le C{\Vert J_{AB}\Vert }_{\text{Lip}}\rho ,$$

$${\left|{\tilde{J}}_{AB,i(x)}-J_{AB}(x)\right|}_2\le \frac{{\sigma }_{\max }(S_{i(x)})}{{\sigma }_{\min }^2(S_{i(x)})}\sqrt{\text{card}(\mathcal{N}(S_{i(x)}))}{\Vert J_{AB}\Vert }_{\text{Lip}}\rho ,$$

$${\left|{\tilde{J}}_{AB,i(x)}-J_{AB}(x)\right|}_2$$

Full text

1 Introduction

In many applications, one is often interested in understanding the rare transitions of a stochastic dynamical system between two metastable subsets of state space. Here, ‘metastable set’ refers to a set in which the system spends a long time on average before exiting, and ‘rare’ means that the probability of observing a transition between any two metastable sets is small. The multiscale behaviour of such a system is captured by the fact that the system spends very short times outside the union of metastable sets and spends long times inside this union.

Transition path theory (TPT) is a framework for the analysis of transition events of multiscale ergodic diffusion processes between two metastable subsets of state space [7, 8]. TPT has played an important role in the research area of molecular dynamics, where the conformational changes of a molecule occur over much longer time scales (milliseconds or longer) and fluctuations in bond lengths occur over very short time scales (nano- or femtoseconds). The conformational changes are of interest because they play an important role in chemical reactions, e.g. in cellular processes. A far from exhaustive list of references in this area of research is [14, 18, 31, 3, 27, 2, 20, 21, 28, 29, 16]. In addition to molecular dynamics, TPT has been applied to drug design in pharmacology [26], to study flows in complex networks [5], to study the temporal behavior of online topics on social networks [12], and to model atmospheric circulation in climate events [25]. Versions of TPT for Markov jump processes and Markov chains [17] and non-ergodic jump processes [32] have been developed. TPT has also been studied using the tools of stochastic analysis [13].

We now briefly outline the main concepts and objects of TPT for ergodic diffusion processes, and refer the reader to [7, 8, 16] for further details. For brevity, we will refer to ergodic diffusion processes simply as ‘processes’ in what follows. Let $A$ and $B$ be two disjoint, simply connected, open subsets of the state space $S$ of a process $X$. Given an infinitely long trajectory of $X$, a finite segment of the trajectory is said to be ‘reactive’ if, after leaving $A$, it enters $B$ before re-entering $A$. The probability that a trajectory of $X$ is reactive, conditioned on its current state being $X_{0}=x$, is given by the (forward) committor $q(x):=\mathbb{P}(X_{\tau_{A\cup B}(X)}\in B|X_{0}=x)$. Here, $\tau_{A\cup B}(X)$ denotes the first hitting time of $X$ with respect to the set $A\cup B$.

Under certain hypotheses on the drift and diffusion coefficient of the process, one can show that the committor is the unique solution of a Feynman-Kac or backward Kolmogorov partial differential equation (PDE) defined by the infinitesimal generator of the process and certain Dirichlet boundary conditions [8, Eq. (10)]. Consequently, the committor is twice continuously differentiable in the interior of $S\setminus(A\cup B)$. This exemplifies a general principle in probability theory, i.e. that objects that describe the averaged behaviour of the process exhibit certain regularity or smoothness properties.

Using the above-mentioned principle, we can ask whether there exist smooth curves that represent the average behavior of reactive trajectories of the process, in the sense that for any $\epsilon>0$, the $L^{\infty}$ ball of radius $\epsilon$ centered at such a curve contains a reactive trajectory. In TPT these smooth curves are called ‘streamlines’; they are integral curves of a vector field called the ‘probability current’. Using the streamlines, one can construct so-called ‘transition tubes’, by choosing a subset $A^{\prime}$ on the boundary of $A$ and taking the union of the streamlines with initial conditions in $A^{\prime}$. Given a probability measure on state space, e.g. the invariant Boltzmann-Gibbs measure, the probability that the process will transition from $A$ to $B$ via this transition tube can be computed from the measure of $A^{\prime}\subset\partial A$ with respect to the surface measure on the boundary of $A$ that is induced by the measure on state space. The ‘dominant transition pathways’ from $A$ to $B$ are the transition tubes with the highest probability. These pathways are important in applications, because they describe the most likely routes that will be taken by a multiscale stochastic dynamical system during transitions between metastable sets.

If the process is the solution to a Smoluchowski equation, i.e. a stochastic differential equation (SDE) where the drift coefficient is given by the gradient of a potential and the diffusion coefficient is proportional to the temperature, then one can additionally show that the probability current is proportional to the product of the Boltzmann-Gibbs density of the invariant measure with the gradient of the committor. If the state space is a subset of $\mathbb{R}^{d}$ for $d\in\{1,2,3\}$, then one can in principle completely avoid the task of collecting statistics of reactive trajectories, and apply deterministic numerical methods to compute the TPT objects described above. This can lead to useful visualisations of the behaviour of reactive trajectories, see e.g. [16]. However, in many applications, the high dimension of the state space makes deterministic numerical methods for solving PDEs impractical. As a result, Monte Carlo-based methods are often used when computing TPT objects in high dimensions.

In many applications of TPT, one often uses a Markovian stochastic process to approximate the underlying diffusion process, where the state space of the Markovian approximation is obtained by discretising the state space of the diffusion process [17]. The transition probabilities of the Markov approximation are estimated using short trajectory data. In molecular dynamics, these approximations are known as ‘Markov state models’; see e.g. [24] for a mathematical treatment of this subject and [18] for an application of TPT to a Markov state model. Since projecting the trajectories of the diffusion process onto the discrete state space results in a loss of the Markov property of the diffusion process, a correction step is needed to produce the Markov state model. The correction is often done by identifying a suitable ‘lag time’, i.e. an observation time window which yields an approximately Markovian process on the discrete state space; see [22] for an illustration. Although error analyses of Markov state models have been developed, e.g. in [23, 24], the justification for the lag time selection is in practice heuristic, and the approximating process is non-Markovian. To the best of our knowledge, there is no error analysis for methods that use a non-Markovian process on a discrete state space to approximate TPT objects for a diffusion process.

1.1 Contributions and outline

The goal of this paper is to develop a new computational method for approximating certain objects of TPT for an ergodic diffusion process $X$. The inputs to this method are a discretisation of the state space $S$ of the underlying process $X$ and sample trajectories of $X$. The outputs are approximations of the committor, probability current, and streamlines of $X$. We develop error bounds for each approximation and use these error bounds to prove convergence of each approximate TPT object to the corresponding TPT object of $X$ in the ‘continuum limit’ where the discretisation becomes arbitrarily uniformly fine.

We highlight some important features of our method. First, because our method computes the approximate TPT objects using only sample trajectories of $X$, it is not necessary to know or estimate the drift and diffusion coefficients of the SDE that defines the underlying process. It suffices to have a ‘black box’ that generates the trajectory data. In this sense, our approach is ‘data-driven’.

Second, our method does not involve constructing a Markovian stochastic process to approximate the underlying process $X$. In particular, we can avoid the task of choosing a suitable lag time, and the use of heuristic arguments to justify the choice of a lag time. This feature differentiates our method from methods that perform TPT that use Markov processes on discrete state spaces to approximate the underlying diffusion process of interest.

Third, our method comes with a rigorous error analysis. While there exist rigorous error analyses for Markov state models in the framework of transfer operators on continuous state spaces [23], there are no rigorous error analyses for TPT for Markov processes. We prove bounds on the approximation error of each TPT object in terms of how fine the discretisation is.

The paper proceeds as follows. We describe our discretisation methods in Section 2 and describe how these lead to a non-Markovian jump process. In Section 3 we define the approximate committor, the approximate probability current, and the approximate streamlines of our method. The approximate committor and the approximate probability current are piecewise constant functions on the continuous state space $S$ of the underlying process $X$, while the approximate streamline is a piecewise linear trajectory in $S$. For each approximate object, we prove a bound on the approximation error with respect to the corresponding TPT object for the underlying process $X$, in terms of the largest diameter of the sets in the discretisation of the continuous state space $S$. We then use these error bounds to prove convergence as the sets in the discretisation shrink to points; this is the above-mentioned ‘continuum limit’. In Section 4 we present numerical results, in which we compare the performance of our approach with the finite differences method for TPT for diffusion processes. We conclude in Section 5.

2 Setup

In what follows, we follow the setup in [8]. Let $X=\{X_{t}\}_{t\geq 0}$ denote the underlying diffusion process of interest. The state space of $X$ is an open, bounded, $d$-dimensional, simply connected set $S\subset\mathbb{R}^{d}$, and the boundary of $S$ is sufficiently regular for us to impose reflecting boundary conditions. We assume that $X$ is ergodic with respect to a probability measure $\mu$, that $\mu$ admits a density with respect to Lebesgue measure on $S$, and that the Lebesgue density of $\mu$ is strictly positive on $S$. We write $A$ and $B$ to denote two open, simply connected subsets of $S$ whose closures are disjoint.

Given an arbitrary set $U\subset S$, let $\textup{int}({U})$, $\textup{cl}({U})$, $\partial U$, $U^{\complement}$ and $\mathbf{1}_{U}$ denote the interior, closure, boundary, complement with respect to $S$ (i.e. $S\setminus U$), and indicator function of $U$ respectively. The cardinality of an arbitrary set $V$ is written $\textup{card}({V})$, and the restriction of a function $f$ to $V$ is written $f|_{V}$. The transpose of a vector $v\in\mathbb{R}^{d}$ is written $v^{\top}$; the $\ell_{p}$ norm for $1\leq p\leq\infty$ of $v$ is written $\left|v\right|_{p}$. The Lipschitz constant of an arbitrary function $f$ by $\left\|f\right\|_{\textup{Lip}}$. Given a finite measure $\nu$ on $S$, we denote the $L^{p}(\nu)$ norm on a space of functions defined on $S$ or a subset of $S$ by $\left\|\cdot\right\|_{L^{p}(\nu)}$. We fix a probability space $(\Omega,\mathcal{F},\mathbb{P})$ that is assumed to be large enough to admit all random variables considered below.

2.1 Discretisation

We review the relevant concepts from Voronoi tessellations and introduce some important terms concerning the discretisations that we will use in the sequel. Our choice of Voronoi tessellations as a discretisation method is motivated by their successful application in not only molecular dynamics [11, 30] but many other applications; see [1, 6] and the references therein. The key property of Voronoi tessellations that we use are that Voronoi tessellations partition a set into polytopes that overlap at most at their boundaries, and that adjacency relations can be computed without computing the full Voronoi tessellation, i.e. without computing the vertices of every polytope in the tessellation.

A Voronoi tessellation of $S$ associated to a finite set of generators $\{g_{1},\ldots,g_{n}\}\subset S$ for some $n\in\mathbb{N}$ is a collection $\{S_{1},\ldots,S_{n}\}$ of nonempty subsets of $S$, where each subset is the Voronoi cell defined by

[TABLE]

Thus, $S_{i}$ is the closed set consisting of all points in state space that are closer in the Euclidean metric to the generator $g_{i}$ than to any other generator. It can be shown that each Voronoi cell $S_{i}\subset\mathbb{R}^{d}$ is a $d$-polytope, i.e. a bounded, closed, convex subset such that $\dim{\textup{int}({S_{i}})}=d$, that is defined as the intersection of finitely many half-spaces or equivalently as the convex hull of finitely many extreme points. The boundedness of each cell $S_{i}$ follows from the boundedness of the state space $S$. Since every Voronoi cell $S_{i}$ is a neighbourhood of its generator $g_{i}$, it has strictly positive Lebesgue measure. The Voronoi cells form a partition of the state space $S$, i.e.

[TABLE]

so that Voronoi cells intersect at most at their boundaries, which have zero Lebesgue measure. We declare two distinct Voronoi cells $S_{i}$ and $S_{j}$ to be adjacent if they share a common facet, i.e. if $\dim(S_{i}\cap S_{j})=d-1$. Given a Voronoi tessellation $\{S_{i}\}_{i\in I}$, $I=\{1,\ldots,n\}$, its dual Delaunay graph is the graph $G=(I,E)$ with vertex set $I$ and edge set $E$ consisting of all pairs $(i,j)$ such that $S_{i}$ and $S_{j}$ are adjacent.

Recall that, given a nonempty set $A\subset\mathbb{R}^{d}$, the Euclidean diameter of $A$ is defined by $\textup{diam}({A})=\sup\{\left|x-y\right|_{2}\ :\ x,y\in A\}$. This leads to the next definition.

Definition 2.1**.**

The width $\rho$ of a Voronoi tessellation $\{S_{i}\}_{i\in I}$ is the largest Euclidean diameter of the Voronoi cells, i.e.

[TABLE]

When there is no risk of confusion, we will omit the argument $\{S_{i}\}_{i\in I}$ of the width and simply write $\rho$. The width provides uniform control over the Euclidean diameter of all the cells in the tessellation. The smaller (respectively larger) the width, the finer (resp. coarser) the tessellation. We shall be interested in obtaining error bounds in the continuum limit, i.e. as $\rho$ decreases to [math]. In the continuum limit, the cells in the Voronoi tessellation shrink to points.

Remark 2.2**.**

An important property of Voronoi tessellations is that, if one knows that two Voronoi cells $S_{i}$ and $S_{j}$ that are generated respectively by $g_{i}$ and $g_{j}$ are adjacent, then the outer unit normal to $S_{i}$ on the facet $S_{i}\cap S_{j}=\partial S_{i}\cap S_{j}$ is given by $(g_{j}-g_{i})/\left|g_{i}-g_{j}\right|_{2}$. This is because the facet $S_{i}\cap S_{j}$ is contained in the hyperplane of points that are equidistant from $g_{i}$ and $g_{j}$ and $g_{j}-g_{i}$ is normal to this hyperplane.

2.2 Definition of approximating jump process

Recall that many applications of TPT involve the construction of a Markovian stochastic process on a discrete state space that approximates the underlying diffusion process on the continuous state space. In this section, we construct a continuous-time process on the Delaunay graph $G=(I,E)$ mentioned above. We explain why this process is not Markovian and show that this process preserves an important property of the underlying process $X$ on $S$.

Let $\{S_{i}\}_{i\in I}$ be a Voronoi tessellation with finite index set $I$, and let $G=(I,E)$ be the Delaunay graph associated to $\{S_{i}\}_{i\in I}$. We define a continuous time jump process $Y=\{Y_{t}\}_{t\geq 0}$ on $G$ by ‘projecting’ the trajectory of the underlying process $X$ to $I$ according to

[TABLE]

The second case in (2.1) can be interpreted as follows: if $X_{t}$ lies on the boundary of a Voronoi cell, then we assign to $Y_{t}$ the state $j\in I$, where $j$ is the index of the set whose interior contained the trajectory of $Y$ in the most recent past, i.e. up to but not including the current time $t$. The memory effect implied by the second case in (2.1) of $Y$ implies that $Y$ is not Markovian.

In order to define the analogues of the sets $A$ and $B$ in the discrete state space $I$, we make the following assumption.

Assumption 2.3**.**

Given the sets $A,B\subset S$, there exist disjoint subsets $J,K\subset I$ satisfying $\textup{cl}({A})=\cup_{j\in J}S_{j}$ and $\textup{cl}({B})=\cup_{k\in K}S_{k}$.

The assumption implies that $J$ and $K$ are the ‘metastable subsets’ for the jump process $Y$. The assumption effectively amounts to making the definition of $A$ and $B$ dependent on the partition $\{S_{i}\}_{i\in I}$. Although this may appear to be problematic at first, we can justify it by both theoretical and practical considerations. The theoretical consideration is that we are interested in the continuum limit, in which case the cells $\{S_{i}\}_{i\in I}$ shrink to points. Thus, even if the assumption is not satisfied for a particular value of the width $\rho$, it will be satisfied in the continuum limit, regardless of the geometry of $A$ or $B$. The practical consideration is that in many applications, there is some flexibility in the definition of the metastable sets $A$ and $B$ in state space. For example, in molecular dynamics, one may choose $A$ and $B$ to each be a set of sufficiently small diameter that contains a local minima of the energy landscape. The precise geometry of each set is less important relative to the preservation of the metastability property. Thus it is reasonable to adapt the definitions of $A$ and $B$ according to the sets in the discretisation, since this facilitates the determination of when a process has entered a metastable set.

Recall the definition of the first hitting time of $X$ with respect to the set $A\cup B$,

[TABLE]

Define the first hitting time of the jump process $Y$ with respect to the set $J\cup K$ as

[TABLE]

In practice, it is not possible to detect exactly when $Y$ hits $J$ or $K$, because numerical methods for generating trajectories of the underlying process $X$ will produce only approximations of the trajectories of $X$, e.g. by linear interpolation between points generated by the Euler-Maruyama method. Thus, using a numerically generated trajectory of $X$ in (2.1) will generate only an approximation of $Y$. However, for the sake of simplicity, we shall not take into account the error incurred by these approximations, and we shall assume that we can detect when $Y$ hits $J$ or $K$.

The preceding definitions lead to the following result.

Lemma 2.4**.**

Suppose that Assumption 2.3 holds, and let $\tau_{A\cup B}(X)$ and $\tau_{J\cup K}(Y)$ be defined as in (2.2) and (2.3) respectively. Then $\tau_{A\cup B}(X)=\tau_{J\cup K}(Y)$.

Proof.

Fix an arbitrary sample trajectory $X(\omega)$ of $X$, and let $t:=\tau_{A\cup B}(X(\omega))$. Given that $A$ and $B$ are open, and given the definition (2.2) of $\tau_{A\cup B}(X)$, it follows that there exist $\epsilon_{1},\epsilon_{2}>0$ such that $X_{s}(\omega)\notin A\cup B$ for $s\in(t-\epsilon_{1},t)$, $X_{t}(\omega)\in\partial A\cup\partial B$, and $X_{s}(\omega)\in A\cup B$ for $s\in(t,t+\epsilon_{2})$.

Suppose first that $X_{t}(\omega)\in\partial A$ and $X_{s}(\omega)\in A$ for $s\in(t,t+\epsilon_{2})$. Let $i\in I$ be such that $X_{s}(\omega)\in\textup{int}({S_{i}})$ for $s\in(t-\epsilon_{1},t)$. By the definition (2.1) of $Y$, it follows that $Y_{s}(\omega)=i$ for $s\in(t-\epsilon_{1},t)$ and $Y_{t}(\omega)=i$. Furthermore, given Assumption 2.3, $X_{s}(\omega)\in A$ for $s\in(t,t+\epsilon_{2})$ implies that $Y_{s}(\omega)\in J$ for $s\in(t,t+\epsilon_{2})$. Since $\inf\{s\ :t<s<t+\epsilon_{2}\}=t$, it follows that $\tau_{J\cup K}(Y)=t=\tau_{A\cup B}(X)$. If instead $X_{t}(\omega)\in\partial B$ and $X_{s}(\omega)\in B$ for $s\in(t,t+\epsilon_{2})$, then the same argument holds after replacing $J$ with $K$. This completes the proof. ∎

Except for $\tau_{J\cup K}(Y)$, our method does not involve computing or estimating any properties of the non-Markovian jump process $Y$, such as its transition probabilities. We introduce the process $Y$ only to help the interpretation of the approximate committor that we define in §3.1.

3 Definition and convergence analysis of approximate TPT objects

In this section we define the approximate TPT objects that are the outputs of our method. We prove error bounds for each object and use these bounds to show that in the continuum limit, each approximate TPT object we define converges to the corresponding object for the underlying process $X$. For simplicity, in the error bounds we do not take into account statistical errors, i.e. errors due to the use of finite-sample Monte Carlo averages as approximations for expected values.

3.1 Approximate committor

Recall the definition of the committor of $X$,

[TABLE]

Recall that $K$ is the subset of discrete state space $I$ corresponding to the metastable subset $B$ of the continuous state space $S$. We define the approximate committor in an analogous way, using the jump process $Y$ defined in (2.1). For $i\in I$, define $\tilde{q}_{i}$ by

[TABLE]

In order to compute each $\tilde{q}_{i}$ using sample trajectories of $X$, we use a Monte Carlo approach: we sample finitely many trajectories, where the initial condition of each trajectory is an independent draw from the uniform distribution on $S_{i}$, and compute the empirical probability that a trajectory starting from $S_{i}$ enters $B$ before $A$.

We define the approximate committor on the state space $S$ of the underlying process $X$ as a function that is piecewise constant on the interiors of the cells $\{S_{i}\}_{i\in I}$ of the tessellation:

[TABLE]

where we may choose an arbitrary $a\in\mathbb{R}$, since $(\cup_{i\in I}\textup{int}({S_{i}}))^{\complement}$ has Lebesgue measure zero and since we will measure the error of $\hat{q}$ with respect to $q$ in the $L^{p}(\mu)$ norm. Up to a redefinition of the values on $(\cup_{i\in I}\textup{int}({S_{i}})))^{\complement}$, the approximate committor $\tilde{q}$ can be interpreted as a committor for $Y$, because the definition (3.3) involves the $\tilde{q}_{i}$’s from (3.2).

Recall that the ergodic measure $\mu$ is assumed to have a strictly positive Lebesgue density on $S$. Define the $\mu$-weighted inner product on $L^{2}(S,\mu;\mathbb{R})$ by $\langle v,w\rangle_{\mu}=\int_{S}v(x)w(x)\mu(\mathrm{d}x)$. Given a tessellation $\{S_{i}\}_{i\in I}$, we use the weighted inner product to define $\hat{q}_{i}$ for each $i\in I$ by

[TABLE]

Since each $S_{i}$ has positive Lebesgue measure, it follows that $\mu(S_{i})>0$ for each $i\in I$. We can define the projected committor function $\hat{q}:S\to[0,1]$ analogously to (3.3):

[TABLE]

for a fixed, arbitrary $b\in\mathbb{R}$. With these preparations, we may state the following proposition. The proof uses some technical lemmas that will not be used in the sequel, so we defer the proof to §A.1.

Proposition 3.1**.**

Suppose that Assumption 2.3 holds. Let $\hat{q}_{i}$ and $\tilde{q}_{i}$ be defined as in (3.4) and (3.2) respectively. If $\mathbb{P}\circ(X_{0})^{-1}=\mu$, then $\hat{q}_{i}=\tilde{q}_{i}$, for all $i\in I$.

The following lemma uses the convexity of Voronoi cells and the continuity of the committor function to prove that every Voronoi cell $S_{i}$ contains a point $x_{i}\in S_{i}$ such that $q$ attains the value $\hat{q}_{i}$ at $x_{i}$. We shall use this lemma later to prove Theorem 3.4.

Lemma 3.2**.**

Let $\{S_{i}\}_{i\in I}$ be a Voronoi tessellation of $S$, and let $\hat{q}_{i}$ be defined as in (3.4). For every $i\in I$, there exists some $x_{i}\in S_{i}$ such that $q(x_{i})=\hat{q}_{i}$.

Proof.

If $q$ is constant on $S_{i}$, then it must equal $\hat{q}_{i}$, and any $x_{i}\in S_{i}$ satisfies the desired property. Therefore, suppose that $q$ is not constant on $S_{i}$, and partition $S_{i}$ into the disjoint subsets $S^{-}_{i}:=\{x\in S_{i}\ :\ q(x)<\hat{q}_{i}\}$, $S^{+}_{i}:=\{x\in S_{i}\ :\ q(x)>\hat{q}_{i}\}$ and $S^{0}_{i}:=\{x\in S_{i}\ :\ q(x)=\hat{q}_{i}\}$. Since $q$ is continuous and not constant on $S_{i}$, it follows that $S^{-}_{i}$ and $S^{+}_{i}$ are nonempty. Let $x^{\prime}\in S^{-}_{i}$ and $x^{\prime\prime}\in S^{+}_{i}$. It follows from the intermediate value theorem that there exists a $t\in(0,1)$ such that $x_{i}(t):=(1-t)x^{\prime}+tx^{\prime\prime}$ satisfies $q(x_{i}(t))=\hat{q}_{i}$. Since $x^{\prime},x^{\prime\prime}\in S_{i}$ and since any Voronoi cell $S_{i}$ is convex, it follows that $x_{i}(t)$ belongs to $S_{i}$. ∎

Remark 3.3**.**

The conclusion of Lemma 3.2 holds for more general partitions $\{S_{i}\}_{i\in I}$ of state space $S$ for which each $S_{i}$ is path connected, since in this case we can replace the line segment $x_{i}(t):=(1-t)x^{\prime}+tx^{\prime\prime}$ with a curve in $S_{i}$ with endpoints $x^{\prime}$ and $x^{\prime\prime}$ and follow the same reasoning thereafter.

We prove an error bound for the error incurred when we approximate the true committor $q$ with the projected committor $\hat{q}$ defined in (3.5).

Theorem 3.4** (Error bound for projected committor).**

Suppose that the committor $q:S\to[0,1]$ is globally Lipschitz continuous and has bounded second-order derivatives. If $\{S_{i}\}_{i\in I}$ is a Voronoi tessellation with width $\rho\leq 1$, then for any $p\geq 1$,

[TABLE]

where $C>0$ depends only on $q$.

Proof.

Recall that the committor $q$ is continuously differentiable in $S\setminus(A\cup B)$, since it solves a backward Kolmogorov equation. Given the assumption that $q$ is globally Lipschitz continuous, it follows that $\sup_{x\in S}\left|\nabla q(x)\right|\leq\left\|q\right\|_{\textup{Lip}}$. In addition, given the assumption of bounded second-order derivatives, there exists some $H>0$ that depends only on $q$, such that for all $x,y\in S$ with $\left|x-y\right|_{2}\leq\rho$ and $\rho\leq 1$,

[TABLE]

Fix an arbitrary $p\in[1,\infty)$, and fix an arbitrary $i\in I$. By Lemma 3.2, there exists an $x_{i}\in S_{i}$ such that $q(x_{i})=\hat{q}_{i}=\hat{q}|_{S_{i}}$. We obtain

[TABLE]

Above, we used (3.6) and $(a+b)^{p}\leq 2^{p-1}(a^{p}+b^{p})$ in the first inequality, the global Lipschitz continuity of $q$ and the Cauchy-Schwarz inequality in the second inequality, and the assumption that $\{S_{i}\}_{i\in I}$ has width $\rho$ and the fact that $x,x_{i}\in S_{i}$ with $\textup{diam}({S_{i}})\leq\rho$ in the third inequality. Summing over $i\in I$, we obtain

[TABLE]

since $\mu(S)=1$, and since $a^{p}+b^{p}\leq(a+b)^{p}$, for $a,b\geq 0$ and $p\geq 1$. ∎

Corollary 3.5** (Error bound for discrete committor).**

Suppose that Assumption 2.3 and the assumptions of Theorem 3.4 hold. Then for any Voronoi tessellation $\{S_{i}\}_{i\in I}$ of $S$ with sufficiently small width $\rho$, the function $\tilde{q}$ defined in (3.3) satisfies

[TABLE]

for the same $C>0$ as in the conclusion of Theorem 3.4.

Proof.

The result follows from Theorem 3.4 and Proposition 3.1. ∎

3.2 Approximate probability current

In this section we define the second output of our method, namely the approximate probability current $\tilde{J}_{AB}$. We describe how $\tilde{J}_{AB}$ is obtained from sample reactive trajectories of $X$.

To define the approximate probability current, we first recall the definition of the probability current from TPT for diffusion processes, from the presentation in [15, Section 3.1.4]. The probability current is a function $J_{AB}:S\setminus(A\cup B)\rightarrow\mathbb{R}^{d}$ such that at any point $x\in S\setminus(A\cup B)$, $J_{AB}(x)$ represents the net flux of reactive trajectories from $A$ to $B$ through that point. For any $d$-dimensional $C\subset S\setminus(A\cup B)$ with boundary $\partial C$, $J_{AB}$ is defined implicitly via

[TABLE]

where $n_{C}(y)$ denotes the outward facing unit normal to $C$ at some point $y\in\partial C$, $\sigma_{C}$ denotes the surface measure on $\partial C$ and $\mathcal{R}$ indicates the set of reactive times, i.e. the union of time intervals in which the trajectory of the process $X$ is reactive.

To motivate our definition of the approximate probability current, we first observe that if we set the region $C$ in the formula above to be a Voronoi cell or more generally, a polytope $S_{i}$, then the unit normal $n_{C}$ in the integral on the right-hand side of (3.7) takes only finitely many values that are given by the outer unit normals to the facets of $S_{i}$. In particular, we have

[TABLE]

where $\mathcal{N}(S_{i})$ denotes the set of indices of cells adjacent to $S_{i}$, i.e.

[TABLE]

and $n_{ik}$ is the outer unit normal vector to $S_{i}$ in the relative interior of the facet $\partial S_{i}\cap\partial S_{k}$. Thus we arrive at a decomposition over facets of the right-hand side of (3.7).

We can obtain a similar decomposition over facets of the left-hand side of (3.7). Observe that if in a time interval of vanishing length a transition out of $C=S_{i}$ takes place, then almost surely with respect to normalised surface measure on $\partial S_{i}$, the transition is into an adjacent cell $S_{k}$. This is due to the almost sure continuity of sample paths of diffusion processes. Therefore, in the $s\to 0^{+}$ limit, $\mathbf{1}_{C^{\complement}}(X_{t+s})=\sum_{k\in\mathcal{N}(S_{i})}\mathbf{1}_{S_{k}}(X_{t+s})$. Thus, if $C=S_{i}$ for arbitrary $i\in I$, then we can write (3.7) as a system of $\textup{card}({\mathcal{N}(S_{i})})$ equations: for $k\in\mathcal{N}(S_{i})$,

[TABLE]

Following the interpretation of the probability current $J_{AB}(y)$ at $y$ as the net flux of reactive trajectories from $A$ to $B$ through $y$, we may interpret $\alpha_{ik}$ as the average net flux of reactive trajectories across the facet $S_{i}\cap S_{k}$.

We now describe how to obtain the approximate probability current. Under reasonable conditions, the probability current is continuous on its domain [15, Section 3.1.4, Eq. (3.15)]. If the probability current is continuous, then for sufficiently small $S_{i}$, we may replace the probability current $J_{AB}$ in the right-hand side of (3.9) with a constant vector $\tilde{J}_{AB,i}\in\mathbb{R}^{d}$. Since the choice of $\tilde{J}_{AB,i}$ should be consistent with (3.9), the resulting expression yields the system of equations for the unknown $\tilde{J}_{AB,i}$:

[TABLE]

Note that in order to compute the volume $\sigma(\partial{S_{i}}\cap\partial{S_{k}})$ of the facet $\partial{S_{i}}\cap\partial{S_{k}}$, we need the vertices of the facet. We will return to this observation in §5.

Let $\hat{\alpha}_{i}\in\mathbb{R}^{\textup{card}({\mathcal{N}(S_{i})})}$ denote the vector with entries given by $\hat{\alpha}_{ik}$, and let $N_{i}\in\mathbb{R}^{\textup{card}({\mathcal{N}(S_{i})})\times d}$ denote the matrix whose $k$-th row is given by $n_{ik}^{\top}$. Then we may rewrite the preceding system of equations as a matrix-vector product

[TABLE]

For (3.11) to admit a solution $\tilde{J}_{AB,i}$, we need that $\hat{\alpha}_{i}$ belongs to the column space of $N_{i}$. Since $\textup{card}({\mathcal{N}(S_{i})})\geq d+1$ (see Lemma A.6 below) this need not hold in general. Therefore we consider the system of normal equations

[TABLE]

and use the fact that for any $d$-dimensional polytope in $\mathbb{R}^{d}$, the matrix $N_{i}$ of outer unit normals to $S_{i}$ has full rank (see Corollary A.7). Thus $M_{i}$ is invertible, and we may solve (3.12) using standard methods from numerical linear algebra. By solving (3.12) for every $i\in I$, we obtain the approximate probability current $\tilde{J}_{AB}$ as a piecewise constant function on $S\setminus(A\cup B)$, by setting

[TABLE]

That is, if $x\in\textup{int}({S_{i}})$ then $\tilde{J}_{AB}(x)=\tilde{J}_{AB,i}$, and if $x\in\partial S_{i}$ for certain $i\in I$, then we set $\tilde{J}_{AB}(x)$ to be the $\tilde{J}_{AB,i}$ with largest $\ell_{2}$ norm and choose the vector with smallest index if there are two or more vectors with the same $\ell_{2}$ norm.

Remark 3.6**.**

For the error bound that we prove in this section, the values of the approximate probability current on the boundaries of cells are not important. This explains the arbitrary nature of the choice we made above. However, the property that the value is chosen from the set of approximate probability current vectors of cells whose boundaries contain $y$ ensures that the boundary values are reasonable.

We outline how to compute each entry of $\hat{\alpha}_{ik}$. Suppose we use a numerical method for generating trajectories of $X$ with time step $0<\Delta t\ll 1$, and suppose that we have generated one long trajectory of $X$ of duration $T:=M\Delta t$ for some $M\in\mathbb{N}$, $M\gg 1$. For a given $i\in I$ and $k\in\mathcal{N}(S_{i})$, we compute $\hat{\alpha}_{ik}$ by taking the difference between the total number of observed one-step transitions from $S_{i}$ to $S_{k}$ and the total number of observed one-step transitions from $S_{k}$ to $S_{i}$, where we count only those one-step transitions that occur within reactive trajectory segments of the long trajectory. Since we cannot observe state transitions that occur over intervals of length shorter than $\Delta t$, we set $s=\Delta t$ and divide this difference by $T$ and $s$ as in (3.9). Note that the $\Delta t$ term in $T=M\Delta t$ cancels the $\Delta t$ component arising from the $\mathrm{d}t$ term in the integral. The procedure just described yields $\alpha_{ik}$. We can compute the surface measure $\sigma(\partial S_{i}\cap\partial S_{k})$ using methods described in [4], and then use (3.10).

Modulo statistical error due to insufficient sampling and the error due to the approximate volume computations needed for $\sigma(\partial S_{i}\cap\partial S_{k})$, the computation method described above incurs errors due to the approximation of the $s\to 0^{+}$ limit and of the $T\to\infty$ limit by a fixed time step and duration respectively. These constraints are inevitable, given a finite computational budget. In particular, the error due to the approximation of the $T\to\infty$ limit is the finite sampling error that is common to all Monte Carlo methods.

For any real matrix $G$, we denote the smallest and largest nonzero singular value of $G$ respectively by $\sigma_{\min}(G)$ and $\sigma_{\max}(G)$. We define the smallest and largest singular value of a Voronoi cell $S_{i}$ as the smallest and largest singular value of the matrix $N_{i}$ of outer unit normals respectively:

[TABLE]

With this notation in mind, we can state and prove the following result.

Theorem 3.7**.**

Let $\{S_{i}\}_{i\in I}$ be a Voronoi tessellation of $S$ with width $\rho$ that satisfies Assumption 2.3. If the probability current $J_{AB}$ is globally Lipschitz, then

[TABLE]

for $C=C(\{S_{i}\}_{i\in I}):=\sup_{i\in I}\sigma_{\max}(S_{i})\sigma_{\min}^{-2}(S_{i})\textup{card}({\mathcal{N}(S_{i})})^{1/2}$.

Proof.

We first prove a pointwise error bound. For an arbitrary $x\in S\setminus(A\cup B)$,

[TABLE]

for the function $i:S\to I$ defined in (3.13). Recall that the matrix $M_{i}$ from (3.12) is invertible. Using (3.12) and (3.14), we obtain

[TABLE]

From (3.11) we obtain

[TABLE]

For any $i\in I$, $x\in S_{i}$, and $j\in\mathcal{N}(S_{i})$,

[TABLE]

Above, we used (3.9) and (3.11) for the equation, the Cauchy-Schwarz inequality for the first inequality, the Lipschitz continuity assumption of $J_{AB}$ and the fact that for every $i\in I$ and $j\in\mathcal{N}(S_{i})$ each $n_{ij}$ has unit norm in the second inequality, and the fact that $\{S_{i}\}_{i\in I}$ has width $\rho$ implies that for $x\in S_{i}$ and $y\in\partial S_{i}$, $|y-x|_{2}\leq\rho$ in the last inequality. Combining the preceding inequalities yields (3.16).

By Assumption 2.3, the union of the $S_{i}$ with $i\in I\setminus(J\cup K)$ covers $S\setminus(A\cup B)$, so

[TABLE]

Above, we use the fact that the boundaries of the cells have Lebesgue measure zero and hence $\mu$-measure zero, the fact the $\{\textup{int}({S_{i}})\}_{i\in I}$ are disjoint, and the fact that $x\in\textup{int}({S_{i}})$ implies $i(x)=i$. Thus (3.15) follows from (3.16), since $\mu(S)=1$. ∎

To establish a convergence result in the continuum limit, we need to specify a sequence $(\rho_{k})_{k\in\mathbb{N}}$ of widths that decreases to zero, and to impose uniform control over the corresponding sequence of tessellation-dependent coefficients $C$ from (3.15).

Corollary 3.8**.**

Let $(\rho_{k})_{k\in\mathbb{N}}\subset(0,\infty)$ satisfy $\rho_{k}\to 0$. For each $k\in\mathbb{N}$, let $\{S^{k}_{i}\}_{i\in I(k)}$ be a partition of $S$ with index set $I(k)$, width $\rho_{k}$, and corresponding approximate probability current $\tilde{J}^{k}_{AB}$, that satisfies Assumption 2.3. If there exists $0<C<\infty$ that does not depend on $k\in\mathbb{N}$ such that

[TABLE]

and if the probability current $J_{AB}$ is globally Lipschitz, then

[TABLE]

Proof.

For each $k\in\mathbb{N}$, Theorem 3.7 yields the upper bound on the error

[TABLE]

Since $\sup_{k\in\mathbb{N}}C_{k}\leq C$ by (3.17), the conclusion follows. ∎

A sequence of Voronoi tessellations that satisfies (3.17) is the sequence of tessellations where each sets in the tessellation is a translation and scaling of the hypercube in $\mathbb{R}^{d}$. In this case, for every $i\in I(k)$ we have that $\textup{card}({\mathcal{N}(S_{i})})=2d$, and the matrix of outer unit normals $N_{i}$ for each $S^{k}_{i}$ satisfies $N_{i}^{\top}N_{i}=2I_{d}$, where $I_{d}$ is the $d\times d$ identity matrix. However, since the number of sets in the resulting hypercubic partition of $S$ increases exponentially with the dimension $d$, such a tessellation would not be suitable when $d\gg 1$. The question of designing algorithms that produce discretisations satisfying (3.17) is interesting, but beyond the scope of this paper.

In this section, the only objects that require computing the vertices of every Voronoi cell are the volumes $\sigma(\partial S_{i}\cap\partial S_{k})$ of the facets. These volumes are needed in (3.10). By Remark 2.2, the full Voronoi tessellation is not necessary to compute the matrices $\{N_{i}\}_{i\in I}$ of outer unit normals.

3.3 Approximate streamlines

In §1, we described the streamlines as integral curves of the probability current $J_{AB}$. That is, a streamline starting at $s_{0}\in\partial A$ and ending at $\partial B$ is the solution $\{s(t)\}_{t\in[0,T(s_{0})]}$ of a differential equation with $J_{AB}$ as the driving vector field. We can write this streamline as

[TABLE]

Note that for every $s_{0}\in\partial A$, $T(s_{0})$ is finite. If this were not the case, then one can obtain a contradiction by using (3.7) to show that there exist reactive trajectories that do not reach $B$.

For a given tessellation $\{S_{i}\}_{i\in I}$ with index subsets $K$ and $J$ from Assumption 2.3 with corresponding approximate probability current $\tilde{J}_{AB}$ defined in (3.13), we define the approximate streamline with initial condition $\tilde{s}_{0}\in\partial(\cup_{i\in J}S_{i})$ by

[TABLE]

Note that $\tilde{T}(\tilde{s}_{0})$ need not be finite for every $\tilde{s}_{0}\in\partial(\cup_{i\in J}S_{i})$: we cannot apply the same reasoning as we did for showing that $T(s_{0})$ is finite, because $\tilde{J}_{AB}$ is only an approximation of $J_{AB}$.

Theorem 3.9**.**

Let $(\rho_{k})_{k\in\mathbb{N}}\subset(0,\infty)$ satisfy $\rho_{k}\to 0$. For each $k\in\mathbb{N}$, let $\{S^{k}_{i}\}_{i\in I(k)}$ be a partition of $S$ with width $\rho_{k}$ that satisfies Assumption 2.3. Fix $s_{0}\in\partial A$. For each $k\in\mathbb{N}$, let $\tilde{s}^{k}$ denote the approximate streamline corresponding to $\{S^{k}_{i}\}_{i\in I(k)}$ with $\tilde{s}^{k}(0)=s_{0}$. If $J_{AB}$ is globally Lipschitz continuous and there exists $0<C<\infty$ that does not depend on $k\in\mathbb{N}$ such that (3.17) holds, then there exists $C^{\prime}>0$ that depends only on $C$, $\|J_{AB}\|_{\textup{Lip}}$ and $s_{0}$, such that for all $k\in\mathbb{N}$,

[TABLE]

Note that the error bound for the approximate probability current in Theorem 3.7 is computed using the $L^{2}$ norm, while the error bound for the approximate streamline in Theorem 3.9 is computed using the $L^{\infty}$ norm. This is because we are interested in the maximal deviation of the approximate streamlines $\tilde{s}^{k}$ from the true streamline, given that they have a common initial condition.

Proof.

Fix an arbitrary $k\in\mathbb{N}$, and let $\tilde{J}^{k}_{AB}$ be the approximate probability current corresponding to $\{S^{k}_{i}\}_{i\in I(k)}$ that generates the approximate streamline $\tilde{s}^{k}$ via (3.20a) and the initial condition $\tilde{s}^{k}(0)=s_{0}$. It follows from (3.19a) that

[TABLE]

By the triangle inequality and Lipschitz continuity of $J_{AB}$, we have

[TABLE]

for all $0\leq r\leq T(s_{0})$. We bound the second term on the right-hand side of the last inequality: for all $0\leq r\leq T(s_{0})$,

[TABLE]

where the equation follows from (3.13), the first inequality follows from the inequality $\left|x\right|_{\infty}\leq\left|x\right|_{2}$ that is valid for any finite-dimensional vector $x$, and the second inequality follows from (3.16) with $C_{k}$ defined in (3.18). Thus for $0\leq t\leq T(s_{0})$,

[TABLE]

and by the Gronwall-Bellman inequality, we obtain the error bound

[TABLE]

Since there exists $0<C<\infty$ such that (3.17) holds, the proof is complete. ∎

Note that the definition of the approximate streamlines does not permit an interpretation in terms of $Y$. This definition allows for two approximate streamlines that lie in the same cell $S_{i}$ at possibly different times to subsequently enter different cells, e.g. if the streamlines exit $S_{i}$ via different facets. In terms of the Delaunay graph $G=(I,E)$ of the Voronoi tessellation, this means that it is possible for $\tilde{s}_{1}(t_{1})=\tilde{s}_{2}(t_{2})=i$ for possibly different $t_{1},t_{2}$ and $\tilde{s}_{1}(t^{\prime}_{1})\neq\tilde{s}_{2}(t^{\prime}_{2})$ for $t^{\prime}_{1}>t_{1}$, $t^{\prime}_{2}>t_{2}$. This non-uniqueness property prevents an interpretation of the discrete streamlines in $G$ for $Y$ as being analogues of the streamlines in $S$ for $X$, because standard results from the theory of ordinary differential equations guarantee that the streamlines of $X$ are uniquely determined by the probability current and initial condition. Since the approximate streamlines are integral curves of the approximate probability current, this implies in addition that the approximate probability current does not permit an interpretation as a probability current for $Y$ on $G$.

4 Numerical results

In this section we implement our method on a well-known example. We focus on the approximate committor and the approximate probability current only, and study the dependence of the errors of these objects as a function of the tessellation width parameter $\rho$.

Let the underlying diffusion process $X$ be the solution to the Smoluchowski SDE

[TABLE]

where $V:\mathbb{R}^{d}\rightarrow\mathbb{R}$ denotes the potential function, $\beta$ is an inverse temperature parameter, i.e. $\beta=(k_{B}\mathcal{T})^{-1}$ with Boltzmann constant denoted by $k_{B}$ and temperature $\mathcal{T}$, and $\Gamma\in\mathbb{R}^{d\times d}$ is a diagonal matrix with friction coefficients on the diagonal.

4.1 Setup, discretisation, trajectory sampling

In our example, $d=2$, and $V$ is the triple well potential from [19] given by

[TABLE]

on the state space $S=\{x=(x_{1},x_{2})\in\mathbb{R}^{2}\hskip 4.30554pt|\hskip 4.30554pt-2\leq x_{1}\leq 2,\hskip 4.30554pt-1.5\leq x_{2}\leq 2.5\}$, as shown in Figure 1(a). We define the metastable sets to be the sublevel sets of $V$ that contain the minima at $(\pm 1,0)$, i.e. $A:=\{(x_{1},x_{2})\in S\hskip 4.30554pt|\hskip 4.30554ptV(x_{1},x_{2})\leq-3,\hskip 4.30554ptx_{1}\leq 0\}$ and $B:=\{(x_{1},x_{2})\in S\hskip 4.30554pt|\hskip 4.30554ptV(x_{1},x_{2})\leq-3,\hskip 4.30554ptx_{1}\geq 0\}$.

For the discretisation of $S$ we used a mesh of side length $h$, where the values of $h$ were chosen from $\{0.5,0.4,0.25,0.1,0.05\}$. Letting $[N]:=\{1,2,\ldots,N\}$, it follows that the corresponding discrete state spaces $I$ were $\{[64],[100],[256],[1600],[6400]\}$, that any cell in any tessellation has at most four adjacent cells, and that the resulting Voronoi tessellations each have width $h\sqrt{2}$. With regards to Assumption 2.3, for each tessellation, we defined $i\in J$ (respectively $i\in K$) if at least one vertex of the square $S_{i}$ belongs to $A$ (resp. $B$). This results in a representation of the sets $A$ and $B$ by the sets $\cup_{j\in J}S_{j}$ and $\cup_{k\in K}S_{k}$ respectively. An example of such a tessellation and the sets corresponding to $J$ and $K$ is shown in Figure 1(b).

We simulated trajectories of $X$ with $\beta=1.67$ and $\Gamma=\mathbb{I}$ in (4.1) using the Euler-Maruyama method, with time step $\Delta t=0.001$. This resulted in the scalar prefactor $\sqrt{2\Delta t\beta^{-1}}\approx 0.034$ for each standard normal random variable in the Euler-Maruyama method.

For the results below, we chose $\mu$ to be the invariant Boltzmann-Gibbs measure with density proportional to $\exp(-\beta V(\cdot))$ on $S$, and used trapezoidal quadrature to compute the $\mu$-measure of cells.

4.2 Approximate committor

To obtain the discrete committor function described in Section 3.1, we sampled $10^{4}$ trajectories for each cell $S_{i}$ in the partition. For simplicity, we drew the initial condition of each trajectory independently from the uniform distribution on $S_{i}$. Since the uniform distribution on $S_{i}$ and the normalised restriction of the ergodic measure $\mu$ to $S_{i}$ in general do not coincide, this choice incurs an error; however, we shall show that for sufficiently small mesh size $h$ this error does not play a significant role relative to other errors. Each sample trajectory was terminated when the state of the trajectory reached $A$ or $B$. The approximate committor $\tilde{q}$ was computed from the computed trajectories using (3.3).

To obtain the reference values of the committor $q$ for each mesh of side length $h$, we used that $q$ solves the Feynman-Kac or backward Kolmogorov PDE

[TABLE]

where $L$ denotes the infinitesimal generator of (4.1); see [8, Eq. (10)]. We used the finite differences method to solve the PDE on each mesh, for the values of $h$ stated above.

To compute the error $\|q-\tilde{q}\|_{L^{2}(\mu)}$, we computed the difference of the vector versions of $q$ and $\tilde{q}$, weighted by the vector of $\mu$-measures of each cell. We plot the results in Figure 2(a) and observe that the observed values of the $L^{2}(\mu)$ error of the approximate committor decreased linearly with the partition width for the considered values of $\rho$, with an observed slope of approximately 1. Thus, the numerical results for the approximate committor are consistent with the bound in Corollary 3.5.

4.3 Approximate probability current

To compute the approximate probability current, we need to compute for each $i\in I$ the vector $\tilde{J}_{AB,i}\in\mathbb{R}^{d}$, which in turn requires the vector $\hat{\alpha}_{i}\in\mathbb{R}^{\textup{card}({\mathcal{N}(S_{i})})}$ of surface area-normalised net flows of reactive trajectories; see (3.11). To do this, we followed (3.9) and sampled a long trajectory starting from the origin until $10^{5}$ reactive trajectory segments were collected. We then followed the procedure described in the paragraph immediately after Remark 3.6. The duration of this long trajectory was computed as $T:=M\Delta t$, for $M$ the number of Euler-Maruyama steps. The observation time window parameter was set to be the Euler-Maruyama time step, i.e. $s:=\Delta t$. In particular, both limits were approximated by taking fixed values of $s$ and $T$. Using the reactive trajectory segments, we computed $\alpha_{ik}$ for each $i\in I$ and for each $k\in\mathcal{N}(S_{i})$. Since each set $S_{i}$ in the partition is a square of side length $h$, it follows that $\sigma(\partial S_{i}\cap\partial S_{k})=h$ for any adjacent $S_{i}$ and $S_{k}$. These sets of numbers then yielded the $\hat{\alpha}_{ik}$ by (3.10). For each $i\in I$, we solved the system (3.12) for $\tilde{J}_{AB,i}$ by using the fact that $M_{i}=2I_{2}$, where $I_{2}$ is the $2\times 2$ identity matrix; see the paragraph immediately before §3.3. We then constructed the approximate probability current according to (3.13).

To obtain the reference values of the probability current $J_{AB}$ for each mesh of side length $h$, we used the fact that for diffusion processes given by Smoluchowski SDEs, the probability current is proportional to the product of the Boltzmann-Gibbs density with the gradient of the committor

[TABLE]

see [16, Eq. (6)]. For each mesh of side length $h$, we used the finite differences approximation of $q$ to approximate the gradient $\nabla q$ and computed the reference values of the probability current using the equation above.

We plot the $L^{2}(\mu)$ error of the approximate probability current relative to the finite differences solution in Figure 2(b). For the three largest values of $h$, i.e. for $h\in\{0.5,0.4,0.25\}$, the observed values of the $L^{2}(\mu)$ error of the approximate committor decrease linearly with the partition width; the slope of the line is between $1/2$ and $1$. Subsequently, the observed values of the $L^{2}(\mu)$ error increase for $h\in\{0.1,0.05\}$. These observed results are not consistent with the convergence behaviour described by Corollary 3.8.

The increase in $L^{2}(\mu)$ error for $h\in\{0.1,0.05\}$ is a consequence of under-counting of transitions, which in turn is due to an observation time window that is too large for the mesh side length parameter $h$. For $\Delta t=0.001$, recall that $\sqrt{2\Delta t\beta^{-1}}\approx 0.034$. Hence, with more than 95% probability, the magnitude of each random perturbation in the Euler-Maruyama method is less than or equal to $0.1$. As a result, we expect to significantly underestimate the number of crossings of $\partial S_{i}\cap\partial S_{k}$ when each $S_{i}$ has side length approximately equal to or less than $\sqrt{2\Delta t\beta^{-1}}$, i.e. for $h\in\{0.05,0.1\}$.

To check the validity of the preceding explanation for the increase in $L^{2}(\mu)$ errors for $h\in\{0.1,0.05\}$, we sampled for each $h\in\{0.25,0.1,0.05\}$ a single long trajectory with $\Delta t=0.00025$ and computed the corresponding approximate probability current. For $\Delta t=0.00025$, $\sqrt{2\Delta t\beta^{-1}}\approx 0.017$ and thus with more than 95% probability the magnitude of each random perturbation in the Euler-Maruyama method is less than or equal to $0.051$. The results of the additional experiments are plotted in Figure 2(b): for $h=0.1$ and $\Delta t=0.00025$ the result approximately follows the same linear trend of the observed $L^{2}(\mu)$ error for $h\in\{0.5,0.4,0.25\}$ and $\Delta t=0.001$. For $h=0.05$, $\sqrt{2\Delta t\beta^{-1}}\approx h$, and we observe that the $L^{2}(\mu)$ error deviates significantly from this trend. These results support our proposed explanation.

We further analysed the error between the approximate probability current $\tilde{J}_{AB}$ and the probability current $J_{AB}$ obtained from finite differences by computing two error metrics at each point: the scaling ratio $R(x)$ and the direction error $D(x)$,

[TABLE]

where we choose units so that $0\leq D(x)\leq\pi$. For $\Delta t=0.00025$ and for each value of $h\in\{0.05,0.1,0.25\}$, we computed the scaling ratio and direction error for cells $S_{i}$ on which the value of the potential $V$ is less than or equal to 1. We restrict to these cells because the sampled reactive trajectories spend significantly less time in the region $[V>1]:=\{x\in S\ :\ V(x)>1\}$ compared to its complement $[V\leq 1]:=\{x\in S\ :\ V(x)\leq 1\}$. As a result, we do not expect the observed statistics for cells in $[V>1]$ to provide reasonable approximations of their true values in $[V>1]$, compared to the observed statistics for cells in $[V\leq 1]$. In addition, we thresholded the values of the scaling ratio to the interval $0\leq R(x)\leq 2$ to reduce the influence of large outliers and to make the comparison symmetric about the desired scaling ratio of 1. We plot the frequency histograms for the direction error $D$ and scaling ratio $R$ in Figures 3(a) and 3(b).

In Figure 3(a) we observe that the frequency histograms of the direction error $D(x)$ are approximately the same for $h\in\{0.05,0.1,0.25\}$. In contrast, Figure 3(b) shows that the frequency histograms for the scaling ratio $R(x)$ differ in some aspects for varying $h$: as $h$ decreases from $h=0.25$ to $h=0.05$, the frequency histogram of the scaling ratio becomes more sharply peaked, and in addition the observed mean of the frequency histogram decreases away from 1. In addition, for $h=0.25$, the frequency of large outliers is at least twice the frequency for $h=0.1$ and $h=0.05$.

Together, Figures 2(b), 3(a), and 3(b) suggest that, relative to the direction error $D(x)$, the observed distribution of the scaling ratio $R(x)$ of the approximate probability current plays a more important role in explaining the discrepancy between the observed and predicted slopes of the numerical $L^{2}(\mu)$ error of the approximate probability current. To support this claim, we plot the direction error $D(x)$ and the absolute deviation $|R(x)-1|$ of the scaling ratio for the numerical results obtained with the parameters $h=0.25$ and $\Delta t=0.00025$ in Figures 4(a) and 4(b) respectively.

In Figure 4(a) we observe that the direction error is less than 0.5 over much of the region of interest. In contrast, in Figure 4(b) the absolute deviation $|R(x)-1|$ is more than 0.9 over much of the region of interest. These figures suggest that, given enough reactive trajectory segments, our method can identify the direction of the probability current $J_{AB}$ better than the magnitude of the probability current. The errors in the estimation of the probability current tend to be larger over areas that the reactive trajectories visit less often, although there are exceptions. Remarkably, the plot analogous to Figure 4(b) for $\Delta t=0.001$ (not shown) is very similar to the plot for $\Delta t=0.00025$ shown above; this suggests that the observation time window $s$ does not play a significant role in the observed deviations in scaling ratio. These observations suggest that too short a trajectory duration $T$, or equivalently, insufficient sampling of reactive trajectories explains the discrepancy between the observed slope and the predicted slope of the numerical $L^{2}(\mu)$ errors in the approximate probability current.

5 Conclusion

In this paper, we presented a method for approximating key objects from TPT for ergodic diffusion processes, namely the committor, probability current, and streamlines, using a Voronoi tessellation of state space and sample trajectory data. In the case of the approximate committor, only sample trajectory data are needed, and we can interpret the approximate committor as the committor of a non-Markovian jump process on the Delaunay graph $G$. For the approximate probability current and the approximate streamlines, sample reactive trajectory data are needed, and we cannot interpret either object as a discrete analogue on $G$ of the corresponding object from TPT for diffusion processes.

For each approximate object, we proved a bound on the approximation error, while neglecting for simplicity both statistical error due to finite sample sizes and errors due to the numerical method used for generating trajectories of the diffusion process. We used these bounds to prove that the approximate objects converge to the continuous objects from TPT in the continuum limit. We presented numerical results for the approximate committor that are consistent with our theoretical prediction of convergence behaviour. We presented numerical results for the approximate probability current, and explained their deviations from the predicted values in terms of the approximation of the $s\to 0^{+}$ and $T\to\infty$ limits in (3.9) with fixed values of $s$ and $T$.

The approach we have taken here is statistical and data-driven in nature, in the sense that our method relies only on computing observed statistics from sampled reactive trajectories and the Voronoi tessellation. This may be a disadvantage when reactive trajectory data is difficult to obtain. Furthermore, the numerical results in §4.3 showed that statistical error or under-sampling in regions of low probability can lead to slow convergence. On the other hand, the data-driven feature may be advantageous when the task of obtaining trajectory data is easier compared to the task of learning the drift and diffusion coefficients of the SDE.

Our method differs from methods for TPT that involve constructing a discrete state-space Markov process for approximating the underlying diffusion process. This feature comes at the cost of interpretability: the approximate probability current and streamlines cannot be interpreted as discrete TPT objects on the Delaunay graph associated to the Voronoi tessellation. On the other hand, we do not need to ensure that an approximating process obtained from discretisation of a continuous state space possesses the Markov property. In particular, we do not need to identify a suitable lag time, and can therefore avoid the preprocessing steps and heuristic reasoning that are often used to identify the lag time. In addition, the error and convergence analysis of the approximation error is simpler, because the definitions of the approximate objects of our method are not based on preserving certain properties on a graph, but on the continuous objects from TPT for diffusion processes.

We now mention some aspects that should be addressed in future work in order to improve the applicability of our method.

In order to compute the approximate probability current, one needs to compute the $(d-1)$-dimensional volumes of the facets between adjacent Voronoi cells; see (3.10). Computing these volumes requires computing the full Voronoi tessellation, which will be expensive in high dimensions. One possible approach to circumvent this would be to use Voronoi tessellations in which every Voronoi cell is a rigid transformation of one Voronoi cell. However, this would prevent the user from incorporating a priori information about the SDE such as the energy landscape when placing generators for the tessellation. Thus, it would be desirable to modify the approach so that it is not necessary to compute the volume of every facet between adjacent cells. A possible approach would be to avoid a full tessellation of state space, e.g. by using core sets.

In §4.3, we observed that the simple Monte Carlo approach that we used above could accurately estimate the committor probability, but not the probability current, over all regions in state space. One important reason for this is that it is easier to obtain the trajectory data needed for estimating committor probabilities than to obtain reactive trajectory data for estimating net fluxes: for the former, one starts a trajectory from a point or region of choice and waits until the trajectory reaches $A$ or $B$; for the latter, one starts a reactive trajectory segment on $\partial A$ and waits until it reaches $\partial B$ without re-entering $A$. Another important reason is that to construct the approximate probability current, one needs to estimate net fluxes across all facets. The accuracy of the net flux estimates depends strongly on how well the sampled reactive trajectories explore state space: the net flux estimates in regions that are explored less will be less accurate. It would be desirable to modify the approach so that sampling reactive trajectories is easier, and so that the sample reactive trajectories explore state space more uniformly. This could be achieved by using enhanced sampling techniques, for example.

For applications where trajectory data is difficult to obtain, it would be of interest to extend the method so that it could be applied to sparse trajectory data, where ‘sparse’ means that subsequent measurement times are far apart from each other, or where the trajectory data degenerates to a cloud of point data. In full generality, such a task appears intractable. However, under certain assumptions, e.g. the assumption that the points lie on an energy landscape that determines the dynamics, it might be possible to develop such an extension. For example, a recent analysis [9] has shown that applying the so-called square root approximation [11] to generate an approximate transition rate matrix from point cloud data on an energy landscape produces an approximation of the associated infinitesimal generator. It would be of interest to modify the definition of the approximate probability current, in order to make our approach applicable to such data.

Acknowledgments

NC and TC were supported by project grant CH14 of the Einstein Center for Mathematics Berlin (ECMath), and by the DFG Research Center Matheon “Mathematics for key technologies” in Berlin. TC is supported by the German Ministry of Research and Education (BMBF) project grant 3FO18501 (Forschungscampus MODAL). The research of NC and HCL has been partially funded by Deutsche Forschungsgemeinschaft (DFG) - SFB1294/1 - 318763901.

The authors thank Christof Schütte and Gabriel Stoltz for their helpful feedback on early versions of this work.

Appendix A Auxiliary results and proofs

A.1 Approximate committor

The purpose of this section is to prove Proposition 3.1, i.e. that for a given tessellation $\{S_{i}\}_{i\in I}$, $\hat{q}_{i}=\tilde{q}_{i}$ for all $i\in I$. We accomplish this by using the notion of a regular conditional distribution. To define a regular conditional distribution, we first recall some definitions; see [10, Section 8.3].

Definition A.1** (Stochastic kernel).**

Let $(\Omega_{1},\mathcal{A}_{1})$ and $(\Omega_{2},\mathcal{A}_{2})$ be measurable spaces. A map $\kappa:\Omega_{1}\times\mathcal{A}_{2}\rightarrow\left[0,\infty\right]$ is called a stochastic kernel from $\Omega_{1}$ to $\Omega_{2}$ if

• (i)

$\kappa(\cdot,A_{2})$ is $\mathcal{A}_{1}$-measurable for any $A_{2}\in\mathcal{A}_{2}$, and

• (ii)

$\kappa(\omega_{1},\cdot)$ is a $\sigma-$finite probability measure on $(\Omega_{2},\mathcal{A}_{2})$ for any $\omega_{1}\in\Omega_{1}$.

Definition A.2** (Regular conditional distribution).**

Let $(\Omega_{1},\mathcal{A}_{1},\mathbb{P})$ be a probability space, $(E,\mathcal{E})$ be a measurable space, and $Y$ be an $E$-valued random variable on $(\Omega_{1},\mathcal{A}_{1},\mathbb{P})$. Let $\mathcal{F}\subset\mathcal{A}_{1}$ be a $\sigma$-algebra on $\Omega_{1}$. A stochastic kernel $\kappa_{Y,\mathcal{F}}$ from $(\Omega_{1},\mathcal{F})$ to $(E,\mathcal{E})$ is called a regular conditional distribution of $Y$ given $\mathcal{F}$ if

[TABLE]

for $\mathbb{P}-$almost all $\omega\in\Omega_{1}$ and for all $C\in\mathcal{E}$.

If $\mathcal{F}$ is generated by a random variable $X$ defined on $(\Omega_{1},\mathcal{A}_{1},\mathbb{P})$ that takes values in some measurable space $(E^{\prime},\mathcal{E}^{\prime})$ then the stochastic kernel $\kappa_{Y,X}$ from $(E^{\prime},\mathcal{E}^{\prime})$ to $(E,\mathcal{E})$ that is defined by

[TABLE]

is called a regular conditional distribution of $Y$ given $X$.

The existence of the regular conditional distribution of $Y$ given $\mathcal{F}$ for a random variable $Y$ taking values in a Borel space is given in [10, Theorem 8.37].

Lemma A.3** (Committors are regular conditional probabilities).**

The committor function $q$ defined in (3.1) is a regular conditional probability.

Proof.

Set the random variables $X,Y,$ and the set $C$ in (A.2) to be equal $X_{0}$, $X_{\tau_{A\cup B}(X)}$, and $B$, respectively. Then

[TABLE]

where we used the definition (3.1) of the committor function in the second equation. ∎

The following result is reproduced from [10, Theorem 8.38].

Theorem A.4**.**

Let $Y$ be a random variable on $(\Omega_{1},\mathcal{A}_{1},\mathbb{P})$ with values in some set $E$, and equip $E$ with the Borel $\sigma$-algebra $\mathcal{E}$. Let $\mathcal{F}\subset\mathcal{A}_{1}$ be a $\sigma$-algebra and let $\kappa_{Y,\mathcal{F}}$ be a regular conditional distribution of $Y$ given $\mathcal{F}$. Further, let $f:E\rightarrow\mathbb{R}$ be measurable and $\mathbb{E}\left[|f(Y)|\right]<\infty.$ Then

[TABLE]

$\mathbb{P}$-almost surely.

We use Theorem A.4 to prove the following lemma, which expresses the joint probability of a pair of random variables $X$ and $Y$ in terms of the regular conditional probability of $Y$ given $X$.

Lemma A.5**.**

Let $X$ and $Y$ be random variables on $(\Omega_{1},\mathcal{A}_{1},\mathbb{P})$, taking values in the measurable spaces $(E^{\prime},\mathcal{E}^{\prime})$ and $(E,\mathcal{E})$ respectively. Then for any $C\in\mathcal{E}$ and $D\in\mathcal{E}^{\prime}$,

[TABLE]

Proof.

Let $\sigma(X)$ denote the $\sigma$-algebra generated by $X$. Then

[TABLE]

Observe that by Theorem A.4, (A.2), and Definition A.1,

[TABLE]

Therefore,

[TABLE]

where we used (A.3) and (A.4) in the first equation and the change of variables formula in the second equation. This completes the proof. ∎

We now use Lemma A.5 to prove Proposition 3.1, which states that for every $i\in I$, the quantities $\hat{q}_{i}$ and $\tilde{q}_{i}$ defined in (3.4) and (3.2) are equal, under the assumption that $\mathbb{P}\circ(X_{0})^{-1}=\mu$.

Proof of Proposition 3.1.

Let $i\in I$ be arbitrary. Recall from (3.4) that

[TABLE]

The definition (3.2), the construction of $Y$ and Lemma 2.4, and the hypothesis that $\mathbb{P}\circ(X_{0})^{-1}=\mu$ imply that

[TABLE]

Thus, it suffices to show that

[TABLE]

By Lemma A.3, the left-hand side can be rewritten as

[TABLE]

Using that $\mu=\mathbb{P}\circ(X_{0})^{-1}$ and Lemma A.5, we obtain

[TABLE]

yielding the desired conclusion. ∎

A.2 Approximate probability current

In this section, we consider $d\in\mathbb{N}$, $d\geq 2$. By ‘$d$-polytope’, we mean a bounded convex polytope $P$ with $\dim{P}=d$.

Lemma A.6**.**

If $P\subset\mathbb{R}^{d}$ is a $d$-polytope, then $P$ at least $d+1$ facets.

Proof.

Any such polytope admits a decomposition into $d$-simplices. Since every $d$-simplex has $d+1$ facets, the conclusion follows. ∎

Corollary A.7**.**

If $P$ is a $d$-polytope in $\mathbb{R}^{d}$, then $P$ has at least $d$ linearly independent outer normals. In particular, the matrix $N$ whose rows are the outer unit normals to $P$ has full rank.

Proof.

The last statement follows from the first, so it suffices to prove the first statement. By Lemma A.6, every $d$-polytope has at least $d+1$ outer normals. We will prove by contradiction that there exist $d$ linearly independent outer normals. Suppose $P$ has no more than $d-1$ linearly independent outer normals. Then the normals to the facets of $P$ span at most a $(d-1)$-dimensional space, which implies that there exists a hyperplane $H$ in $\mathbb{R}^{d}$ containing all the outer normals of $P$. Let $v$ be normal to $H$, and let $n$ be an arbitrary outer normal associated to to some facet $F$ of $P$. Then $v$ and $n$ are orthogonal, which implies that $v$ is parallel to $F$, and thus that $F$ is unbounded along the direction of $v$. This implies that $P$ is unbounded, which contradicts the boundedness hypothesis on $P$. ∎