2. MHI samplers and Wasserstein distances

As they will be considered here, MHI samplers simulate a Markov chain with invariant distribution $\Pi$ supported on a nonempty set $\Theta \subseteq \mathbb{R}^d$ using a proposal distribution Q which, to avoid trivialities, is assumed throughout to be different than $\Pi$ . We also assume throughout that $\Pi$ has Lebesgue density $\pi$ with support $\Theta$ , and Q has Lebesgue density q with support $\Theta$ . Define

\begin{equation*} a(\theta, \theta^{\prime}) = \begin{cases} \min\bigg\{ \dfrac{ \pi(\theta^{\prime}) q(\theta) }{\pi(\theta) q(\theta^{\prime})}, 1 \bigg\} & \text{ if }\pi(\theta) q(\theta^{\prime}) > 0, \\ 1 & \text{ if } \pi(\theta) q(\theta^{\prime}) = 0. \end{cases} \end{equation*}

We consider MHI samplers initialized at a point $\theta_0 \in \Theta$ . MHI proceeds as follows: for $t \in \{1, 2, \ldots \}$ , given $\theta_{t-1}$ , draw $\theta^{\prime}_t \sim Q({\cdot})$ and $U_t \sim \text{Unif}(0, 1)$ independently so that

\begin{equation*} \theta_{t} = \begin{cases} \theta^{\prime}_t & \text{if } U_t \le a( \theta_{t-1}, \theta^{\prime}_t ), \\ \theta_{t-1} & \text{otherwise}. \end{cases} \end{equation*}

If $\delta_\theta$ denotes the Dirac measure at the point $\theta$ , the MHI Markov kernel P is defined for $\theta \in \mathbb{R}^d$ and $B \subseteq \mathbb{R}^d$ by

\begin{equation*} P(\theta, B) = \int_{B} a(\theta, \theta^{\prime}) q(\theta^{\prime}) \,\mathrm{d}\theta^{\prime} + \delta_\theta(B) \bigg( 1 - \int a(\theta, \theta^{\prime}) q(\theta^{\prime}) \,\mathrm{d}\theta^{\prime} \bigg) .\end{equation*}

For $\theta \in \mathbb{R}^d$ , define the Markov kernel at iteration time $t \ge 2$ recursively by

\begin{equation*} P^t(\theta, B) = \int P(\theta, \mathrm{d}\theta^{\prime}) P^{t-1}(\theta^{\prime}, B).\end{equation*}

Let $\mathcal{C}(P^t(\theta, \cdot), \Pi)$ be the set of all joint probability measures with marginals $P^t(\theta, \cdot)$ , and $\Pi$ and $\rho$ be a lower semicontinuous metric. The $L_1$ -Wasserstein distance [Reference Kantorovich and Rubinstein24, Reference Villani50, Reference Villani51], which we will call simply the Wasserstein distance, is

\begin{equation*} \mathcal{W}_{\rho}( P^t(\theta, \cdot), \Pi ) = \inf_{\xi \in \mathcal{C}(P^t(\theta, \cdot), \Pi)} \int \rho(\theta^{\prime}, \omega) \,\mathrm{d}\xi(\theta^{\prime}, \omega) .\end{equation*}

Notice that when the metric $\rho$ is $\rho(\theta, \omega) = I_{\theta \not= \omega}$ , then the Wasserstein distance is the total variation distance. More generally, for a lower semicontinuous function $V \ge 1$ , $\rho(\theta, \omega) = [V(\theta) + V(\omega)] I_{\theta \not= \omega}$ defines a weighted total variation distance. Another example is $\rho(\theta, \omega) = \min\{ \left\lVert \theta - \omega \right\rVert, 1 \}$ , which is always less than the Hamming metric used in total variation. More general $L_p$ -Wasserstein distances with $p \ge 2$ are not studied in this work.

Wasserstein distances control the bias of integrals of $P^t(\theta, \cdot)$ over all $\rho$ -Lipschitz functions, whereas total variation controls the bias for bounded, measurable functions. Often in applications, various types of Lipschitz functions are of interest. Consider, for example, functions such as the identity function and Winsorized functions such as $g(\omega) = ({-}R) \vee (R \wedge \omega)$ where $R \in (0, +\infty)$ , which reduce sensitivity to extreme values. An alternative motivation for using Wasserstein distances is that convergence analysis in certain Wasserstein distances may improve the scaling to high dimensions [Reference Hairer, Stuart and Vollmer15, Reference Qin and Hobert39, Reference Qin and Hobert40]. Serious problems can arise with existing convergence analysis in total variation [Reference Tierney48] in high dimensions even for seemingly trivial examples such as the one below.

At a specific initialization, we may look for a smaller convergence rate which may scale to high dimensions in a Wasserstein distance weaker than total variation. We will not only develop an exact convergence analysis which controls the bias for Lipschitz functions, but, roughly speaking, we also discover that the convergence rate in Example 1 is exact for Wasserstein distances and the same for every initialization. Nevertheless, we show that convergence rates can scale to large problem sizes using a novel proposal and exact convergence analysis.

3. Exact convergence rates for MHI samplers

When the ratio of the proposal and target densities is bounded below by a positive number, i.e. $\varepsilon^* = \inf_{\theta \in \Theta} \{ q(\theta)/\pi(\theta) \} > 0$ , the MHI sampler is uniformly ergodic in total variation with convergence rate upper bounded by $1 - \varepsilon^*$ [Reference Tierney48, Corollary 4]. Unlike in accept–reject sampling, $\varepsilon^*$ does not need to be known explicitly or computed in order to implement MHI. However, this requirement was shown to be necessary for uniform ergodicity in total variation [Reference Mengersen and Tweedie29, Theorem 2.1]. More recently, it was shown that the convergence rate cannot be improved [Reference Wang52, Theorem 1]. We show this to be the case even in weaker Wasserstein distances where the lower bound does not follow trivially from that of the total variation lower bound [Reference Wang52, Theorem 1].

Theorem 1. Suppose $\rho(\cdot, \cdot) \le 1$ . Then

\begin{align*} \sup_{\theta \in \Theta} \mathcal{W}_\rho(P^t(\theta, \cdot), \Pi) \le (1 - \varepsilon^*)^t \sup_{\theta \in \Theta} \int \rho(\cdot, \theta) \,\mathrm{d}\Pi. \end{align*}

If, in addition, q is lower semicontinuous on $\Theta$ , $\pi$ is upper semicontinuous on $\Theta$ , and $\Theta$ can be expressed as a countable union of compact sets, then

\begin{align*} (1 - \varepsilon^*)^t \inf_{\theta \in \Theta} \int \rho(\cdot, \theta) \,\mathrm{d}\Pi \le \sup_{\theta \in \Theta} \mathcal{W}_\rho(P^t(\theta, \cdot), \Pi) \le (1 - \varepsilon^*)^t \sup_{\theta \in \Theta} \int \rho(\cdot, \theta) \,\mathrm{d}\Pi. \end{align*}

