Parameter and dimension dependence of convergence rates to stationarity for Reflecting Brownian Motions

TL;DR

This paper derives explicit convergence rates to stationarity for multidimensional reflected Brownian motions, improving existing bounds and providing new dimension-dependent estimates applicable to various rank-based diffusions.

Contribution

It provides the first explicit dimension and parameter-dependent convergence rates for RBMs and rank-based diffusions under Wasserstein distance, without requiring stationarity or reversibility.

Findings

Improved relaxation time bounds from O(d^4(log d)^2) to O((log d)^2) for certain RBMs.

First explicit rates for rank-based diffusions including the Atlas model.

Relaxation time bound of O(d^6(log d)^2) for the standard Atlas model.

Abstract

We obtain rates of convergence to stationarity in L^1-Wasserstein distance for a d-dimensional reflected Brownian motion (RBM) in the nonnegative orthant that are explicit in the dimension and the system parameters. The results are then applied to a class of RBMs considered in Blanchet-Chen (2016) and to rank-based diffusions including the Atlas model. In both cases, we obtain explicit rates and bounds on relaxation times. In the first case we improve the relaxation time estimates of O(d^4(log d)^2) obtained in Blanchet-Chen (2016) to O((log d)^2). In the latter case, we give the first results on explicit parameter and dimension dependent rates under the Wasserstein distance. The proofs do not require an explicit form for the stationary measure or reversibility of the process with respect to this measure, and cover settings where these properties are not available. In the special case of the standard Atlas model, we obtain a bound on the relaxation time of O(d^6(log d)^2).

Equations

$$X(t;\mathbf{x})=\mathbf{x}+\mathbf{D}B(t)+\mu t+RL(t),$$

$$L(0)=0,{\int }_0^t{X}_i(s;\mathbf{x})d{L}_i(s)=0\text{mbox}forallt>0\text{mbox}and1\le i\le d.$$

$${\eta }_i^k(\mathbf{x})=\inf \{t\ge {\eta }^{k-1}(\mathbf{x})+1:X_i(t,\mathbf{x})=0\},{\eta }^k(\mathbf{x})=\sup \{{\eta }_i^k(\mathbf{x}):1\le i\le d\}.$$

$$\mathcal{N}(t;\mathbf{x})=\sup \{k\ge 0:{\eta }^k(\mathbf{x})\le t\}.$$

$$n(R):=\inf \{n\ge 1:\Vert P^{n}1{\Vert }_{\infty }\le 1/2\},$$

$$W_1(\mu ,\nu )=\inf \left\{{\int }_{{\mathbb{R}}_{+}^d\times {\mathbb{R}}_{+}^d}\Vert \mathbf{x}-\mathbf{y}{\Vert }_1\gamma (d\mathbf{x},d\mathbf{y}):\gamma \text{ is a coupling of }\mu \text{ and }\nu \right\},$$

$$t_{rel}(\mathbf{x}):=\inf \{t\ge 0:W_1\left(\mathcal{L}(X(t;\mathbf{x})),\mathcal{L}(\mathbf{X}(\infty ))\right)\le 1/2\}.$$

$$a(\Theta )$$

$$R_1(\Theta ,d)$$

$$C_1(\mathbf{x},\Theta )$$

$$C_2(\mathbf{x},\Theta ,\kappa )$$