The semicontinuity assumption is not required when working with the total variation distance [Reference Wang52, Theorem 1], but it is a mild assumption that holds in many practical applications. The upper bound constant can improve upon upper bounds in total variation [Reference Tierney48] if, for example, $\rho$ is continuous and $\Theta$ is compact. If $\varepsilon^* = 0$ , Theorem 1 also gives the lower bound

\begin{equation*} \inf_{\theta \in \Theta} \int \rho(\cdot, \theta) \,\mathrm{d}\Pi \le \sup_{\theta \in \Theta} \mathcal{W}_\rho(P^t(\theta, \cdot), \Pi),\end{equation*}

which shows that MHI cannot converge uniformly from any starting point for many Wasserstein distances. Thus, under mild assumptions, Theorem 1 gives a complete characterization of the worst-case convergence of the MHI sampler in many Wasserstein distances.

Exact convergence expressions are available when the Markov chain is initialized at $\theta^* = \text{argmin}\{ q(\theta) / \pi(\theta) \,:\, \theta \in \Theta\}$ using techniques from [Reference Wang52].

Proposition 1. Suppose there exists a solution $\theta^* = \mathrm{argmin}\{ q(\theta) / \pi(\theta) \,:\, \theta \in \Theta\}$ . Then

\begin{equation*} \mathcal{W}_{\rho}( P^t(\theta^*, \cdot), \Pi ) = ( 1 - q(\theta^*)/\pi(\theta^*))^t \int \rho(\theta, \theta^*) \,\mathrm{d}\Pi(\theta) . \end{equation*}

Proof. Define $\varepsilon_{\theta^*} = q(\theta^*)/\pi(\theta^*) = \varepsilon^*$ . Under our assumptions, $P^t(\theta^*, \cdot)$ can be represented as a convex combination of the target distribution and the Dirac measure at the point $\theta^*$ [Reference Wang52, Remark 1, Theorem 2], that is,

\begin{equation*} P^t(\theta^*, \cdot) = (1 - (1 - \varepsilon_{\theta^*})^t) \Pi + (1 - \varepsilon_{\theta^*})^t \delta_{\theta^*} . \end{equation*}

Let $\psi \,:\, \Theta \to \mathbb{R}$ be a function such that $\int_\Theta |\psi| \,\mathrm{d}\Pi < \infty$ . We have the identity

\begin{equation*} \int_\Theta \psi \,\mathrm{d} P^t(\theta^*, \cdot) = (1 - (1 - \varepsilon_{\theta^*})^t) \int_\Theta \psi \,\mathrm{d}\Pi + (1 - \varepsilon_{\theta^*})^t \psi(\theta^*). \end{equation*}

Since the only coupling between $\Pi^*$ and the Dirac measure $\delta_{\theta^*}$ is the independent coupling [Reference Giraudo14], the Wasserstein distance takes the simple form $\mathcal{W}_{\rho}( \delta_{\theta^*}, \Pi) = \int \rho(\theta, \theta^*) \,\mathrm{d}\Pi(\theta)$ .

Since q is not exactly $\pi$ , $\varepsilon_{\theta^*} \in (0, 1)$ . Let $M_b(\mathbb{R}^d)$ be the set of bounded measurable functions on $\mathbb{R}^d$ and, for real-valued functions $\varphi$ , let $\left\lVert\varphi\right\rVert_{\text{Lip}(\rho)} = \sup_{x, y, x \not= y} \{ |\varphi(x) - \varphi(y)| / \rho(x, y) \}$ denote the Lipschitz norm with respect to the distance $\rho$ . Applying the Kantorovich–Rubinstein theorem [Reference Villani50, Theorem 1.14],

\begin{align*} \mathcal{W}_{\rho}(P^t(\theta^*, \cdot), \Pi) & = \sup_{\substack{\varphi \in M_b(\mathbb{R}^d) \\ \left\lVert\varphi\right\rVert_{\text{Lip}(\rho)} \le 1}} \int_{\Theta} \varphi \,\mathrm{d}(P^t(\theta^*, \cdot) - \Pi) \\ & = \sup_{\substack{\varphi \in M_b(\mathbb{R}^d) \\ \left\lVert\varphi\right\rVert_{\text{Lip}(\rho)} \le 1}} \bigg\{(1 - \varepsilon_{\theta^*})^t \int_{\Theta} \varphi \,\mathrm{d}(\delta_{\theta^*} - \Pi) \bigg\} \\ & = (1 - \varepsilon_{\theta^*})^t \sup_{\substack{\varphi \in M_b(\mathbb{R}^d) \\ \left\lVert\varphi\right\rVert_{\text{Lip}(\rho)} \le 1}} \int_{\Theta} \varphi \,\mathrm{d} \left( \delta_{\theta^*} - \Pi \right) \\ & = (1 - \varepsilon_{\theta^*})^t \mathcal{W}_{\rho}(\delta_{\theta^*}, \Pi) = (1 - \varepsilon_{\theta^*})^t \int_{\Theta} \rho(\theta, \theta^*) \,\mathrm{d}\Pi(\theta). \end{align*}

3.1. Application: Bayesian quantile regression

Fix $r \in (0, 1)$ and suppose, for $i=1,\ldots,n$ , that $\varepsilon_i$ are independent and identically distributed (i.i.d.) with density $p_r(\varepsilon) = r(1 - r) (\exp\!((1 - r) \varepsilon ) I_{\varepsilon < 0} + \exp\!({-}r\varepsilon) I_{\varepsilon \ge 0})$ . Let $v_0, s_0 \in (0, \infty)$ and $C \in \mathbb{R}^{d \times d}$ be symmetric positive-definite. We parameterize the inverse gamma distribution so that if $\sigma \sim \text{IG}(v, s)$ for some $v, c \in (0, \infty)$ , then it has a density proportional to $\sigma^{-v - 1} \exp\!({-}s/\sigma)$ . Assume the Bayesian quantile regression model for $i \in 1, \ldots, n$ where $X_i \in \mathbb{R}^d$ is fixed and

\begin{equation*} \sigma \sim \text{IG}(v_0, s_0), \qquad \beta \mid \sigma \sim \text{N}_d(0, \sigma C) , \qquad Y_i = \beta^\top X_i + \sigma \varepsilon_i.\end{equation*}

Let $\Pi(\cdot \mid X, Y)$ denote the posterior and $\pi(\cdot \mid X, Y)$ denote the density for this Bayesian model with normalizing constant $Z_{\Pi(\cdot \mid X, Y)}$ .