$$W_1\left(\mathcal{L}(X(t;\mathbf{x}),\mathcal{L}(\mathbf{X}(\infty ))\right)$$

$$\le C_1(\mathbf{x},\Theta )\left(2e^{-\frac{D_{1}t}{R_1(\Theta ,d)}}+e^{-\frac{t}{16D_2{R}_2(\Theta )}}\right)+C_2(\mathbf{x},\Theta ,D_2)e^{-\frac{t}{8D_2{R}_2(\Theta )}}.$$

$$t_{rel}(\mathbf{x})\le \max \{D_1^{-1}{R}_1(\Theta ,d)\log (8C_1(\mathbf{x},\Theta ))+16D_2{R}_2(\Theta )\log [4(C_1(\mathbf{x},\Theta )+C_2(\mathbf{x},\Theta ,D_2))],$$

$$t_0\left(1+(a(\Theta ){)}^2\log (2d)\right)\}.$$

$$\begin{gathered} W_1\left(\mathcal{L}(X(t;\mathbf{x}),\mathcal{L}(\mathbf{X}(\infty ))\right)\le \mathbb{E}(\Vert X(t;\mathbf{x})-X(t;\mathbf{X}(\infty )){\Vert }_1) \\ \le C_1(\mathbf{x},\Theta )\left(e^{-\frac{D_{1}n(R)t}{R_1(\Theta ,d)}}+e^{-\frac{32(\log 2)D_{1}t}{R_1(\Theta ,d)}}+e^{-\frac{t}{16D_2{R}_2(\Theta )}}\right)+C_2(\mathbf{x},\Theta ,D_2)e^{-\frac{t}{8D_2{R}_2(\Theta )}}. \end{gathered}$$

$$\Vert X(t;\mathbf{x})-X(t;0){\Vert }_1\le 2\Vert \mathbf{x}{\Vert }_1{2}^{-\mathcal{N}(t;\mathbf{x})/n(R)}.$$

$$\Vert X(t;\mathbf{x})-X(t;0){\Vert }_1\le \Vert P^{\mathcal{N}(t;\mathbf{x})}1{\Vert }_{\infty }\Vert \mathbf{x}{\Vert }_1.$$

$$\Vert P^{n}1{\Vert }_{\infty }\le C(R,d)2^{-n/n^{\prime }(R)},n\ge 0.$$

$$\Vert X(t;\mathbf{x})-X(t;0){\Vert }_1\le C(R,d)\Vert \mathbf{x}{\Vert }_1{2}^{-\mathcal{N}(t;\mathbf{x})/n^{\prime }(R)}.$$

$$R_1^{\prime }(\Theta ,d)$$

$$C_1^{\prime }(\mathbf{x},\Theta ,d)$$

$$\mathbb{E}(\Vert X(t;\mathbf{x})-X(t;\mathbf{X}(\infty )){\Vert }_1)\le 2\left(2\Vert \mathbf{x}{\Vert }_1+E_1{d}^2\right)e^{-E_{2}t/\log (2d)}+\left(4\Vert \mathbf{x}{\Vert }_1+E_1{d}^2\right)e^{-E_{4}t/2}+E_3{d}^2{e}^{-E_{4}t}.$$

$$t_{rel}(\mathbf{x})\le \max \{E_2^{-1}\log \left[8\left(2\Vert \mathbf{x}{\Vert }_1+E_1{d}^2\right)\right]\log (2d)+E_4^{-1}\left[2\log \left[8\left(4\Vert \mathbf{x}{\Vert }_1+E_1{d}^2\right)\right]+\log (8E_3{d}^2)\right],$$

$$t_1\max \{\Vert \mathbf{x}{\Vert }_{\infty },\log (2d)\}\}.$$

$$\Vert P^{n}1{\Vert }_{\infty }\le 1^T{P}^{n}1\le d\Vert 1^T{P}^n{\Vert }_{\infty }\le d\kappa (1-\beta {)}^n.$$

$$C(R,d)=\kappa d,n^{\prime }(R)\equiv n^{\prime }:=\frac{\log (2)}{\log (1-\beta {)}^{-1}}+1.$$

$$a(\Theta )$$

$$b(\Theta )$$

$$R_1^{\prime }(\Theta ,d)$$

$$R_2(\Theta )$$

$$C_1^{\prime }(\mathbf{x},\Theta ,d)$$

$$C_2(\mathbf{x},\Theta ,D_2)e^{-\frac{t}{8D_2{R}_2(\Theta )}}$$

$$+a(\Theta ){\left[2d(1+d)\left(\sum _{i,j}(R^{-1}{)}_{ij}^2\right)\left(\sum _{j=1}^d{\sigma }_j^2\right)\right]}^{1/2}{e}^{-\frac{t}{8D_2{R}_2(\Theta )}}$$

$$\le 2\Vert \mathbf{x}{\Vert }_1{e}^{-\frac{t}{16D_2{R}_2(\Theta )}}+\frac{\kappa \sigma }{\beta \delta }{\left[2d(1+d)d{\left(\Vert R^{-1}1{\Vert }_{\infty }\right)}^2\left(d{\sigma }^2\right)\right]}^{1/2}{e}^{-\frac{t}{8D_2{R}_2(\Theta )}}$$

$$\le 2\Vert \mathbf{x}{\Vert }_1{e}^{-\frac{t}{16D_2{R}_2(\Theta )}}+\frac{2{\kappa }^2{\sigma }^2}{{\beta }^2\delta }{d}^2{e}^{-\frac{t}{8D_2{R}_2(\Theta )}}$$

$$\le 2\Vert \mathbf{x}{\Vert }_1{e}^{-\frac{{\beta }^3{\delta }^{2}t}{16D_2{\kappa }^3{\sigma }^4}}+\frac{2{\kappa }^2{\sigma }^2}{{\beta }^2\delta }{d}^2{e}^{-\frac{{\beta }^3{\delta }^{2}t}{8D_2{\kappa }^3{\sigma }^4}}.$$

$$\begin{gathered} \mathbb{E}(\Vert X(t;\mathbf{x})-X(t;\mathbf{X}(\infty )){\Vert }_1)\le C_1^{\prime }(\mathbf{x},\Theta ,d)\left(2e^{-\frac{D_{1}t}{R_1^{\prime }(\Theta ,d)}}+e^{-\frac{t}{16D_2{R}_2(\Theta )}}\right)+C_2(\mathbf{x},\Theta ,D_2)e^{-\frac{t}{8D_2{R}_2(\Theta )}} \\ \le \left(2\Vert \mathbf{x}{\Vert }_1+E_1{d}^2\right)\left(2e^{-E_{2}t/\log (2d)}+e^{-E_{4}t/2}\right)+2\Vert \mathbf{x}{\Vert }_1{e}^{-E_{4}t/2}+E_3{d}^2{e}^{-E_{4}t}. \end{gathered}$$

Full text

60J60,

60H10,

60J55,

60K25,

Reflected Brownian motion,

local time,

coupling,

Wasserstein distance,

relaxation time,

heavy traffic,

Atlas model,

1 Introduction

A $d$-dimensional obliquely reflected Brownian motion with drift in the nonnegative orthant plays a central role in Queuing Theory where it arises as a diffusion limit of scaled queue length processes when the system is in the heavy traffic regime (namely the arrival rate and the service rate are approximately equal) [25, 15, 14, 29, 5]. Such a process is also used to describe the behavior of rank-based diffusions, namely a system of particles whose trajectories are given by Brownian motions with drift, where the drift and diffusion coefficients of a given particle at any given time depend on its relative rank in the system at that time. These models appear frequently in mathematical finance, eg. the Atlas model [23, 18, 2]. There has been extensive work in the study of stability of such reflected Brownian motions (RBM) that gives explicit sufficient conditions for positive recurrence for the RBM and the corresponding queuing systems [16, 9, 1, 8, 27]. In this work, we obtain explicit exponential convergence rates (in Wasserstein distance) to equilibrium for multidimensional reflected Brownian motion (RBM) under a key stability condition identified in [16] (see Assumption (A2)). This assumption is known to be ‘almost necessary’ for stability (see Remark 1 for a precise statement). The convergence rates obtained in this work are explicit (up to some universal constants) in the dimension and system parameters. The system parameters are given by the drift vector, the covariance matrix of the Brownian motion, and the reflection matrix. Stationary distributions of RBM are rarely explicit and the convergence rates of the form obtained in this work provide important information for the construction of numerical schemes that sample from these stationary distributions.

There has been some prior work in this area. Exponential ergodicity was proved in [7] for semimartingale reflecting Brownian motions under the stability condition of [9]. This class includes RBM of the form considered in this work. The paper [7] also established exponential ergodicity of certain reflected diffusions with state dependent drift and diffusion coefficients. The key ingredient in the proof was the construction of a suitable Lyapunov function along with establishing a minorization condition on a sufficiently large compact set (referred to as a ‘small set’). The Lyapunov function provides good control on the exponential moments of the return times to the small set while the minorization condition implies the existence of abstract couplings of two copies of the process (via construction of ‘pseudo-atoms’ as described in Chapter 5 of [22]) which have a positive chance of coalescing inside the small set. Together, they furnish exponential rates of convergence (in a weighted total variation distance). However, due to the somewhat implicit treatment of the process inside the small set, the rates obtained by this method shed little light on how they qualitatively depend on the system parameters or the state dimension. The paper [18] obtained explicit convergence rates for a class of reversible rank-based diffusions with explicit stationary measures using Dirichlet form techniques (which crucially use reversibility). See also the discussion in Section 5.2. The convergence considered in [18] corresponds to that of time averages of bounded functionals of the state process to the corresponding stationary values in probability (see Theorem 1 of [18]), which is considerably weaker than the $L^{1}$-Wasserstein distance considered in the current work. The setting of one-dimensional RBM was considered in [28] where (among other results) an estimate on the spectral gap was provided as a function of the drift and the diffusion coefficient. In a recent work, [3] obtained dimension dependent bounds on rates of Wasserstein convergence for a class of RBM. Under conditions on the drift vector, the covariance matrix of the Brownian motion, and the reflection matrix (see Conditions (BC1)-(BC3) in Section 5), [3] analyzed the behavior of the RBM inside the small set explicitly by considering synchronous couplings (namely, couplings where the RBM starting from different points are driven by the same Brownian motion). Using explicit couplings to obtain better convergence rate estimates is a relatively recent but developing area. See [4, 10, 12, 11] for such results for other classes of diffusions. In this work, we revisit the idea of constructing synchronous couplings for RBM. Under quite general conditions (specifically, the ones introduced in [16] that guarantee the existence of strong solutions and positive recurrence), we construct a suitable Lyapunov function and identify (an appropriate analogue of) a small set that both depend crucially on the process parameters and the state dimension. This, along with a careful treatment of excursions from the small set, enables us to quantify contraction rates in $L^{1}$-distance for synchronous couplings starting from distinct points and thereby obtain rates of Wasserstein convergence that are given explicitly in terms of the system parameters, the state dimension, and some constants (that do not depend on dimension or model parameters). These convergence rates, together with bounds on relaxation times of the RBM that follow from it, are the main results of this work and are given in Theorem 1. In Section 5 we apply these results to the class of RBMs considered in [3] and rank-based diffusions considered in [18]. In the former case, we substantially improve the relaxation time estimates from $O(d^{4}(\log d)^{2})$ obtained in [3] to $O((\log d)^{2})$. In the latter case, we give the first results on explicit parameter and dimension dependent rates under the Wasserstein distance. The proofs do not require an explicit form for the stationary measure or reversibility of the process with respect to this measure, and cover settings where these properties are not available. In the special case of the standard Atlas model [13], we give a bound on the relaxation time of $O(d^{6}(\log d)^{2})$ (see Remark 6).

2 Model, notation and assumptions

Let $B$ be a $d$-dimensional standard Brownian motion and let $\mu\in\mathbb{R}^{d}$ and $\mathbf{D},R\in\mathbb{R}^{d\times d}$. Consider for $\mathbf{x}\in\mathbb{R}^{d}_{+}:=[0,\infty)^{d}$ the $\mathbb{R}^{d}_{+}$-valued continuous stochastic process given by the equation

[TABLE]

where $L$, referred to as the local time process, is a non-decreasing continuous process satisfying

[TABLE]

We will make the following basic assumptions.

**Assumptions: **

• (A1)

The matrix $P:=I-R^{T}$ is substochastic (non-negative entries and row sums bounded above by 1) and transient ($P^{n}\rightarrow 0$ as $n\rightarrow\infty$).

• (A2)

$\mathbf{b}:=-R^{-1}\mu>0$.

• (A3)

The matrix $\Sigma=\mathbf{D}\mathbf{D}^{T}$ is positive definite.

The paper [15] shows that under (A1) there is a unique strong solution to (1) - (2), namely for each $\mathbf{x}\in\mathbb{R}^{d}_{+}$ there is a unique pair of continuous stochastic processes $(X,L)$ satisfying the above equations. This assumption is satisfied by the routing matrix of any single-class open queueing network [15] and consequently diffusion limits of such networks can be characterized by (1) - (2). The collection $\{X(\cdot;\mathbf{x})\}_{\mathbf{x}\in\mathbb{R}^{d}_{+}}$ defines a strong Markov process (see [16]) which we denote as $\operatorname{RBM}(\mu,\Sigma,R)$ and refer to simply as the reflected Brownian motion (RBM). The matrix $R$ describes the reflection mechanism, specifically, the $i$-th column of $R$ gives the direction of reflection on the $i$-th face of the orthant. The conditions on $P$ in particular say that its spectral radius is strictly less than $1$. The matrix $\Sigma=\mathbf{D}\mathbf{D}^{T}$ gives the covariance matrix associated with the diffusion term of (1).

**Notation: ** Although $\mu,\Sigma$ and $R$ depend on the dimension $d$, this dependence is suppressed to avoid cumbersome notation. We will write $\mathbf{b}=-R^{-1}\mu$. The entries of $\mathbf{b}$ will be denoted by $b_{i}$, and the diagonal entries of $\Sigma$ will be denoted by $\sigma_{i}^{2}$, where $1\leq i\leq d$. *All constants appearing in the statements of lemmas and theorems will be universal in that they do not depend on model parameters or the dimension $d$, unless noted otherwise.

Remark 1**.**

Unique strong solutions of the RBM that follow from (A1) imply that any coupling of the driving Brownian motions translate into a coupling of the processes themselves. Throughout this work we will take the family* $\{X(\cdot;\mathbf{x})\}_{\mathbf{x}\in\mathbb{R}^{d}_{+}}$ to be * driven by the same Brownian motion, namely we will consider a synchronous coupling of the processes starting from different initial conditions. Assumption (A2) is the well known ‘stability condition’ which is sufficient for the existence of a stationary measure [16]. The condition is almost necessary for stability in that if $\mathbf{b}_{i}<0$ for some $i$ then the RBM is transient [6]. Assuming (A3) in addition to (A1)-(A2) gives that the strong Markov process $\operatorname{RBM}(\mu,\Sigma,R)$ has a unique stationary probability distribution [16].

3 Main Result

Following [3], define the following stopping times: $\eta^{0}(\mathbf{x})=0$ and

[TABLE]

Define

[TABLE]

Also define the contraction coefficient

[TABLE]

where $\mathbf{1}$ is a $d$-dimensional vector of ones and for $u\in\mathbb{R}^{d}$, $\|u\|_{\infty}:=\sup_{1\leq i\leq d}|u_{i}|$. By Assumption (A1), $n(R)<\infty$. This quantity plays a key role in quantifying the convergence rate to equilibrium.

We now present the main result of this work. Given probability measures $\mu$ and $\nu$ on $\mathbb{R}^{d}_{+}$, a probability measure $\gamma$ on $\mathbb{R}^{d}_{+}\times\mathbb{R}^{d}_{+}$ is said to be a coupling of $\mu$ and $\nu$ if $\gamma(\cdot\times\mathbb{R}^{d}_{+})=\mu(\cdot)$ and $\gamma(\mathbb{R}^{d}_{+}\times\cdot)=\nu(\cdot)$. The $L^{1}$-Wasserstein distance between two probability measures $\mu$ and $\nu$ on $\mathbb{R}^{d}_{+}$ is given by

[TABLE]

where for a vector $z\in\mathbb{R}^{d}$, $\|z\|_{1}=\sum_{i=1}^{d}|z_{i}|$. We will denote the law of a random variable $X$ by $\mathcal{L}(X)$. Recall that from [16], under Assumptions (A1)-(A3), there is a unique stationary distribution of the RBM. Denote by $\mathbf{X}(\infty)$ a random vector sampled from this stationary distribution. Define the relaxation time, $t_{rel}(\mathbf{x})$ for the RBM starting from $\mathbf{x}\in\mathbb{R}^{d}_{+}$ as

[TABLE]

We will abbreviate the parameters of the RBM as $\Theta:=(\mu,\Sigma,R)$. Recall that these parameters are required to satisfy (A1)-(A3). We will quantify rate of convergence to equilibrium in terms of the following functions of $\Theta,d$. Fix $\kappa\in(0,\infty)$. Let

[TABLE]

For any $\mathbf{x}\in\mathbb{R}^{d}_{+}$, define $\|\mathbf{x}\|_{\infty}^{*}:=\sup_{1\leq i\leq d}\sigma_{i}^{-1}\mathbf{x}_{i}$.

Theorem 1**.**

There exist a $t_{0}\in(0,\infty)$ and $D_{1},D_{2}\in(0,\infty)$ such that for every $d\in\mathbb{N}$, $\mathbf{x}\in\mathbb{R}^{d}_{+}$, every parameter choice $\Theta$, and $t\geq t_{0}\left(1+(a(\Theta))^{2}\log(2d)\right)$,

[TABLE]

In particular, the relaxation time satisfies

[TABLE]

Remark 2**.**

The universal constants $t_{0},D_{1}$ and $D_{2}$ will be identified in Sections 6 and 7. Specifically, $t_{0}$ and $D_{1}$ are introduced in Lemma 8 (see (43)) and $D_{2}:=\max\{A_{0},9\}$, where $A_{0}$ is introduced in Lemma 5.

Remark 3**.**

The proof of Theorem 1 (see Remark 8) will show that one can, in fact, obtain a better bound of the form

[TABLE]

This bound leads to a better choice of the universal constants appearing in the exponents of the bound when $n(R)$ is large. As a consequence the bounds on relaxation times and the bounds given in the examples of Section 5 can be slightly improved using the above estimate. However, this improved bound leads to cumbersome expressions in the bounds and relaxation time estimates in Section 5. Moreover, our main goal is to highlight the dependence of the convergence rates on system parameters which is completely captured by Theorem 1. Hence, we do not give details on how the improved bound can be obtained, however see Remark 8 for some additional comments.

An important ingredient in the proof is the following analogue of Lemma 3 from [3] which shows that the synchronous coupling gives an a.s. contraction of the $L^{1}$-distance $\|X(t;\mathbf{x})-X(t;\mathbf{0})\|_{1}$ which can be quantified as follows. The proof is similar to that in [3] and so only a sketch is provided.

Lemma 2** ( see [3]).**

For $\mathbf{x}\in\mathbb{R}^{d}_{+}$ and $t\geq 0$,

[TABLE]

Proof.

The main idea is to associate the substochastic matrix $P$ with a Markov chain on states $\{0,1,\dots,d\}$ absorbed at [math] and show that $\|\mathbf{x}\|_{1}^{-1}\|X(t;\mathbf{x})-X(t;\mathbf{0})\|_{1}$ (assuming $\|\mathbf{x}\|_{1}\neq 0$) is bounded above by the maximum over the initial state $i$ of the probability that, starting from $i$, the Markov chain is not absorbed by time $\mathcal{N}(t;\mathbf{x})$. Using this idea, Lemma 2 in [3] and the proof of Lemma 3 in [3] establish

[TABLE]

The lemma now follows from the definition of $n(R)$ given in (3) above. ∎

Remark 4**.**

The quantity $n(R)^{-1}$ defined in (3) gives an explicit bound on the exponential decay rate of $\|P^{n}\mathbf{1}\|_{\infty}$ with $n$. Note that $n(R)$ possibly depends on the dimension $d$, but the dependence is solely through $R$. Sometimes (as we will see in the first example of Section 5) it is possible to get a better bound in the sense that we can obtain positive constants $C(R,d)$ and $n^{\prime}(R)<n(R)$ such that

[TABLE]

In this case, we can replace the bound in Lemma 2 by

[TABLE]

Remark 5**.**

In cases where we can obtain the better bound (5), the constants $R_{1}(\Theta,d)$ and $C_{1}(\mathbf{x},\Theta)$ appearing in the bounds on Wasserstein distance and relaxation time in Theorem 1 can be replaced by $R_{1}^{\prime}(\Theta,d)$ and $C_{1}^{\prime}(\mathbf{x},\Theta,d)$ respectively, where

[TABLE]

4 Outline of Approach

We now give an outline of our approach.

• (i)

We use a key idea from [3] which shows that, under the synchronous coupling, the $L^{1}$-distance between the two processes $X(\cdot;\mathbf{0})$ and $X(\cdot;\mathbf{x})$ decreases with time. Using this idea, we provide an estimate on the rate of decay of this $L^{1}$-distance in terms of a ‘contraction coefficient’ which quantifies the decay rate of $\|P^{n}\mathbf{1}\|_{\infty}$ with $n$. The precise statement was formulated as Lemma 2 in Section 3.

• (ii)

We use the fact that for any $v>0$ in $\mathbb{R}^{d}$ satisfying $R^{-1}v\leq\mathbf{b}$, one can dominate the process $X(\cdot;\mathbf{x})$ in an appropriate manner by a normally reflected Brownian motion with drift $-v$ in $\mathbb{R}^{d}_{+}$. This process, written as $X^{+}_{v}(\cdot;\mathbf{x})$, is technically simpler to analyze. The idea of dominating an $\operatorname{RBM}(\mu,\Sigma,R)$ by a normally reflected RBM is due to [16]. Next, we choose an appropriate compact set (which plays a role similar to the ‘small set’ in the terminology of [22]) such that one can obtain a tight control over return times to this set (this is done via Lyapunov function techniques in Lemma 5) and, loosely speaking, is such that the $L^{1}$-distance between the synchronously coupled processes $X(\cdot;\mathbf{0})$ and $X(\cdot;\mathbf{x})$ decreases by a constant factor each time the process $X^{+}_{v}(\cdot;\mathbf{x})$ visits this set (this result is formulated in Lemma 7). A crucial ingredient here is the introduction of a suitable weighted norm (see (16)) whose sub-level sets are the appropriate ‘small sets’ with the desired contraction property. The definition of this norm is guided by an analysis of how the maximum process for each coordinate scales with the system parameters. This weighted norm is used to construct the small set and also an appropriate Lyapunov function. These constructions and their properties are studied in Section 6.

• (iii)

In Section 7, we obtain the rate of decay of $\|X(t;\mathbf{x})-X(t;\mathbf{0})\|_{1}$ with time $t$, in terms of the parameter $v$ of the dominating normally reflected RBM, by decomposing the path of $X^{+}_{v}(\cdot;\mathbf{x})$ into excursions from the small set obtained in (ii) and using the estimates from Section 6 for probabilities of certain events associated with these excursions.

• (iv)

Finally in Section 8 we prove our main result, namely Theorem 1, where we obtain explicit parameter and dimension dependent rates of decay in $L^{1}$-Wasserstein distance between the processes $X(\cdot;\mathbf{0})$ and $X(\cdot;\mathbf{x})$ with time $t$ by optimizing the rates derived in (iii) over the parameter $v>0$ of the dominating RBM.

Before proceeding to the proof we apply Theorem 1 in two settings, the first is that of RBM satisfying the assumptions of [3] and the second corresponds to that of rank-based diffusions such as the Atlas model.

5 Examples

We will use Theorem 1 (and Remark 5) to obtain bounds on the rate of convergence to equilibrium in two examples that are discussed in Sections 5.1 and 5.2 below.

5.1 Blanchet-Chen RBM

This refers to the class of RBM under the set of assumptions in [3], namely:

• (BC1)

The matrix $P$ is substochastic and there exist $\kappa>0$ and $\beta\in(0,1)$ not depending on the dimension $d$ such that $\|\mathbf{1}^{T}P^{n}\|_{\infty}\leq\kappa(1-\beta)^{n}$ for all $n\geq 0$.

• (BC2)

There exists $\delta>0$ independent of $d$ such that $R^{-1}\mu<-\delta\mathbf{1}$.

• (BC3)

There exists $\sigma>0$ independent of $d$ such that $\sigma_{i}:=\sqrt{\Sigma_{ii}}$ satisfies $\sigma^{-1}\leq\sigma_{i}\leq\sigma$ for every $1\leq i\leq d$.

Under the above conditions [3] give a polynomial bound of $O(d^{4}(\log d)^{2})$ on the relaxation time of the RBM. As shown in the following theorem, Theorem 1 gives a substantial improvement by establishing a polylogarithmic relaxation time of $O((\log d)^{2})$.

Theorem 3**.**

Under Assumptions (BC1), (BC2) and (BC3), there exist positive constants $E_{1},E_{2},E_{3},E_{4},t_{1}$ such that for any $\mathbf{x}\in\mathbb{R}^{d}_{+},t\geq t_{1}\max\{\|\mathbf{x}\|_{\infty},\log(2d)\}$,

[TABLE]

In particular, the relaxation time satisfies

[TABLE]

Proof.

Observe that

[TABLE]

Thus, the hypothesis of Remark 4 is satisfied with

[TABLE]

Now we will use Theorem 1 in conjunction with Remark 5.

Under Assumptions (BC1), (BC2) and (BC3), we have the following bounds:

[TABLE]

where we have used the observation that under Assumption (BC1), one has $\|R^{-1}\mathbf{1}\|_{\infty}\leq\kappa/\beta$. Next, observe that $\|\mathbf{x}\|_{\infty}^{*}\leq\sigma\|\mathbf{x}\|_{\infty}$. This, along with the bound on $a(\Theta)$ obtained above, implies that for $t\geq\frac{48\kappa\sigma^{2}}{\beta\delta}\|\mathbf{x}\|_{\infty}$, $3(D_{2}a(\Theta)b(\Theta))^{-1}\|\mathbf{x}\|_{\infty}^{*}\leq t/(16D_{2}R_{2}(\Theta))$. Hence, for such $t$,

[TABLE]

Take $E_{1}=\frac{\kappa^{3}\sigma^{2}}{2\beta^{2}\delta},E_{2}=D_{1}\left[n^{\prime}\left(2+\frac{\kappa^{2}\sigma^{2}}{\beta^{2}\delta^{2}}\right)\right]^{-1},E_{3}=\frac{2\kappa^{2}\sigma^{2}}{\beta^{2}\delta},E_{4}=\frac{\beta^{3}\delta^{2}}{8D_{2}\kappa^{3}\sigma^{4}}$. Using the above bounds in Theorem 1 (modified as in Remark 5), for any $\mathbf{x}\in\mathbb{R}^{d}_{+}$, $t\geq\max\left\{t_{0}\left(1+\frac{\kappa^{2}\sigma^{2}}{\beta^{2}\delta^{2}}\log(2d)\right),\frac{48\kappa\sigma^{2}}{\beta\delta}\|\mathbf{x}\|_{\infty}\right\}$,

[TABLE]

This proves the first part of the theorem upon taking $t_{1}=\max\left\{t_{0}\left(1+\frac{\kappa^{2}\sigma^{2}}{\beta^{2}\delta^{2}}\right),\frac{48\kappa\sigma^{2}}{\beta\delta}\right\}$. The bound on the relaxation time follows immediately from the first part. ∎

5.2 Gap process of rank-based diffusions

Rank based diffusions are interacting particle systems where the drift and diffusion coefficient of each particle depends on its rank. Mathematically, they are represented by the SDE:

[TABLE]

for $1\leq i\leq d+1$, where $\{X_{(j)}(t):t\geq 0\}$ denotes the trajectory of the rank $j$ particle as a function of time $t$ ($X_{(1)}(t)\leq\dots\leq X_{(d+1)}(t)\text{ for all }t\geq 0$), $\delta_{j},\sigma_{j}$ denote the drift and diffusion coefficients of the rank $j$ particle, and $W_{i}$, $1\leq i\leq d$, are mutually independent standard one dimensional Brownian motions. We will assume throughout that $\sigma_{i}>0$ for all $1\leq i\leq d+1$. Rank-based diffusions have been proposed and extensively studied as models for problems in finance and economics. A special case is the Atlas model [13] where the minimum particle (i.e. the particle with rank $1$) is a Brownian motion with positive drift and the remaining particles are Brownian motions without drift (i.e. $\delta_{i}=0$ for all $i>1$). The general setting considered in (6) was introduced in [2]. In order to study the long time behavior, it is convenient to consider the gap process $Y=(Y_{1},\dots,Y_{d})$, given by $Y_{i}=X_{(i+1)}-X_{(i)}$ for $1\leq i\leq d$. The process $Y\equiv Y(t;\mathbf{y})$ is a RBM in $\mathbb{R}_{+}^{d}$ given as

[TABLE]

where $\mathbf{y}$ is the initial gap sequence, $B$ is a standard $d$-dimensional Brownian motion, $\mu_{i}=\delta_{i+1}-\delta_{i}$ for $1\leq i\leq d$, $\mathbf{D}\in\mathbb{R}^{d\times d}$, $L$ is the local time process associated with $Y$ and $R$ satisfies Assumption (A1). The covariance matrix $\Sigma=\mathbf{D}\mathbf{D}^{T}$ has entries $\Sigma_{ii}=\sigma_{i}^{2}+\sigma_{i+1}^{2}$ when $1\leq i\leq d$, $\Sigma_{i(i-1)}=-\sigma_{i}^{2}$ for $2\leq i\leq d$, $\Sigma_{i(i+1)}=-\sigma_{i+1}^{2}$ for $1\leq i\leq d-1$ and $\Sigma_{ij}=0$ otherwise. In particular, (A3) is satisfied, namely $\Sigma$ is positive definite. Moreover, $R$ is given explicitly as $R=I-P^{T}$, where $P$ is the substochastic matrix given by $P_{i(i+1)}=P_{i(i-1)}=1/2$ for all $2\leq i\leq d-1$, $P_{12}=P_{d(d-1)}=1/2$ and $P_{ij}=0$ if $|i-j|\geq 2$. From [16] the process is positive recurrent and has a unique stationary distribution if Assumption (A2) is satisfied, namely $\mathbf{b}=-R^{-1}\mu>0$, which is same as the following condition.

[TABLE]

In the special case where

[TABLE]

the stationary distribution is explicit and takes the form $\mathcal{L}(\mathbf{Y}(\infty))=\otimes_{k=1}^{d}\operatorname{Exp}(2b_{k}\left(\sigma_{k}^{2}+\sigma_{k+1}^{2}\right)^{-1})$ (see Section 5 of [19]). For the general case (i.e. $\sigma_{i}$ are strictly positive and (7) is satisfied) explicit formulas for stationary distribution are not available, however from [7], the law of $Y(t;\mathbf{y})$ converges to the unique stationary distribution in (weighted) total variation distance at an exponential rate. As noted previously, this result does not provide information on parameter or dimension dependence of the rate of convergence. The paper [18] provides explicit rate of convergence to stationarity, that shows a clear parameter dependence, under the stability condition (7) and the assumption that $\sigma_{i}=1$ for all $1\leq i\leq d$. In this case the stationary measure takes an explicit form and the process is reversible with respect to the stationary measure. The proofs in [18], which are based on Dirichlet form techniques, crucially make use of these properties. The explicit representation of the stationary measure is available only under the skew-symmetry condition (see [17]) guaranteed by (8) and the reversibility of the process with respect to this measure is not available if the $\sigma_{i}$ are not all equal. The convergence considered in [18] corresponds to that of time averages of bounded functionals of the state process to the corresponding stationary values in probability (see Theorem 1 of [18]), which is considerably weaker than the $L^{1}$-Wasserstein distance or total variation convergence.

From Theorem 1 we have the following bound on the rate of $L^{1}$-Wasserstein convergence of the gap process to $\mathbf{Y}(\infty)$. Note that we do not require reversibility or an explicit expression for the stationary measure.

Two key quantities appearing in the rate of convergence are

[TABLE]

where $b_{i}$ are defined in (7) and $\sigma_{i}$ is the standard deviation of the rank $i$ particle (see (6)).

Theorem 4**.**

There exist positive constants $F_{1},F_{2},F_{3},F_{4},t_{2}$ such that for any $\mathbf{y}\in\mathbb{R}^{d}_{+}$ and any $t\geq t_{2}\max\{\sigma^{2}a^{*}\|\mathbf{y}\|_{\infty},1+\sigma^{2}a^{*2}\log(2d)\}$,

[TABLE]

In particular, the relaxation time satisfies

[TABLE]

Proof.

Direct calculation shows that $R^{-1}$ takes the form

[TABLE]

Therefore,

[TABLE]

Using (10) and recalling (9), we obtain

[TABLE]

To compute $R_{1}(\Theta,d)$, we need to estimate $n(R)$. To do this, let $\{S^{*}_{n}\}_{n\geq 0}$ denote a simple, symmetric random walk on $\mathbb{Z}$ starting from $S^{*}_{0}\in\{1,2,\dots,d\}$ and absorbed when it hits [math] or $d+1$. Then for any $n\geq 0$,

[TABLE]

For $j\in\{0,1,\dots,d+1\}$, define $\tau^{S^{*}}_{j}:=\inf\{k\geq 0:S^{*}_{k}=j\}$. By Chapter 10, Example 10.17 of [21], for any $i\in\{1,\ldots,d\}$, $\mathbb{E}\left(\tau^{S^{*}}_{d+1}\wedge\tau^{S^{*}}_{0}\mid S^{*}_{0}=i\right)=i(d+1-i)$. Using this observation and Markov’s inequality, for all $d\in\mathbb{N}$ and $i\in\{1,\ldots,d\}$,

[TABLE]

and consequently,

[TABLE]

Using (11) and (12), we obtain

[TABLE]

For $t\geq 48\sigma^{2}a^{*}\|\mathbf{y}\|_{\infty}$, using the bound on $a(\Theta)$ obtained in (11), and noting $b(\Theta)\geq 1$ and $\|\mathbf{y}\|_{\infty}^{*}\leq\sigma\|\mathbf{y}\|_{\infty}$,

[TABLE]

Moreover, from the explicit form of $R^{-1}$, $\sum_{i,j}(R^{-1})^{2}_{ij}\leq 4d^{4}.$ Using the above two bounds along with (11), for $t\geq 6a^{*}(d+1)^{2}\|\mathbf{y}\|_{\infty}$,

[TABLE]

Take $F_{1}=1,F_{2}=D_{1}/2,F_{3}=4,F_{4}=(2D_{2})^{-1}$. Using the bounds obtained in (11), (13) and (14) in Theorem 1, for any $\mathbf{y}\in\mathbb{R}^{d}_{+},t\geq\max\left\{t_{0}\left(1+\sigma^{2}a^{*2}\log(2d)\right),48\sigma^{2}a^{*}\|\mathbf{y}\|_{\infty}\right\}$,

[TABLE]

This proves the first part of the theorem upon taking $t_{2}=\max\left\{t_{0},48\right\}$. The bound on the relaxation time follows from the first part. ∎

Remark 6**.**

The standard Atlas model [13] is a special case of (6) with $\delta_{1}=1$, $\delta_{i}=0$ for all $i\geq 2$ and $\sigma_{i}=1$ for all $i$. For this model, using (7), for any $k\geq 1$,

[TABLE]

and

[TABLE]

Using these in Theorem 4, we obtain positive constants $G_{1},G_{2},G_{3},G_{4},t_{3}$ such that for any $\mathbf{y}\in\mathbb{R}^{d}_{+}$ and any $t\geq t_{3}\{d^{2}\|\mathbf{y}\|_{\infty},1+d^{2}\log(2d)\}$,

[TABLE]

In particular, the relaxation time for the standard Atlas model is $O(d^{6}(\log d)^{2})$ as $d\rightarrow\infty$.

6 Bounding processes, small sets and return times

Fix a vector $v>0$ satisfying $R^{-1}v\leq\mathbf{b}$ and consider the collection $\{X_{v}^{+}(\cdot;\mathbf{x})\}_{\mathbf{x}\in\mathbb{R}^{d}_{+}}=\operatorname{RBM}(-v,\Sigma,I)$, where $I$ is the identity matrix, given as

[TABLE]

where $B$ is the same Brownian motion as used in the synchronous coupling of $\{X(\cdot;\mathbf{x})\}_{\mathbf{x}\in\mathbb{R}^{d}_{+}}$, and $L^{+}$ is the local time process associated with $X_{v}^{+}$. Observe that $X_{v}^{+}(\cdot;\mathbf{x})$ can be written as

[TABLE]

where $L^{*}(t)=R^{-1}L^{+}(t)+(\mathbf{b}-R^{-1}v)t$ is a non-decreasing process. By minimality of the local time process (see [24, Appendix]), $L^{*}(t)\geq L(t)$ for all $t\geq 0$ implying $R^{-1}X(t;\mathbf{x})\leq R^{-1}X_{v}^{+}(t;\mathbf{x})$ for every $t\geq 0$. Since in this section $v$ will be fixed, we abbreviate $X_{v}^{+}(\cdot;\mathbf{x})$ as $X^{+}(\cdot;\mathbf{x})$. An optimal choice of $v$ will be made later in Section 8. We will hereby refer to $X^{+}(\cdot;\mathbf{x})$ as the bounding process.

We now introduce an appropriate compact set that depends on system parameters and for which one can obtain useful bounds on exponential moments of return times to the set. In order to motivate the choice of the set consider a one dimensional Brownian motion $W_{a,b}(t)=bW(t)-at$ with variance $b^{2}$ and drift $-a$ (here $W$ is a standard one dimensional Brownian motion). Standard techniques using scale functions (see [26, V.46]) show that for any $a>0$ and $b\in\mathbb{R}$, $ab^{-2}\sup_{t<\infty}W_{a,b}(t)$ has an Exponential distribution with mean $1/2$. This result says that the maximum of the $i$-th co-ordinate of $X^{+}(\cdot,\mathbf{x})$ scales like $v_{i}\sigma_{i}^{-2}$. This scaling property suggests considering the following weighted supremum norm :

[TABLE]

This weighted norm will play a central role in our analysis. Also define

[TABLE]

Note that $\phi(v)\geq 2d$. For $A>0$, consider the compact set

[TABLE]

and define the following stopping time for the process $X^{+}(\cdot;\mathbf{x})$:

[TABLE]

The following lemma gives bounds on the exponential moments of the hitting time of the compact set $K_{A}$, namely $\tau_{A}^{+}(\mathbf{x})$.

Lemma 5**.**

There exists $A_{0}>0$ such that for any $A\geq A_{0}$ and any $\mathbf{x}\in\mathbb{R}^{d}_{+}$,

[TABLE]

where $\Lambda(v):=\inf_{1\leq i\leq d}\frac{v_{i}^{2}}{\sigma_{i}^{2}}$.

Proof.

Fix $A>0$ and without loss of generality assume that $\|\mathbf{x}\|_{\infty,v}>A\log\phi(v)$. Consider the ‘Lyapunov function’

[TABLE]

where $g$ is any non-negative, non-decreasing $C^{2}$ function defined on $\mathbb{R}_{+}$ such that $g^{\prime}(0)=0$, $g(u)\leq u,g^{\prime}(u)\leq 2,g^{\prime\prime}(u)\leq 9$ for all $u\geq 0$ and $g(u)=u$ for all $u\geq\log 2$. An example of such a function is $g(u)=(\log 2)h(u/\log 2)\mathbf{1}_{[u\leq\log 2]}+u\mathbf{1}_{[u>\log 2]}$ where $h(u)=u^{4}-3u^{3}+3u^{2}.$ The definition of the Lyapunov function is motivated by a similar function introduced in [3]. The main difference is that here different coordinates are weighted differently depending on system parameters. We will prove that for sufficiently large $A$,

[TABLE]

where $\nabla$ denotes the gradient and $\nabla^{2}$ denotes the Hessian. By Itô’s formula, this will imply that $M(t):=\exp\left(V(X^{+}(t;\mathbf{x}))+\frac{\Lambda(v)}{2A}t\right)$ is a positive supermartingale and therefore, by the optional sampling theorem, for such $A$,

[TABLE]

Since $\|\mathbf{x}\|_{\infty,v}>A\log\phi(v)>A\log d$, we have

[TABLE]

Combining the two displays we have that for $A$ that satisfy (19)

[TABLE]

Thus in order to prove the lemma, it suffices to establish (19) for sufficiently large $A$. Let $w_{i}(\mathbf{y},A)=\frac{e^{g(2A^{-1}v_{i}\sigma_{i}^{-2}y_{i})}}{\sum_{k=1}^{d}e^{g(2A^{-1}v_{k}\sigma_{k}^{-2}y_{k})}}$. By similar calculations as in the proof of Lemma 4 of [3], it follows that

[TABLE]

using $g^{\prime\prime}(u)\leq 9$ and $g^{\prime}(u)\leq 2$ for all $u\geq 0$. Moreover,

[TABLE]

where we have used $\Sigma_{ij}\leq\sigma_{i}\sigma_{j}$ in the first inequality on the second line and the Cauchy-Schwarz inequality, the fact that $g^{\prime}(u)\leq 2$ for all $u\geq 0$, and the fact that $\sum_{j=1}^{d}w_{i}(\mathbf{y},A)=1$ in the last inequality. From the above bounds, we obtain

[TABLE]

Using the definition of $w_{i}$ and the monotonicity of $g$,

[TABLE]

Next, note for any $\|\mathbf{y}\|_{\infty,v}>A\log\phi(v)$, there is $1\leq j\leq d$ such that $v_{j}\sigma_{j}^{-2}y_{j}>A\log\left(\frac{2\sum_{i=1}^{d}v_{i}^{2}\sigma_{i}^{-2}}{\inf_{1\leq i\leq d}v_{i}^{2}\sigma_{i}^{-2}}\right)$. Hence, since $\frac{2\sum_{i=1}^{d}v_{i}^{2}\sigma_{i}^{-2}}{\inf_{1\leq i\leq d}v_{i}^{2}\sigma_{i}^{-2}}\geq 2d\geq 2$, we obtain,

[TABLE]

From (21) and (6),

[TABLE]

Hence, for any $A\geq 68$, we obtain

[TABLE]

proving (19) and hence the lemma holds with $A_{0}=68$. ∎

The next lemma gives an estimate of the running maximum of a reflected Brownian motion with drift.

Lemma 6**.**

Let $X_{t}=x+\sigma^{\prime}B_{t}-\mu^{\prime}t-\min\{\inf_{s\leq t}\left(x+\sigma^{\prime}B_{s}-\mu^{\prime}s\right),0\}$, where $x\geq 0$, $\sigma^{\prime},\mu^{\prime}>0$ and $B$ is a one dimensional standard Brownian motion. Then for any $A,T>0$ and any $x\in[0,A/2]$,

[TABLE]

Proof.

Fix $A,T>0$ and $x\in[0,A/2]$. We define the following stopping times: $\tau_{0}=0$, and for $k\geq 0$,

[TABLE]

Let $\mathcal{N}:=\inf\{k\geq 0:\sup_{t\in[\tau_{2k},\tau_{2k+1}]}X_{t}\geq A\}$. By the strong Markov property, $\{\tau_{2k+1}-\tau_{2k}\}_{k\geq 1}$ are i.i.d, each being distributed as the hitting time of the level $-A/2$ by the process $\sigma^{\prime}B_{t}-\mu^{\prime}t$. By [20, Exercise 5.10], for any $\alpha>0,k\geq 1$,

[TABLE]

Thus, for any $n\geq 0$,

[TABLE]

Optimizing the above bound in $\alpha$ yields the following bound

[TABLE]

Moreover, recalling that the scale function for the process $t\mapsto\sigma^{\prime}B_{t}-\mu^{\prime}t$ is given by $s(z)=e^{2\mu^{\prime}z/\sigma^{\prime 2}}$,

[TABLE]

and hence, for $n\geq 1$,

[TABLE]

From (23) and (24), for any $n\in\mathbb{N}$,

[TABLE]

The result follows on taking $n=\lfloor 4\mu^{\prime}TA^{-1}\rfloor+1$ in (25). ∎

Recall the quantities $\Lambda(v)$ defined in the statement of Lemma 5 and $\phi(v)$ defined in (17). Define

[TABLE]

The next lemma shows that for any $C_{0}\in(0,\infty)$ there are positive constants $C_{1},C_{2}$ such that whenever $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{d}_{+}$ satisfy $\|\mathbf{x}_{1}\|_{\infty,v}\leq C_{0}M(v)$ and $R^{-1}\mathbf{x}_{2}\leq R^{-1}\mathbf{x}_{1}$, with (uniform) positive probability, all the coordinates of $X(\cdot;\mathbf{x}_{2})$ hit zero by time $C_{2}T(v)$ and the weighted supremum norm $\|\cdot\|_{\infty,v}$ of $X^{+}(\cdot;\mathbf{x}_{1})$ is bounded by $C_{1}M(v)$ over the time interval $[0,C_{2}T(v)]$.

Lemma 7**.**

For any $C_{0}>0$, there exists $C_{1}>C_{0}$ and $C_{2}>0$ such that for any $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{d}_{+}$ satisfying $\|\mathbf{x}_{1}\|_{\infty,v}\leq C_{0}M(v)$ and $R^{-1}\mathbf{x}_{2}\leq R^{-1}\mathbf{x}_{1}$,

[TABLE]

Proof.

Let $\{e_{i}\}_{1\leq i\leq d}$ denote the unit coordinate vectors in $\mathbb{R}^{d}$ and let $S_{i}:=\{R^{-1}\mathbf{y}:\mathbf{y}\geq 0,\ y_{i}=0\}$ for $1\leq i\leq d$. Let $U(t;\mathbf{x})=R^{-1}\mathbf{x}+R^{-1}\mathbf{D}B(t)+R^{-1}\mu t$. We first claim that for any $1\leq i\leq d$ and $T>0$,

[TABLE]

To see this, suppose $R^{-1}X(t;\mathbf{x})\notin S_{i}$ for all $0\leq t\leq T$. Then $X(t;\mathbf{x})$ is strictly positive over $[0,T]$. Since $(R^{-1}X)_{i}(t;\mathbf{x})=U_{i}(t,\mathbf{x})+L_{i}(t)$, we have from (2) that $(R^{-1}X)_{i}(t;\mathbf{x})=U_{i}(t,\mathbf{x})$ for all $0\leq t\leq T$. For any $\mathbf{y}\in S\setminus S_{i}$, there exists $\mathbf{z}\in\mathbb{R}^{d}_{+}$ with $z_{i}>0$ such that $\mathbf{y}=R^{-1}\mathbf{z}$. Hence,

[TABLE]

as $(R^{-1})_{ii}=(I+P^{T}+(P^{T})^{2}+\dots)_{ii}\geq 1$. Therefore, $U_{i}(t;\mathbf{x})=(R^{-1}X)_{i}(t;\mathbf{x})>0$ for all $0\leq t\leq T$. This proves (26).

Note that $U(\cdot;\mathbf{x})$ is a Brownian motion with drift in $\mathbb{R}^{d}$ with covariance matrix $R^{-1}\Sigma(R^{-1})^{T}$ and drift vector $-b$. Write $\hat{\sigma}_{i}^{2}:=(R^{-1}\Sigma(R^{-1})^{T})_{ii}$ for the variance of the $i$-th coordinate process $U_{i}$ of $U$. Define for each $1\leq i\leq d$, $\tau^{U}_{i}(\mathbf{x})=\inf\{t\geq 0:U_{i}(t;\mathbf{x})=0\}.$ Also define the vector $\hat{w}$ given by $\hat{w}_{i}=\sigma_{i}^{2}v_{i}^{-1}$ for $1\leq i\leq d$.

For any $i$, recalling $R^{-1}v\leq\mathbf{b}$, note that

[TABLE]

Moreover, using $\Sigma_{jk}\leq\sigma_{j}\sigma_{k}$ for all $1\leq j,k\leq d$,

[TABLE]

From the above two bounds, we conclude from the definition of $T(v)$ that for any $i$,

[TABLE]

Fix $C^{\prime}>0$ and take $\mathbf{y}\in\mathbb{R}^{d}_{+}$ satisfying $R^{-1}\mathbf{y}\leq C^{\prime}(R^{-1}\hat{w})M(v)$. Using (27) and writing $N(0,1)$ for a standard normal random variable, we obtain that for any $C^{\prime\prime}>\max\{2C^{\prime},1\}$,

[TABLE]

where on the last line we have used (27) in the first inequality and $\phi(v)\geq 2d$ in the second inequality.

Recall the upper bounding process $X^{+}$ from (15). Note that the $i$-th coordinate process $X_{i}^{+}$ is a one dimensional reflected Brownian motion with variance $\sigma_{i}^{2}$ and drift $-v_{i}$. Now let $C_{0}>0$ be arbitrary and consider any $C^{\prime}>\max\{2C_{0},1\}$, any $C^{\prime\prime}>\max\{C^{\prime},2\}$. Then, from Lemma 6, for any $\mathbf{x}\in\mathbb{R}^{d}_{+}$ satisfying $\|\mathbf{x}\|_{\infty,v}\leq C_{0}M(v)$,

[TABLE]

where we have used $T(v)\geq(\sigma_{i}^{2}/v_{i}^{2})M(v)\geq(\sigma_{i}^{2}/v_{i}^{2})\log\phi(v)$ in the last step. From the definition of $\phi(v)$ and $T(v)$ respectively,

[TABLE]

Using these observations in (29), we obtain

[TABLE]

where once more we have used $\phi(v)\geq 2d$ to obtain the last bound.

Note that for any $C^{\prime},C^{\prime\prime}>0$, for any $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{d}_{+}$ satisfying $\|\mathbf{x}_{1}\|_{\infty,v}\leq C_{0}M(v)$ and $R^{-1}\mathbf{x}_{2}\leq R^{-1}\mathbf{x}_{1}$,

[TABLE]

Note that for $\mathbf{z}_{1},\mathbf{z}_{2}\in\mathbb{R}_{+}^{d}$, if $\|\mathbf{z}_{1}\|_{\infty,v}\leq CM(v)$ for some $C>0$, and $R^{-1}\mathbf{z}_{2}\leq R^{-1}\mathbf{z}_{1}$, then $R^{-1}\mathbf{z}_{2}\leq C(R^{-1}\hat{w})M(v)$. Using (26) and (28), choosing any $C^{\prime\prime}>1+\max\{2C^{\prime},1\}$, we obtain by the Markov property applied at time $1$,

[TABLE]

Using the estimates (30) and (32) in (31), we obtain for $C^{\prime}>\max\{2C_{0},8\}$, $C^{\prime\prime}=2+\max\{2C^{\prime},33\}$, and any $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{d}_{+}$ satisfying $\|\mathbf{x}_{1}\|_{\infty,v}\leq C_{0}M(v)$ and $R^{-1}\mathbf{x}_{2}\leq R^{-1}\mathbf{x}_{1}$,

[TABLE]

The lemma follows on taking $C_{1}=C^{\prime},C_{2}=C^{\prime\prime}$. ∎

Remark 7**.**

Recall the quantity $A_{0}$ from Lemma 5 and consider $C_{0}\geq A_{0}$. Let $C_{1}$ be as in Lemma 7 associated with this $C_{0}$. Then the set $S:=\{\mathbf{y}:\|\mathbf{y}\|_{\infty,v}\leq C_{1}M(v)\}$ plays a role similar to that of a ‘small set’ in the theory developed in [22], in the following sense. For any $\mathbf{x}\geq 0$, (i) by Lemma 5, we have tight control over return times of the bounding process $X^{+}(\cdot;\mathbf{x})$ to the set $S^{\prime}=\{\mathbf{y}:\|\mathbf{y}\|_{\infty,v}\leq C_{0}M(v)\}\subset S$, and (ii) by Lemma 7, given that the bounding process $X^{+}(t;\mathbf{x})$ lies in $S^{\prime}$ for some $t\geq 0$, then with probability at least a half, all the co-ordinates of $X(\cdot;\mathbf{x})$ hit zero at least once in the time interval $[t,t+T(v)]$ without $X^{+}(\cdot;\mathbf{x})$ leaving $S$. This, in view of Lemma 2, says that $\|X(\cdot;\mathbf{x})-X(\cdot;\mathbf{0})\|_{1}$ is reduced by a factor $2^{-1/n(R)}$ over this time interval.

7 Excursions from the small set

In the following lemma, we combine the estimates from Sections 3 and 6 to decompose the path of $X^{+}(\cdot;\mathbf{x})$ into excursions from the small set (described in Remark 7) and quantify the rate of decay of $\|X(t;\mathbf{x})-X(t;\mathbf{0})\|_{1}$ as $t$ increases.

Lemma 8**.**

For any $A\geq A_{0}$, where $A_{0}$ is the constant appearing in Lemma 5, there exist positive constants $t_{0},D_{1}$ such that for any $\mathbf{x}\in\mathbb{R}^{d}_{+}$, any $v>0$ with $R^{-1}v\leq\mathbf{b}$, and any $t\geq t_{0}T(v)$,

[TABLE]

Proof.

Fix $A\geq A_{0}$ and consider constants $C_{1},C_{2}$ from Lemma 7 that are associated with $C_{0}=A$. We also consider the following stopping times. Let $\tau_{0}=\inf\{t\geq 0:\|X^{+}(t;\mathbf{x})\|_{\infty,v}\leq C_{1}M(v)\}$. For $k\geq 0$, having defined the stopping times $\tau_{0},\dots,\tau_{2k}$, define

[TABLE]

Define $\mathbf{N}_{t}:=\inf\{k\geq 0:\tau_{2k}\leq t\}.$ For $k\geq 0$, define the event

[TABLE]

On the event $E_{k}$, all the coordinates of $X(\cdot;\mathbf{x})$ hit zero in the time interval $[\eta^{\mathcal{N}(\tau_{2k};\mathbf{x})}+1,\tau_{2k+2}]$ as it contains the interval $[\tau_{2k+1}+1,\tau_{2k+2}]$. Consequently, $\mathcal{N}(\tau_{2k+2};\mathbf{x})-\mathcal{N}(\tau_{2k};\mathbf{x})\geq 1$. Thus, for any $k\geq 0$,

[TABLE]

Hence,

[TABLE]

Let $\mathcal{F}_{t}:=\sigma\{B(s):0\leq s\leq t\}$ be the filtration generated by the Brownian motion. For $k\geq 0$, let

[TABLE]

Then $(M_{n},\mathcal{F}_{\tau_{2n}})_{n\geq 1}$ is a martingale with increments bounded by $1$. By Lemma 7, for every $k\geq 0$, $\mathbb{E}\left(\mathbf{1}_{E_{k}}\mid\mathcal{F}_{\tau_{2k}}\right)\geq 1/2.$ Thus, for any $\delta^{\prime}>0$, using the Azuma-Hoeffding inequality with $t\geq 4T(v)/\delta^{\prime}$,

[TABLE]

Note that for any $\delta^{\prime}\in(0,C_{2}^{-1}/2]$,

[TABLE]

where the last inequality follows because $\tau_{2k+2}-\tau_{2k+1}\leq C_{2}T(v)$ and hence, as $\delta^{\prime}\in(0,C_{2}^{-1}/2]$,

[TABLE]

Since $C_{0}=A$ and $M(v)>\log\phi(v)$, for any $k\geq 0$, conditionally on $\mathcal{F}_{\tau_{2k}}$, $\tau_{2k+1}-\tau_{2k}$ is stochastically dominated by $\tau_{A}^{+}(X^{+}(\tau_{2k};\mathbf{x}))$ where $\tau_{A}^{+}(\cdot)$ is defined in (18).

For any $n\geq 0$ $s>0$, using Lemma 5, we obtain

[TABLE]

where we have used $\|X^{+}(\tau_{2n};\mathbf{x})\|_{\infty,v}\leq C_{1}M(v)$ by definition of $\tau_{2n}$, and we take $\sum_{k=0}^{n-1}(\tau_{2k+1}-\tau_{2k})=0$ when $n=0$. Iteratively using the same argument, we obtain

[TABLE]

From (37), for any positive $\delta^{\prime}\leq\min\{C_{2}^{-1}/2,(64C_{1})^{-1}\}$ and $t\geq 3T(v)/\delta^{\prime}$, taking $n=\lfloor\frac{\delta^{\prime}t}{T(v)}\rfloor$ and $s=t/4$,

[TABLE]

where we have used $T(v)=M(v)/\Lambda(v)$ in the equality above. Moreover, as $\tau_{0}\leq\tau^{+}_{A}(\mathbf{x})$, by Lemma 5,

[TABLE]

Using (36), (38) and (39) in (35), for any $\delta^{\prime}$ and $t$ as above,

[TABLE]

From (34) and (40), for positive $\delta^{\prime}\leq\min\{C_{2}^{-1}/2,(64C_{1})^{-1}\}$ and $t\geq 4T(v)/\delta^{\prime}$,

[TABLE]

Now, using Lemma 2, (33) and (41), for positive $\delta^{\prime}\leq\min\{C_{2}^{-1}/2,(64C_{1})^{-1}\}$ and $t\geq 4T(v)/\delta^{\prime}$,

[TABLE]

This proves the lemma with

[TABLE]

∎

8 Main result: optimizing over $v$

In this section, we state and prove our main theorem. This will involve optimizing the bound obtained in Lemma 8 over all possible choices of $v$ along with making an appropriate choice of $A$.

Proof of Theorem 1. Fix any $A\geq A_{0}$ whose value will be appropriately chosen later. Recall from Lemma 8 that the quantities $\frac{D_{1}}{n(R)T(v)}$ and $\frac{\Lambda(v)}{A}$ govern the rate of decay of $\mathbb{E}(\|X(t;\mathbf{x})-X(t;\mathbf{0})\|_{1})$ for any $\mathbf{x}\in\mathbb{R}^{d}_{+}$. To obtain the result in the theorem, we first obtain a value $v(\Theta)$ of $v$ which simultaneously maximizes $\frac{D_{1}}{n(R)T(v)}$ and $\frac{\Lambda(v)}{A}$ over all vectors $v>0$ satisfying $R^{-1}v\leq\mathbf{b}$. For any such $v$, define the vector $s(v)>0$ by

[TABLE]

Then, from the definition of $\Lambda$ we see that $\Lambda(s(v))=\Lambda(v)$. Moreover, $\phi(s(v))=2d\leq\phi(v)$ from definition of $\phi$. Therefore,

[TABLE]

Thus, for maximizing the rate, it suffices to restrict attention to vectors $v$ of the form $v_{i}=v_{*}\sigma_{i}$ for $1\leq i\leq d$. From the constraint $R^{-1}v\leq\mathbf{b}$, we obtain

[TABLE]

and hence,

[TABLE]

From this observation, it follows that for any such vector $v$, $\Lambda(v)=v_{*}^{2}\leq(a(\Theta))^{-2}$ and $T(v)\geq 1+(a(\Theta))^{2}\log(2d)$. The vector $\widetilde{v}$ given by $\widetilde{v}_{i}:=(a(\Theta))^{-1}\sigma_{i}$ for each $i$ satisfies $R^{-1}\widetilde{v}\leq\mathbf{b}$ and hence, simultaneously maximizes $\frac{D_{1}}{n(R)T(v)}$ and $\frac{\Lambda(v)}{A}$. Also note that $T(\widetilde{v})=1+(a(\Theta))^{2}\log(2d)$ and $\Lambda(\widetilde{v})=(a(\Theta))^{-2}$.

From Lemma 8 with $\widetilde{v}$ in place of $v$, we obtain for any $\mathbf{x}\in\mathbb{R}^{d}_{+}$ and $t\geq t_{0}\left(1+(a(\Theta))^{2}\log(2d)\right)$,

[TABLE]

Consider the dominating process $X^{+}(\cdot;\mathbf{x})=X_{v}^{+}(\cdot;\mathbf{x})$ with $v=\widetilde{v}$. Since $R^{-1}X(t;\mathbf{x})\leq R^{-1}X^{+}(t;\mathbf{x})$, for each $i$,

[TABLE]

and hence,

[TABLE]

Moreover,

[TABLE]

Denote by $X(\infty)$ and $X^{+}(\infty)$ the random vectors sampled from the stationary distribution of $X(\cdot;\mathbf{x})$ and $X^{+}(\cdot;\mathbf{x})$ respectively. By [7], the laws of $X(t;\mathbf{x})$ and $X^{+}(t;\mathbf{x})$ converge in total variation to those of $X(\infty)$ and $X^{+}(\infty)$ respectively. Consequently, $\|X(\infty)\|_{1}$ is stochastically dominated by $\sum_{i,j}(R^{-1})_{ij}X^{+}_{j}(\infty)$ and $\|X(\infty)\|_{\infty,\widetilde{v}}$ is stochastically dominated by $b(\Theta)\|X^{+}(\infty)\|_{\infty,\widetilde{v}}$. As $X_{i}^{+}(\infty)$ is the stationary distribution of a one-dimensional reflected Brownian motion with drift $-\widetilde{v}_{i}$ and variance $\sigma_{i}^{2}$, it is a standard fact that $X_{i}^{+}(\infty)$ follows an exponential distribution with mean $\sigma_{i}^{2}/(2\widetilde{v}_{i})=\sigma_{i}a(\Theta)/2$. This implies

[TABLE]

and

[TABLE]

Using the above estimates we have for any $D_{2}\geq 9$,

[TABLE]

Now fix $D_{2}=\max\{A_{0},9\}$. Using the above estimates and (44) with $A=D_{2}b(\Theta)$ (noting that $b(\Theta)\geq 1$), we obtain for any $\mathbf{x}\in\mathbb{R}^{d}_{+}$ and $t\geq t_{0}\left(1+(a(\Theta))^{2}\log(2d)\right)$,

[TABLE]

which proves the Wasserstein bound in the theorem upon noting that $\|\mathbf{x}\|_{\infty,\widetilde{v}}=(a(\Theta))^{-1}\|\mathbf{x}\|_{\infty}^{*}$. This, in turn, implies the stated bound on the relaxation time. ∎

Remark 8**.**

To obtain the better bound displayed in Remark 3, note that we can replace the bound in Lemma 8 by

[TABLE]

which follows from the second-to-last inequality in the calculation (42). The vector $\widetilde{v}$ in the proof of Theorem 1 still optimizes the above bound over all $v$ and leads to the bound displayed in Remark 3.

Acknowledgements. SB was partially supported by a Junior Faculty Development Award made by UNC, Chapel Hill. Research of AB was partially supported by the National Science Foundation (DMS- 1814894, DMS-1853968). We thank Kavita Ramanan and Andrey Sarantsev for helpful discussions. We also thank two anonymous referees whose comments and suggestions helped to improve the article.