Upper bounds on the convergence rate were previously investigated for Gibbs samplers [Reference Khare and Hobert25] in this setting. We will study the MHI sampler with a normal–inverse-gamma proposal constructed as follows. Define the convex function by $\ell_r(u) = u (r - I_{u < 0})$ and $s_{n, r}\,:\, \mathbb{R}^d \to \mathbb{R}$ by $s_{n, r}(\beta) = \sum_{i = 1}^n \ell_r(Y_i - \beta^\top X_i) + \beta^\top C^{-1} \beta / 2$ . Since $s_{n, r}$ is strongly convex, let $\beta^* \in \mathbb{R}^d$ be the unique minimum of the function $s_{n, r}$ . Now the MHI proposal is given by

\begin{equation*} \sigma \sim \text{IG}(n + v_0, s_0 + s_{n, r}(\beta^*)) , \qquad \beta \mid \sigma \sim \text{N}_d(\beta^*, \sigma C).\end{equation*}

Let $\Gamma \,:\, (0, \infty) \to \mathbb{R}$ be the usual Gamma function and define

\begin{equation*} \varepsilon_{\beta^*} = Z_{\Pi(\cdot \mid X, Y)}(2 \pi)^{-{d}/{2}} \det\!(C)^{-1/2}(s_0 + s_{n, r}(\beta^*))^{n + v_0} \Gamma(n + v_0)^{-1}.\end{equation*}

The following gives an exact convergence rate of this algorithm which completely characterizes its convergence from a specific initialization.

Theorem 2. For any $\sigma_0 \in (0, \infty)$ ,

\begin{equation*} \mathcal{W}_{\rho}(P^t((\beta^*, \sigma_0), \cdot), \Pi(\cdot \mid X, Y)) = (1 - \varepsilon_{\beta^*})^t \int \rho((\beta, \sigma), (\beta^*, \sigma_0)) \,\mathrm{d}\Pi(\beta, \sigma \mid X, Y). \end{equation*}

Proof. We may define the function $f \,:\, \mathbb{R}^d \times (0, \infty) \to \mathbb{R}$ by

\begin{equation*} f(\beta, \sigma) = \frac{s_0 + s_{n, r}(\beta)}{\sigma} + (n + v_0 + 1 + d/2) \log\!(\sigma) \end{equation*}

and write the posterior density as $\pi(\beta, \sigma \mid X, Y) = Z_{\Pi(\cdot \mid X, Y)}^{-1} \exp\!({-}f(\beta, \sigma))$ . Since the function $\beta \mapsto s_{n, r}(\beta) - \beta^\top C^{-1} \beta / 2$ is a convex function on $\mathbb{R}^d$ , by Lemma 4, for every $\beta \in \mathbb{R}^d$ , $s_{n, r}(\beta) \ge s_{n, r}(\beta^*) + \frac{1}{2}(\beta - \beta^*)^\top C^{-1}(\beta - \beta^*)$ . For any $(\beta, \sigma) \in \mathbb{R}^d \times (0, \infty)$ , we then have the lower bound

\begin{align*} f(\beta, \sigma) & = \frac{s_0 + s_{n, r}(\beta)}{\sigma} + (n + v_0 + 1 + d/2) \log\!(\sigma) \\ & \ge \frac{s_0 + s_{n, r}(\beta^*)}{\sigma} + (n + v_0 + 1 + d/2) \log\!(\sigma) + \frac{1}{2 \sigma}(\beta - \beta^*)^\top C^{-1}(\beta - \beta^*). \end{align*}

This implies that

\begin{equation*} f(\beta, \sigma) - \frac{1}{2 \sigma}(\beta - \beta^*)^\top C^{-1}(\beta - \beta^*) - \frac{s_0 + s_{n, r}(\beta^*)}{\sigma} - (n + v_0 + 1 + d/2) \log\!(\sigma) \ge 0. \end{equation*}

Let q denote the proposal’s normal–inverse-gamma density. For any $\sigma_0 \in (0, \infty)$ and for every $(\beta, \sigma) \in \mathbb{R}^d \times (0, \infty)$ , we have shown that

\begin{equation*} \frac{q(\beta, \sigma)}{\pi(\beta, \sigma)} \ge Z_{\Pi(\cdot \mid X,Y)}(2\pi)^{-{d}/{2}}\det\!(C)^{-1/2}(s_0 + s_{n,r}(\beta^*))^{n + v_0}\Gamma(n + v_0)^{-1} = \frac{q(\beta^*, \sigma_0)}{\pi(\beta^*, \sigma_0)}. \end{equation*}

An application of Proposition 1 completes the proof.

Note that $\varepsilon_{\beta^*}$ is difficult to compute since it depends on the normalizing constant, but we give an example later where upper bounding the convergence rate is possible in Bayesian logistic and probit regression.

4. The convergence rate at arbitrary initializations

The previous section studies the worst-case convergence rate and the convergence rate at an individual point for the MHI sampler. We can study the convergence rate at every point, as was done for total variation [Reference Wang52]. The technique needed to prove this relies on the exact representation of the MHI sampler [Reference Smith and Tierney46, Theorem 1], but new techniques are needed to show this in the Wasserstein distance. Similar to the convergence rate in total variation [Reference Wang52], we define the Wasserstein convergence rate for a point $\theta \in \Theta$ as $r_{\rho}(\theta) = \lim_{t \to \infty} \mathcal{W}_{\rho}(P^t(\theta, \cdot), \Pi)^{1/t}$ .

When the distance metric is $\min\{ \left\lVert\cdot\right\rVert, 1 \}$ where $\left\lVert\cdot\right\rVert$ can be any norm, Theorem 3 shows we can obtain the convergence rate at every point under mild conditions. We require $\pi$ and q to be locally $\left\lVert\cdot\right\rVert$ -Lipschitz continuous and bounded, which is stronger than only locally Lipschitz as in total variation [Reference Wang52]. However, this additional condition is satisfied in many practical applications in statistics. Theorem 3 also lower bounds the convergence rate for Wasserstein $L_p$ -distances, and the rate of convergence for these distances cannot be improved.

Proof. If the initialization is at $\theta^*$ , then the result follows from Proposition 1. Fix a point $\theta_0 \in \Theta$ such that $\theta_0 \not= \theta^*$ . Using Lemma 1, we have an upper bound on the convergence rate by

\begin{equation*} \limsup_{t \to \infty}\mathcal{W}_{\min\{\left\lVert\cdot\right\rVert, 1\}}(P^t(\theta_0,\cdot),\Pi)^{1/t} \le 1 - q(\theta^*)/\pi(\theta^*). \end{equation*}

It remains to lower bound the limit inferior.

Denote the standard p-norms for vectors $x \in \mathbb{R}^d$ by $\left\lVert x\right\rVert_p = \big(\sum_{i} |x_i|^p\big)^{1/p}$ . For $h \in (0, 1]$ , define the function $\varphi_h(\theta) = (2h)^{-d}\exp\!({-}h^{-1} \left\lVert\theta - \theta^*\right\rVert_1)$ . This is the probability density function for a Laplace distribution and $\varphi_h$ is $2^{-d} h^{-d - 1}$ $\left\lVert\cdot\right\rVert_1$ -Lipschitz, and so $2^{d} h^{d + 1} \varphi_h(\theta) = h \exp\!({-}h^{-1} \left\lVert\theta - \theta^*\right\rVert_1)$ is nonnegative, $\left\lVert\cdot\right\rVert_1$ -Lipschitz with constant 1 and bounded by 1. In particular, it is readily shown that $2^{d} h^{d + 1} \varphi_h(\theta)$ is $\min\{1, \left\lVert\cdot\right\rVert_1\}$ -Lipschitz. Using Kantorovich–Rubinstein duality [Reference Villani50, Theorem 1.14], we have the lower bound

(1) \begin{equation} \mathcal{W}_{\min\{\left\lVert\cdot\right\rVert_1, 1\}}(P^t(\theta_0, \cdot), \Pi) \ge 2^{d} h^{d + 1} \bigg[\int \varphi_h \,\mathrm{d}\Pi - \int \varphi_h \,\mathrm{d} P^t(\theta_0, \cdot) \bigg] . \end{equation}

We will develop some approximation properties of $\varphi_h$ . Since $\theta_0 \not= \theta^*$ , it is readily shown that $\lim_{h \downarrow 0} \varphi_h(\theta_0) = 0$ . Using a change of variables and since we have assumed $\sup_{\theta} \pi(\theta) < \infty$ ,

(2) \begin{align} \bigg| \int \varphi_h(\theta^{\prime}) \pi(\theta^{\prime}) \,\mathrm{d}\theta^{\prime} - \pi(\theta^*) \bigg| & \le \int_{\left\lVert\theta^{\prime}\right\rVert_2 \le t} |\pi(\theta^* + h \theta^{\prime}) - \pi(\theta^*)|2^{-d} \exp\!({-}\left\lVert\theta^{\prime}\right\rVert_1) \,\mathrm{d}\theta^{\prime}\end{align}

(3) \begin{align}& \qquad\qquad\qquad\qquad + 2 \sup_{\theta} \pi(\theta) \int_{\left\lVert\theta^{\prime}\right\rVert_2 > t} 2^{-d} \exp\!({-}\left\lVert\theta^{\prime}\right\rVert_1) \,\mathrm{d}\theta^{\prime}. \end{align}

If $Y_1, \ldots, Y_d$ are i.i.d. Laplace, then we have the tail bound

\begin{equation*} \mathbb{P}(\left\lVert Y\right\rVert_2 \ge t) \le \mathbb{P}\bigg( \max_{i} |Y_i| \ge \frac{t}{\sqrt{d}\,} \bigg) \le \sum_{i = 1}^d \mathbb{P}\bigg( |Y_i| \ge \frac{t}{\sqrt{d}\,} \bigg) \le d \exp\bigg({-}\frac{t}{\sqrt{d}}\bigg). \end{equation*}

Since norms are equivalent on $\mathbb{R}^d$ , $\pi$ and q are locally Lipschitz with respect to any norm. Choosing $h = h_0/t^2$ for some $h_0 \in (0, 1)$ , since $\pi$ is locally $\left\lVert\cdot\right\rVert_2$ -Lipschitz, we can find a universal constant $L \in (0, \infty)$ such that, for any $\left\lVert\theta^{\prime}\right\rVert_2 \le t$ , $|\pi(\theta^* + h \theta^{\prime}) - \pi(\theta^*)| \le {h_0 L}/{t}$ . Applying these upper bounds to (2) and (3), for large enough t we have

(4) \begin{equation} \bigg|\int \varphi_h(\theta^{\prime}) \pi(\theta^{\prime}) \,\mathrm{d}\theta^{\prime} - \pi(\theta^*)\bigg| \le \frac{2 h_0 L}{t}. \end{equation}

A similar argument with the assumptions on q yields, for large enough t,

(5) \begin{equation} \bigg|\int \varphi_h(\theta^{\prime}) q(\theta^{\prime}) \,\mathrm{d}\theta^{\prime} - q(\theta^*)\bigg| \le \frac{2 h_0 L}{t}. \end{equation}

It remains to lower bound (1). We will use the exact representation of the independence sampler [Reference Smith and Tierney46, Theorem 1, Lemma 3]. Define the importance sampling weight by $w(\theta) = \pi(\theta) / q(\theta)$ and its maximum $w^* = w(\theta^*) = \pi(\theta^*) / q(\theta^*)$ . For $w \in (0, \infty)$ , define

\begin{equation*} \lambda(w) = \int_{w(\theta^{\prime}) \le w} \bigg[1 - \frac{w(\theta^{\prime})}{w}\bigg] q(\theta^{\prime}) \,\mathrm{d}\theta^{\prime}, \qquad T_t(w) = \int_{w}^{\infty} \frac{t \lambda^{t-1}(v)}{v^2} \,\mathrm{d} v, \end{equation*}

and, using the exact representation of the independence sampler [Reference Smith and Tierney46, Theorem 1, Lemma 3], for measurable sets $B \subseteq \mathbb{R}^d$ ,

\begin{equation*} P^t(\theta_0, B) = \int_B T_t(\!\max\{w(\theta_0), w(\theta^{\prime})\}) \pi(\theta^{\prime}) \,\mathrm{d}\theta^{\prime} + \lambda^t(w(\theta_0)) \delta_{\theta_0}(B). \end{equation*}

The proof of existing results [Reference Wang52, Theorem 4] shows that

\begin{equation*} T_t(w) \le 1 + (1 - 1/w^*)^t \bigg[\frac{t}{w^* - 1}\bigg(\frac{w^*}{w} - 1\bigg) - 1\bigg]. \end{equation*}

We now use the exact representation of the independence sampler to lower bound (1). We have $\lambda^t(w(\theta_0)) \le (1 - 1/w^*)^t$ , so we then have the upper bound

\begin{align*} \int\! \varphi_h \,\mathrm{d} P^t(\theta_0, \cdot) - \int\! \varphi_h \,\mathrm{d}\Pi & \le -(1 - 1/w^*)^t \!\int\! \varphi_h \,\mathrm{d}\Pi + \lambda^t(w(\theta_0)) \varphi_h(\theta_0) \\ & \quad + (1 - 1/w^*)^t \frac{t}{w^* - 1} \int \bigg[\frac{w^*}{\max\{w(\theta_0),w(\theta^{\prime})\}} - 1\bigg]\varphi_h(\theta^{\prime})\pi(\theta^{\prime})\,\mathrm{d}\theta^{\prime} \\ & \le (1 - 1/w^*)^t \bigg[\varphi_h(\theta_0) - \int \varphi_h \,\mathrm{d}\Pi\bigg] \\ & \quad + (1 - 1/w^*)^t \frac{t}{w^* - 1} \bigg[ w^* \!\int\! \varphi_h(\theta^{\prime}) q(\theta^{\prime}) \,\mathrm{d}\theta^{\prime} - \int\! \varphi_h(\theta^{\prime}) \pi(\theta^{\prime}) \,\mathrm{d}\theta^{\prime} \bigg]. \end{align*}

Using this upper bound with the approximations (4) and (5) yields

\begin{align*} \int \varphi_h \,\mathrm{d} P^t(\theta_0, \cdot) - \int \varphi_h \,\mathrm{d}\Pi & \le (1 - 1/w^*)^t \bigg[ \varphi_h(\theta_0) - \pi(\theta^*) + \frac{2 h_0 L}{t} \bigg] \\ & \quad + (1 - 1/w^*)^t \frac{t}{w^* - 1} \bigg[ w^* q(\theta^*) - \pi(\theta^*) + (w^* + 1) \frac{2 h_0 L}{t} \bigg] \\ & \le (1 - 1/w^*)^t \bigg[ \varphi_h(\theta_0) - \pi(\theta^*) + \frac{2 h_0 L}{t} + \frac{w^* + 1}{w^* - 1} 2 h_0 L \bigg]. \end{align*}

Therefore, we can choose a small enough $h_0 \in (0, 1)$ independently of t such that we have the upper bound

\begin{equation*} \int \varphi_h \,\mathrm{d} P^t(\theta_0, \cdot) - \int \varphi_h \,\mathrm{d}\Pi \le -\frac{\pi(\theta^*)}{4} (1 - 1/w^*)^t. \end{equation*}

Applying these bounds to (1) and using that we have chosen $h = h_0/t^2$ , we have the lower bound

\begin{align*} \mathcal{W}_{ \min\{ \left\lVert\cdot\right\rVert_1, 1 \}}(P^t(\theta_0, \cdot), \Pi) & \ge 2^{d} h^{d + 1} \bigg[ \int \varphi_h \,\mathrm{d}\Pi - \int \varphi_h \,\mathrm{d} P^t(\theta_0, \cdot) \bigg] \\ & \ge 2^{d - 2} h^{d + 1} \pi(\theta^*) (1 - 1/w^*)^t \\ & \ge 2^{d - 2} \bigg( \frac{h_0}{t^2} \bigg)^{d + 1} \pi(\theta^*) (1 - 1/w^*)^t. \end{align*}

Taking the limit,

\begin{equation*} \liminf_{t \to \infty} \mathcal{W}_{ \min\{ \left\lVert\cdot\right\rVert_1, 1 \}}(P^t(\theta_0, \cdot), \Pi)^{1/t} \ge 1 - 1/w^* = 1 - q(\theta^*) / \pi(\theta^*). \end{equation*}

Since all norms are equivalent on $\mathbb{R}^d$ , this can be extended to any norm $\left\lVert\cdot\right\rVert$ .

5. MHI samplers with centered Gaussian proposals

We look to apply the exact convergence expression from Proposition 1 to practical applications since the convergence rate is the same at every initialization under mild assumptions. Recently, centered drift functions have been used to improve convergence analyses of some Monte Carlo Markov chains [Reference Ekvall and Jones11, Reference Qin and Hobert37, Reference Qin and Hobert40]. Our focus is instead on centering the proposal distribution, that is, matching the optimal points of the proposal and target densities similar to Laplace approximations.

We shall see in the next section that by centering a Gaussian proposal, we may satisfy the assumptions of Proposition 1 for a general class of target distributions with $\theta^*$ being the optimum of the target’s density. While we focus on Gaussian proposals, the technique of centering proposals is in fact more general.

We will assume the target distribution $\Pi$ is a probability distribution supported on $\mathbb{R}^d$ . With $f \,:\, \mathbb{R}^d \to \mathbb{R}$ and normalizing constant $Z_{\Pi}$ , define the density $\pi$ by $\pi(\theta) = Z_{\Pi}^{-1} \exp\!({-}f(\theta))$ . Let $\theta^*$ be the unique maximum of $\pi$ , $\alpha \in (0, +\infty)$ , and $C \in \mathbb{R}^{d \times d}$ be a symmetric, positive-definite matrix. Let the proposal distribution Q with density q correspond to a d-dimensional Gaussian distribution, $\text{N}_{d}(\theta^*, \alpha^{-1} C )$ . In this case, the ratio of the proposal density and target density is

\begin{equation*}\varepsilon_{\theta^*} = (2 \pi)^{-d/2} \alpha^{d/2} \det\!(C)^{-1/2} Z_{\Pi} \exp\!(f(\theta^*)).\end{equation*}

Proposition 2. If $\theta^*$ exists and, for any $\theta \in \mathbb{R}^d$ , $f(\theta) \ge f(\theta^*) + \alpha(\theta - \theta^*)^\top C^{-1}(\theta - \theta^*)/2$ , then

\begin{equation*} \mathcal{W}_{\rho}(P^t(\theta^*, \cdot), \Pi) = (1 - \varepsilon_{\theta^*})^t \int \rho(\theta, \theta^*) \,\mathrm{d}\Pi(\theta). \end{equation*}

Proof. Since the proposal density has been centered at the point $\theta^*$ , it then satisfies $q(\theta^*) = (2 \pi)^{-d/2} \alpha^{d/2} \det\!(C)^{-1/2}$ . For every $\theta \in \mathbb{R}^d$ , we have the lower bound

\begin{align*} \frac{q(\theta)}{\pi(\theta)} & = (2 \pi)^{-d/2} \alpha^{d/2} \det\!(C)^{-1/2} Z_{\Pi} \exp\bigg( f(\theta) - \frac{\alpha}{2} (\theta - \theta^*)^\top C^{-1}(\theta - \theta^*)\bigg) \\ & \ge (2 \pi)^{-d/2} \alpha^{d/2} \det\!(C)^{-1/2} Z_{\Pi} \exp\!(f(\theta^*)) = \frac{q(\theta^*)}{\pi(\theta^*)}. \end{align*}

Since both densities are positive and the proposal is independent of the previous iteration, we have shown that the conditions for Proposition 1 are satisfied and an application of Proposition 1 with the proposal and target distribution Q and $\Pi$ as we have defined them in this section completes the proof.

The assumption in Proposition 2 is sometimes referred to as a quadratic growth condition at $\theta^*$ . This condition does not require convexity of f, but is satisfied and the point $\theta^*$ is guaranteed to exist if the function f satisfies a strong convexity property. A function $h \,:\, \mathbb{R}^d \to \mathbb{R}$ is strongly convex with parameter $\mu$ if there is a $\mu \in (0, +\infty)$ such that $h({\cdot}) - \mu \left\lVert\cdot\right\rVert^2 / 2$ is convex [Reference Hiriart-Urruty and Lemaéchal17, Reference Nesterov32]. The norm in this definition is often taken to be the Euclidean norm, but we will use the norm induced by the matrix $C^{-1}$ . We consider using a Gaussian proposal centered at a point $\theta_0$ which is not necessarily the optimum of the target density. Let $g_{f(\theta_0)} \in \mathbb{R}^d$ be a subgradient of f at $\theta_0$ [Reference Nesterov32]. For a point $\theta_0 \in \mathbb{R}^d$ , we consider the proposal corresponding to a d-dimensional Gaussian distribution, $\text{N}_{d}(\theta_0 - \alpha^{-1} C g_{f(\theta_0)}, \alpha^{-1} C)$ . When f is differentiable, this construction of the proposal uses the gradient of f in a similar way to MALA. The ratio of the proposal and target density evaluated at $\theta_0$ is $\varepsilon_{\theta_0} = (2 \pi)^{-d/2} \det\!(\alpha^{-1} C)^{-1/2} Z_{\Pi} \exp\!(f(\theta_0) - {g_{f(\theta_0)}}^\top C g_{f(\theta_0)} / (2\alpha) )$ . Choosing $\theta_0 \equiv \theta^*$ maximizes the convergence rate and yields the centered Gaussian proposal, but we also have an exact convergence expression in other cases.

Proof. Since the function $f(\theta) - \alpha \theta^\top C^{-1} \theta/2$ is convex for all points on $\mathbb{R}^d$ , for each $\theta \in \mathbb{R}^d$ ,

\begin{align*} f(\theta) & \ge f(\theta_0) + {g_{f(\theta_0)}}^\top (\theta - \theta_0) + \frac{\alpha}{2}(\theta - \theta_0)^\top C^{-1} (\theta - \theta_0) \\ & = f(\theta_0) - \frac{1}{2 \alpha} (C g_{f(\theta_0)})^\top g_{f(\theta_0)} + \frac{\alpha}{2}(\theta - \theta_0 + \alpha^{-1} C g_{f(\theta_0)})^\top C^{-1}(\theta - \theta_0 + \alpha^{-1} C g_{f(\theta_0)}). \end{align*}

This implies that, for every $\theta \in \mathbb{R}^d$ , the ratio of the proposal density q corresponding to the distribution $\text{N}_{d}(\theta_0 - \alpha^{-1} C g_{f(\theta_0)}, \alpha^{-1} C)$ and target density $\pi$ satisfies

\begin{equation*} \frac{q(\theta)}{\pi(\theta)} \ge \frac{q(\theta_0)}{\pi(\theta_0)} = \varepsilon_{\theta_0}. \end{equation*}

An application of Proposition 1 completes the proof.

5.1. Application: Bayesian generalized linear models

We consider Bayesian Poisson and negative-binomial regression for count response regression and Bayesian logistic and probit regression for binary response regression. Suppose there are n discrete-valued responses $Y_i$ with features $X_i \in \mathbb{R}^d$ , and a parameter $\beta \in \mathbb{R}^d$ . For Poisson regression, assume the $Y_i$ are conditionally independent with $Y_i \mid X_i, \beta \sim \text{Poisson}(\exp\!(\beta^\top X_i))$ . Similarly, for negative-binomial regression, if $\xi \in (0, +\infty)$ , assume $Y_i \mid X_i, \beta \sim \text{Negative-Binomial}(\xi,(1 + \exp\!({-}\beta^\top X_i))^{-1})$ . For binary response regression, if $S \,:\, \mathbb{R} \to (0, 1)$ , assume $Y_{i} \mid X_{i}, \beta \sim \text{Bernoulli} (S(\beta^\top X_i))$ . For logistic regression, we will consider $S(x) = (1 + \exp\!(x))^{-1}$ , and for probit regression, we will consider S(x) to be the cumulative distribution function of a standard Gaussian random variable.

Suppose $\beta \sim \text{N}_{d}(0, \alpha^{-1} C)$ , where $\alpha \in (0, +\infty)$ and $C \in \mathbb{R}^{d \times d}$ is a symmetric, positive-definite matrix. Both $\alpha$ and C are assumed to be known. Define the vector $Y = (Y_1, \ldots, Y_n)^\top$ and the matrix $X = (X_1, \ldots, X_n)^\top$ . Let $\Pi(\cdot \mid X, Y)$ denote the posterior with density $\pi(\cdot \mid X, Y)$ . If $\ell_n$ denotes the negative log-likelihood for each model, the posterior density is characterized by

\begin{equation*} \pi(\beta \mid X, Y) = Z_{\Pi(\cdot \mid X,Y))}^{-1}\exp\bigg({-}\ell_n(\beta) - \frac{\alpha}{2}\beta^\top C^{-1}\beta\bigg).\end{equation*}

Observe that the function $\ell_n$ is convex in all four models we consider. Let $\beta^*$ denote the unique maximum of $\pi(\cdot \mid X, Y)$ . For the MHI algorithm, we use an $\text{N}_{d}(\beta^*, \alpha^{-1} C)$ proposal distribution, and Proposition 3 immediately yields the following for each posterior.

Corollary 1. $\mathcal{W}_{\rho}(P^t(\beta^*, \cdot), \Pi(\cdot \mid X, Y)) = (1 - \varepsilon_{\beta^*})^t \int \rho(\beta, \beta^*) \,\mathrm{d}\Pi(\beta \mid X, Y)$ , where $\varepsilon_{\beta^*} = \exp\!(\ell_n(\beta^*) + ({\alpha}/{2}){\beta^*}^\top C^{-1} \beta^*) Z_{\Pi(\cdot \mid X, Y)} ((2 \pi)^{d/2} \det\!(\alpha^{-1} C)^{1/2})^{-1}$ .

5.2. Convergence complexity analysis in binary response regression

For the MHI algorithm, we continue to use a centered proposal $\text{N}_{d}(\beta^*, \alpha^{-1} C)$ and first consider a more general posterior density of the form $\pi(\beta \mid X, Y) \propto \exp\!({-}\ell_n(\beta) - \alpha \beta^\top C^{-1} \beta/ 2)$ depending on data X, Y of size n. We also assume the limit of the trace of the covariance matrix used in our prior to be finite, i.e. $\mathrm{tr}(C) \to s_0 \in (0, +\infty)$ as $d \to +\infty$ . Note that the trace of the covariance being finite is a necessary condition for Gaussian distributions to exist in an infinite-dimensional Hilbert space [Reference Bogachev3].

Theorem 4. Suppose that the following conditions hold for a sequence $d_n, n \to \infty$ :

1. (i) The negative log-likelihood $\ell_n$ is a twice continuously differentiable convex function.

2. (ii) There is a universal constant $H_0 \in (0, +\infty)$ such that the largest eigenvalue of the Hessian of the negative log-likelihood $H_{\ell_n}$ satisfies, for every $\beta \in \mathbb{R}^d$ , $\limsup_{d_n, n \to \infty} \lambda_\mathrm{max}(H_{\ell_n}(\beta)) \le H_0$ .

Then $\limsup_{d_n, n \to \infty}\mathcal{W}_{\rho}(P^t(\beta^*, \cdot), \Pi(\cdot \mid X, Y)) \le M_0 (1 - \exp\!({-}H_0 s_0 / (2 \alpha)))^t$ , where $M_0 = \limsup_{d_n, n \to \infty} \int \rho(\beta, \beta^*) \,\mathrm{d}\Pi(\beta \mid X, Y)$ .

Proof. Define the function $f \,:\, \mathbb{R}^d \to \mathbb{R}$ by $f(\beta) = \ell_n(\beta) + \alpha \beta^\top C^{-1} \beta / 2$ , where $\ell_n$ is the negative log-likelihood loss function, and define $Z_Q = (2 \pi)^{d/2} \det\!(\alpha^{-1} C)^{1/2}$ . We first lower bound the quantity $\exp\!(f(\beta^*)) Z_{\Pi( \cdot \mid X, Y)} / Z_{Q}$ . For $\varepsilon \in (0, 1)$ and sufficiently large $d_n, n$ , we have, for any $\beta \in \mathbb{R}^d$ and any $v \in \mathbb{R}^d$ , $v^\top H_{\ell_n}(\beta)v \le (1 + \varepsilon)H_0\left\lVert v\right\rVert_2^2$ . This implies that, for any $\beta \in \mathbb{R}^d$ and any $v \in \mathbb{R}^d$ , the Hessian of f, denoted by $H_{f}$ , satisfies $v^\top H_{f}(\beta) v \le v^\top ((1 + \varepsilon) H_0 I_d + \alpha C^{-1}) v$ . Since the function $\ell_n$ is twice continuously differentiable, f is also twice continuously differentiable. Since both the gradient $\nabla f$ and Hessian $H_{f}$ are continuous and $\nabla f(\beta^*) = 0$ , we use a Taylor expansion to obtain the upper bound

\begin{align*} f(\beta) & = f(\beta^*) + \int_0^1 \int_0^t (\beta - \beta^*)^\top H_{f}(\beta^* + s (\beta - \beta^*)) (\beta - \beta^*) \,\mathrm{d} s \,\mathrm{d} t \\ & \le f(\beta^*) + \tfrac{1}{2} (\beta - \beta^*)^\top ((1 + \varepsilon) H_0 I_d + \alpha C^{-1}) (\beta - \beta^*). \end{align*}

We then have a lower bound on the normalizing constant of the target posterior,

\begin{equation*} Z_{\Pi(\cdot \mid X, Y)} = \int_{\mathbb{R}^d} \exp\!({-}f(\beta)) \,\mathrm{d}\beta \ge \frac{\exp\!({-}f(\beta^*)) (2 \pi)^{d/2}}{\det\!((1 + \varepsilon) H_0 I_d + \alpha C^{-1})^{1/2}}. \end{equation*}

This in turn yields a lower bound on the ratio

\begin{align*} \frac{Z_{\Pi(\cdot \mid X, Y)}}{Z_Q} \exp\!(f(\beta^*)) \ge \frac{\det\!(\alpha C^{-1})^{1/2}}{\det\!((1 + \varepsilon) H_0 I_d + \alpha C^{-1})^{1/2}}. \end{align*}

The matrix C is symmetric and positive-definite, so its eigenvalues exist and are positive. Let $(\lambda_i(C))_{i = 1}^d$ be the eigenvalues of C. It is readily verified that the eigenvalues of the matrix $(1 + \varepsilon) H_0) I_d + \alpha C^{-1}$ exist and are $((1 + \varepsilon)H_0 + {\alpha}/{\lambda_i(C)})_{i = 1}^d$ . Then

(6) \begin{align} \frac{\det\!(\alpha C^{-1})}{\det\!((1 + \varepsilon) H_0 I_d + \alpha C^{-1})} & = \frac{\prod_{i = 1}^d {\alpha}/{\lambda_i(C)}} {\prod_{i = 1}^d((1 + \varepsilon) H_0 + {\alpha}/{\lambda_i(C)})} \nonumber \\ & = \prod_{i = 1}^d \frac{{\alpha}/{\lambda_i(C)}}{(1 + \varepsilon)H_0 + {\alpha}/{\lambda_i(C)}} \nonumber \\ & = \prod_{i = 1}^d \frac{1}{({1}/{\alpha}) (1 + \varepsilon) H_0 \lambda_i(C) + 1} \nonumber \\ & = \exp\Bigg({-}\sum_{i = 1}^d\log\bigg(\frac{1}{\alpha}(1 + \varepsilon)H_0\lambda_i(C) + 1\bigg)\Bigg). \end{align}

We have the basic inequality $\log\!(x + 1) \le x$ for any $x \in [0, +\infty)$ . Since the eigenvalues of C are positive and $H_0$ is nonnegative, we have the upper bound

\begin{equation*} \sum_{i = 1}^d \log\bigg(\frac{1}{\alpha} (1 + \varepsilon) H_0 \lambda_i(C) + 1 \bigg) \le \frac{1}{\alpha} (1 + \varepsilon) H_0 \sum_{i = 1}^d \lambda_i(C) \end{equation*}

This yields a lower bound on (6). Define the doubly-indexed sequence $(a_{d, n})$ by

\begin{equation*} a_{d, n} = \frac{1}{2\alpha} (1 + \varepsilon) H_0 \sum_{i = 1}^d \lambda_i(C). \end{equation*}

We have then shown that

\begin{equation*} \frac{Z_{\Pi(\cdot \mid X, Y)}}{Z_Q} \exp\!(f(\beta^*)) \ge \exp\!({-}a_{d, n}). \end{equation*}

By our assumption, $\mathrm{tr}(C) \to s_0$ as $d \to \infty$ . That is to say, $\lim_{d \to +\infty} \sum_{i = 1}^d \lambda_i(C) = s_0$ . This implies, by using continuity,

\begin{equation*} \lim_{d_n, n \to \infty}(1 - \exp\!({-}a_{n, d}))^t = (1 - \exp\!({-}(1 + \varepsilon) H_0 s_0 / (2 \alpha)))^t. \end{equation*}

By Corollary 1, we have the upper bound on the Wasserstein distance for each $d_n$ and each n:

\begin{align*} \mathcal{W}_{\rho}(P^t(\beta^*, \cdot), \Pi(\cdot \mid X, Y)) & = \bigg(1 - \exp\!(f(\beta^*)) \frac{Z_{\Pi(\cdot \mid X, Y)}}{Z_Q}\bigg)^t \int_{\mathbb{R}^d} \rho(\beta, \beta^*) \,\mathrm{d}\Pi(\beta \mid X, Y) \\ & \le (1 - \exp\!({-}a_{n, d}))^t \int_{\mathbb{R}^d} \rho(\beta, \beta^*) \,\mathrm{d}\Pi(\beta \mid X, Y). \end{align*}

Suppose that $\limsup_{d_n, n \to \infty}\int_{\mathbb{R}^d}\rho(\beta,\beta^*)\,\mathrm{d}\Pi(\beta\mid X,Y)<\infty$ . Using properties of the limit superior,

\begin{align*} \limsup_{d_n, n \to \infty} \mathcal{W}_{\rho}(P^t(\beta^*, \cdot), & \Pi(\cdot \mid X, Y)) \\ & \le \limsup_{d_n, n \to \infty} (1 - \exp\!({-}a_{n, d}))^t \limsup_{d_n, n \to \infty} \int_{\mathbb{R}^d} \rho(\beta, \beta^*) \,\mathrm{d}\Pi(\beta \mid X, Y) \\ & = \lim_{\substack{d, n \to \infty \\ {d}/{n} \to \gamma}} (1 - \exp\!({-}a_{n, d}))^t \limsup_{d_n, n \to \infty} \int_{\mathbb{R}^d} \rho(\beta, \beta^*) \,\mathrm{d}\Pi(\beta \mid X, Y) \\ & = (1 - \exp\!({-}(1 + \varepsilon) H_0 s_0 / (2 \alpha)))^t \limsup_{d_n, n \to \infty} \int_{\mathbb{R}^d} \rho(\beta, \beta^*) \,\mathrm{d}\Pi(\beta \mid X, Y). \end{align*}

This holds for every $\varepsilon \in (0, 1)$ , so taking the limit completes the proof in this case. The other case, when $\limsup_{d_n, n \to \infty}\int_{\mathbb{R}^d}\rho(\beta, \beta^*)\,\mathrm{d}\Pi(\beta\mid X, Y) = +\infty$ , is trivial.

Our goal now is to obtain an upper bound on the rate of convergence established in Corollary 1 in high dimensions for binary response regression. In this context, it is natural to treat the $(Y_i, X_i)_{i = 1}^n$ as stochastic; each time the sample size increases, the additional observation is randomly generated. Specifically, we assume that $(Y_i, X_i)_{i = 1}^n$ are independent with $Y_{i} \mid X_{i}, \beta \sim \text{Bernoulli}(S(\beta^\top X_i))$ and $X_i \sim N_{d}(0, \sigma^2 n^{-1} I_d)$ with $\sigma^2 \in (0, +\infty)$ . This scaling ensures the variance of the columns of the random matrix of features X is of fixed order; it is often used to ensure nondegeneracy in large-system limits of d, n [Reference Geman12]. Similar scaling assumptions on the data are used for high-dimensional maximum-likelihood theory in logistic regression [Reference Sur and Candès47].

Corollary 2. Suppose the following conditions hold:

1. (i) The negative log-likelihood $\ell_n$ is a twice continuously differentiable convex function.

2. (ii) There is a universal constant $r_0 \in (0, +\infty)$ such that the largest eigenvalue of the Hessian of the negative log-likelihood $H_{\ell_n}$ satisfies, for every $\beta \in \mathbb{R}^d$ , $\lambda_\mathrm{max}(H_{\ell_n}(\beta)) \le r_0 \lambda_\mathrm{max}(X^\top X)$ .

Let $a_0 = r_0 (1 + \gamma^{1/2})^2 \sigma^2 s_0 / (2 \alpha)$ . If $d, n \to +\infty$ in such a way that $d/n \to \gamma \in (0, +\infty)$ , then, almost surely, $\limsup_{d/n \to \gamma}\mathcal{W}_{\rho}(P^t(\beta^*, \cdot), \Pi(\cdot \mid X, Y)) \le M_0 (1 - \exp\!({-}a_0))^t$ , where $M_0 = \limsup_{d/n \to \gamma} \int \rho(\beta, \beta^*) \,\mathrm{d}\Pi( \beta \mid X, Y)$ .

Proof. Under our assumption, we may write the matrix $X = n^{-1/2} G$ , where G is a matrix with i.i.d. Gaussian entries with mean 0 and variance $\sigma^2$ . Denote the largest eigenvalue of the matrix $X^\top X$ by $\lambda_\mathrm{max}(X^\top X)$ . Therefore, as $d, n \to \infty$ in such a way that $d/n \to \gamma \in (0, +\infty)$ ,

\begin{equation*} \lambda_\mathrm{max}(X^\top X) = \lambda_\mathrm{max}\bigg(\frac{1}{n} G^\top G\bigg) = \frac{1}{n} \sup_{v \in \mathbb{R}^d, \left\lVert v\right\rVert_2 = 1} \left\lVert G^\top G v\right\rVert_2 \to (1 + \gamma^{1/2})^2\sigma^2 \end{equation*}

almost surely [Reference Geman12, Theorem 1]. We have, for any $\beta \in \mathbb{R}^d$ and any $v \in \mathbb{R}^d$ , $v^\top H_{\ell_n}(\beta) v \le r_0 \lambda_\mathrm{max}(X^\top X) \left\lVert v\right\rVert_2^2$ . The proof follows from Theorem 4.

Corollary 2 applies to both Bayesian logistic and probit regression. For logistic regression, $\ell_n$ is a twice continuously differentiable convex function and we may choose $r_0 = 4^{-1}$ . Similarly, for probit regression, $\ell_n$ is also a twice continuously differentiable convex function and we may choose $r_0 = 1$ [Reference Demidenko6].

In Figure 1 we plot $(1 - \exp\!({-}a_0))^t$ , the limiting decrease in the Wasserstein distance according to our upper bound, with varying values of the limiting ratio $\gamma$ and the other remaining values in $a_0$ fixed. We observe that as this ratio increases, the number of iterations needed to approximately converge may still increase rather rapidly.

Figure 1. The limiting decrease in the Wasserstein distance using different values of $\gamma$ , the limiting ratio of the dimension and sample size, versus the number of iterations.

6. Final remarks

We have studied the exact convergence behavior of the MHI sampler across general Wasserstein distances. We showed upper and lower bounds on the worst-case convergence rate for Wasserstein distances weaker than the total variation distance. We showed that the exact convergence rate at every initialization for Wasserstein distances weaker than the total variation distance is the same and matches that of the total variation convergence rate [Reference Wang52]. When starting at a certain point, we gave exact convergence expressions. By centering an independent proposal, we directly imbue the Markov chain with a strong attraction to a set where the target distribution has high probability. We showed this technique can provide uniform control over acceptance probability yielding exact convergence rates in Bayesian quantile regression. The centered MHI sampler turns out to have many applications for posteriors that arise in Bayesian generalized linear models where more sophisticated proposals are often used. With additional assumptions on the data and prior, we also showed that this exact convergence rate may be upper bounded when sampling high-dimensional posteriors in Bayesian binary response regression.