Dimension-free Relaxation Times of Informed MCMC Samplers on Discrete Spaces

Approximating probability distributions with unknown normalizing constants is a ubiquitous challenge in data science. In Bayesian statistics, the posterior probability distribution can be obtained through the Bayes rule, but without conjugacy, calculating its normalizing constant is a formidable challenge. In such scenarios, Markov chain Monte Carlo (MCMC) methods are commonly employed to simulate a Markov chain whose stationary distribution coincides with the target distribution. The convergence assessment of MCMC methods usually relies on an empirical diagnosis of Markov chain samples, such as effective sample size (Robert et al., 1999; Gong and Flegal, 2016) and Gelman-Rubin statistic (Gelman and Rubin, 1992), but the results can sometimes be misleading due to the intrinsic difficulty in detecting the non-convergence (Cowles and Carlin, 1996). For this reason, there has been a growing recognition of the importance of taking into account algorithmic convergence when evaluating Bayesian statistical models. A widely used metric for assessing the convergence of MCMC methods is the mixing time, which indicates the number of iterations required for the sampler to get sufficiently close to the target posterior distribution in total variation distance. Most existing studies on the mixing of MCMC algorithms consider Euclidean spaces. In particular, a very rich theory for log-concave target distributions has been developed (Dwivedi et al., 2019; Durmus and Moulines, 2019; Cheng et al., 2018; Dalalyan, 2017; Chewi et al., 2021), providing useful insights into the sampling complexity and guidance on the tuning of algorithm parameters. These results encompass a wide range of statistical models owing to the Bernstein-von Mises theorem (Belloni and Chernozhukov, 2009; Tang and Yang, 2022), which suggests that the posterior distribution becomes approximately normal (and thus log-concave) as the sample size tends to infinity. In contrast, the literature on the complexity of MCMC samplers for discrete statistical models is limited. Further, since every discrete space comes with its own combinatorial structure, researchers often choose to investigate each problem on a case-by-case basis rather than formulating a unified theory for arbitrary discrete spaces. Notable examples include variable selection (Yang et al., 2016), community detection (Zhuo and Gao, 2021), structure learning (Zhou and Chang, 2023), and classification and regression trees (CART) (Kim and Rockova, 2023).

Zhou and Chang (2023) undertook the endeavor to establish a general theoretical framework for studying the complexity of MCMC sampling for high-dimensional statistical models with discrete parameter spaces. They assumed a mild unimodal condition on the posterior distribution, akin to the log-concavity assumption in continuous-space problems, and derived mixing time bounds for random walk Metropolis-Hastings (MH) algorithms that were sharper than existing ones. Here “unimodal” means that other than the global maximizer, every state has a “neighbor” with strictly larger posterior probability, where two states are said to be neighbors if the proposal probability from one to the other is nonzero. Selecting this neighborhood relation often entails a trade-off, and the result of Zhou and Chang (2023) shows that rapid mixing can be achieved if (i) the neighborhood is large enough so that the posterior distribution becomes unimodal, and (ii) the neighborhood is not too large so that the proposal probability of each neighboring state is not exceedingly small. Techniques from high-dimensional statistical theory can be used to establish the unimodality in a way similar to how posterior consistency is proved, but the analysis is often highly complicated and tailored to individual problems. We will not explore these statistical techniques in depth in this paper; instead, our focus will be the mixing time analysis under the unimodal condition or a weaker condition that will be introduced later.

Our first objective is to extend the result of Zhou and Chang (2023) to more sophisticated MH schemes that do not use random walk proposals. We consider the informed MH algorithm proposed in Zanella (2020), which emulates the behavior of gradient-based MCMC samplers on Euclidean spaces (Duane et al., 1987; Roberts and Stramer, 2002; Girolami and Calderhead, 2011) and has gained increasing attention among the MCMC practitioners (Grathwohl et al., 2021; Zhang et al., 2022). The proposal schemes used in the informed algorithms always assess the local posterior landscape surrounding the current state and then tune the proposal probabilities to prevent the sampler from visiting states with low posterior probabilities. However, such informed schemes do not necessarily result in faster mixing, because the acceptance probabilities can be extremely low for proposed states, and it was shown in Zhou et al. (2022) that, for Bayesian variable selection, naive informed MH algorithms can even mix much more slowly than random walk MH algorithms. We will show that this acceptance probability issue can be overcome by generalizing the “thresholding” idea introduced in Zhou et al. (2022). Moreover, if the posterior distribution is unimodal and tails decay sufficiently fast, we prove that there always exists an informed MH scheme whose relaxation time can be bounded by a constant, independent of the problem dimension, which immediately leads to a nearly optimal bound on the mixing time. Our result significantly generalizes the finding of Zhou et al. (2022), which only considered the high-dimensional variable selection problem.

Delving into the more technical aspects, we investigate and compare two different approaches to obtaining sharp mixing time bounds on discrete spaces: the multicommodity flow method (Giné et al., 1996; Sinclair, 1992), and the single-element drift condition (Jerison, 2016). The former bounds the spectral gap of the transition matrix by identifying likely paths connecting any two distinct states. Compared to another path argument, known as “canonical path ensemble” and commonly used in the statistical literature (Yang et al., 2016; Zhuo and Gao, 2021; Kim and Rockova, 2023; Chang et al., 2022), the multicommodity flow method is more flexible and can yield tighter bounds. The drift condition method is based on a coupling argument and stands as the most popular technique for deriving the convergence rates of MCMC algorithms on general state spaces (Rosenthal, 1995; Roy and Hobert, 2007; Fort et al., 2003; Johndrow et al., 2020). A notable difference from the path method is that this approach directly bounds the total variation distance from the stationary distribution without analyzing the spectral gap. The existing literature suggests that both methods have their own unique strengths, and which method yields a better mixing time bound depends on the problem (Jerrum and Sinclair, 1996; Guruswami, 2000; Anil Kumar and Ramesh, 2001). To our knowledge, the two approaches have never been compared for analyzing MCMC algorithms for discrete-space statistical models. Indeed, the only works we are aware of that use drift conditions in these contexts are Zhou et al. (2022) and Kim and Rockova (2023), which considered variable selection and CART, respectively. We will demonstrate how the two methods can be applied under the general framework considered in this paper, and it will be shown that the drift condition approach yields a slightly sharper bound.

As the last major contribution of this work, we further extend our general theory beyond unimodal settings. For multimodal target distributions, while various mixing time bounds have been obtained (Guan and Krone, 2007; Woodard et al., 2009; Zhou and Smith, 2022), it is generally impossible to obtain rapid mixing results where the mixing time grows only polynomially with respect to the problem dimension. This is because the definition of mixing time considers the worst-case scenario regarding the choice of the initial distribution, and the chain can easily get stuck if it is initialized at a local mode. When these local modes (other than the global one) possess only negligible posterior mass and are unlikely to be visited by the sampler given a “warm” initialization, mixing time may provide an overly pessimistic estimate for the convergence rate. One possible remedy was described in the recent work of Atchadé (2021), who proposed to study initial-state-dependent mixing times using restricted spectral gap, a notion that generalizes the spectral gap of a transition matrix, and derived a rapid mixing result for Bayesian variable selection. We extend this technique to our setting and show that it can be integrated with the multicommodity flow method to produce sharp mixing time bounds for both random walk and informed MH schemes.

In Section 2, we review some recent results about the Bayesian variable selection problem, which will be used as an illustrative example throughout this paper; this section can be skipped for knowledgeable readers. Section 3 presents the setup for our theoretical analysis and reviews the mixing time bounds for random walk MH algorithms. In Section 4, we consider unimodal target distributions and prove mixing time bounds for informed MH algorithms via both the path method and drift condition analysis. Section 5 generalizes our results to a potentially multimodal setting. Simulation studies are conducted in Section 6, and Section 7 concludes the paper with some further discussion. A summary of the mixing time bounds obtained in this work is presented in Table 1.

Most technical details are deferred to the appendix. In Appendix A, we review the path method used in Zhou and Chang (2023) for bounding the spectral gaps of MH algorithms. In Appendix B, we review the result about restricted spectral gaps obtained in Atchadé (2021) and develop the multicommodity flow method for bounding restricted spectral gaps. Appendix C gives a brief review of the single-element drift condition of Jerison (2016). Proofs for all results presented in the main text are given in Appendix D.

We consider the standard linear regression model where a design matrix $\bm{X}\in\mathbb{R}^{n\times p}$ and a response vector $\bm{y}\in\mathbb{R}^{n}$ are assumed to satisfy

$\displaystyle\bm{y}=\bm{X}\bm{\beta}^{*}+\bm{z},\quad\bm{z}\sim\mathrm{MVN}(0,\sigma^{2}\bm{I}_{n}),$

(1)

where $\bm{\beta}^{*}\in\mathbb{R}^{p}$ is the unknown vector of regression coefficients and $\mathrm{MVN}$ denotes the multivariate normal distribution. Introduce the indicator vector $\delta=(\delta_{1},\dots,\delta_{p})\in\{0,1\}^{p}$ such that $\delta_{j}=1$ indicates that the $j$-th variable has a non-zero effect on the response. Variable selection is the task of identifying the true value of $\delta$, that is, finding the subset of variables with nonzero regression coefficients. Given $\delta$, we can write $\bm{y}=\bm{X}_{\delta}\bm{\beta}_{\delta}+\bm{z}$ with $\bm{z}\sim\mathrm{MVN}(0,\sigma^{2}\bm{I}_{n})$, where $\bm{X}_{\delta}$ and $\bm{\beta}_{\delta}$, respectively, denote the submatrix and subvector corresponding to those variables selected in $\delta$. We will also call $\delta\in\mathcal{V}$ a model, and $\delta=(0,0,\dots,0)$ will be referred to as the empty model. Consider the following two choices for the space of allowed models:

$$\mathcal{V}=\{0,1\}^{p},\text{ or }\mathcal{V}_{s}=\{\delta\in\mathcal{V}\colon||\delta||_{1}\leq s\},$$

(2)

where $||\cdot||_{1}$ denotes the $L^{1}$-norm and $s$ is a positive integer. The set $\mathcal{V}$ is the unrestricted space, while $\mathcal{V}_{s}$ represents a restricted space with some sparsity constraint. Let $|\cdot|$ denote the cardinality of a set. It is typically assumed in the high-dimensional literature that the sparsity parameter $s$ increases to infinity with $p$, in which case both $|\mathcal{V}|=2^{p}$ and $|\mathcal{V}_{s}|=O(p^{s})$ grow super-polynomially with $p$.

Spike-and-slab variable selection is a Bayesian procedure for constructing a posterior distribution over $\mathcal{V}$ or $\mathcal{V}_{s}$, which we denote by $\pi(\delta)$. In Section 6.1, we recall one standard approach to specifying the prior distribution for $(\delta,\bm{\beta},\sigma^{2})$, which leads to the following closed-form expression for $\pi$ up to a normalizing constant.

$$\quad\pi(\delta)\propto\frac{1}{p^{\kappa||\delta||_{1}}}\cdot\frac{(1+g)^{-(||\delta||_{1}/2)}}{\left(1+g\left(1-\mathsf{r}^{2}(\delta)\right)\right)^{n/2}}.$$

(3)

In (3), $\kappa,g$ are prior hyperparameters and $\mathsf{r}^{2}(\delta)=\bm{y}^{\top}(\bm{X}_{\delta}(\bm{X}_{\delta}^{\top}\bm{X}_{\delta})^{-1}\bm{X}_{\delta}^{\top})\bm{y}/(\bm{y}^{\top}\bm{y})$ is the coefficient of determination. Exact calculation of the normalizing constant is usually impossible, because it would require a summation over the entire parameter space, which involves super-polynomially many evaluations of $\pi$.

To find posterior probabilities of models of interest or evaluate integrals with respect to $\pi$, the most commonly used method is to use an MH algorithm to generate samples from $\pi$. The transition probability from $\delta$ to $\delta^{\prime}$ in an MH scheme can be expressed by

$\displaystyle\mathsf{P}(\delta,\delta^{\prime})=\mathsf{K}(\delta,\delta^{\prime})\mathsf{A}(\delta,\delta^{\prime}),$

where $\mathsf{K}(\delta,\delta^{\prime})$ is the probability of proposing to move from $\delta$ to $\delta^{\prime}$, and the associated acceptance probability $\mathsf{A}(\delta,\delta^{\prime})$ is given by

$$\mathsf{A}(\delta,\delta^{\prime})=\min\left\{1,\frac{\pi(\delta^{\prime})\mathsf{K}(\delta^{\prime},\delta)}{\pi(\delta)\mathsf{K}(\delta,\delta^{\prime})}\right\},$$

(4)

ensuring that $\pi$ is the stationary distribution of $\mathsf{P}$.

The efficiency of the MH algorithm depends on the choice of $\mathsf{K}$, and there are many strategies for selecting the proposal neighborhood and assigning the proposal probabilities. Recall that neighborhood refers to the support of the distribution $\mathsf{K}(\delta,\cdot)$ for each $\delta$, which we denote by $\mathcal{N}(\delta)$. Let us begin by considering random walk proposals such that $\mathsf{K}(\delta,\cdot)$ is simply the uniform distribution on $\mathcal{N}(\delta)$. To illustrate the importance of selecting a proper proposal neighborhood, we first present two naive choices that are bound to result in slow convergence.

The neighborhood size is exponential in $p$ in Example 1 and is a fixed constant in Example 2. Both are undesirable, and it is better to use a neighborhood with size polynomial in $p$. The following choice is common in practice:

$$\mathcal{N}_{1}(\delta)=\{\delta^{\prime}\in\mathcal{V}\colon||\delta-\delta^{\prime}||_{1}=1,\text{ and }\delta\neq\delta^{\prime}\}.$$

(5)

The set $\mathcal{N}_{1}(\delta)$ contains all the models that can be obtained from $\delta$ by either adding or removing a variable. When the design matrix $\bm{X}$ contains highly correlated variables, it is generally considered that $\mathcal{N}_{1}$ is too small and introducing “swap” moves is beneficial, which means to remove one variable and add another one at the same time. The resulting random walk MH algorithm is often known as the add-delete-swap sampler; denote its neighborhood by $\mathcal{N}_{\rm{ads}}$, which is defined by

$\displaystyle\mathcal{N}_{\rm{ads}}(\delta)=\mathcal{N}_{1}(\delta)\cup\mathcal{N}_{\mathrm{swap}}(\delta),\quad\mathcal{N}_{\mathrm{swap}}(\delta)=\left\{\delta^{\prime}\in\mathcal{V}\colon||\delta^{\prime}-\delta||_{1}=2,||\delta^{\prime}||_{1}=||\delta||_{1}\right\}.$

(6)

We will treat $\mathcal{N}_{\rm{ads}}$ and $\mathcal{N}_{1}$ as neighborhood relations defined on $\mathcal{V}$. When implementing the add-delete-swap sampler on the restricted space $\mathcal{V}_{s}$, one can still propose $\delta^{\prime}$ from $\mathcal{N}_{\rm{ads}}(\delta)$ and simply reject the proposal if $\delta^{\prime}\notin\mathcal{V}_{s}$. Note that $|\mathcal{N}_{\rm{ads}}(\delta)\cap\mathcal{V}_{s}|=O(ps)$.

We end this subsection with two comments. First, a popular alternative to MH algorithms is Gibbs sampling. Consider a random-scan Gibbs sampler that randomly picks $j\in\{1,2,\dots,p\}$ and updates $\delta_{j}$ from its conditional posterior distribution given the other coordinates. It is not difficult to see that this updating is equivalent to randomly proposing $\delta^{\prime}\in\mathcal{N}_{1}(\delta)$ and accepting $\delta^{\prime}$ with probability $\pi(\delta^{\prime})/(\pi(\delta)+\pi(\delta^{\prime}))$.¹¹1 In a general sense, this random-scan Gibbs sampler is also an MH scheme according to the original construction in Hastings (1970). By Peskun’s ordering (Mira, 2001), this Gibbs sampler is less efficient than the random walk MH algorithm with proposal neighborhood $\mathcal{N}_{1}$ (George and McCulloch, 1997).

Second, for simplicity, we assume in this work that $\mathsf{K}(\delta,\cdot)$ is a uniform distribution on $\mathcal{N}(\delta)$ for random walk proposals. But in practice, these proposals can be implemented in a more complicated, non-uniform fashion. For example, one can first choose randomly whether to add or remove a variable, and then given the type of move, a proposal of that type is generated with uniform probability. This distinction has minimal impact on all theoretical results we will develop.

The seminal work of Yang et al. (2016) considered a high-dimensional setting with $n,p,s\rightarrow\infty$ and $s\log p=o(n)$, and they proved that, under some mild assumptions, the mixing time of the add-delete-swap sampler on the restricted space $\mathcal{V}_{s}$ has order $O(pns^{2}\log p)$, polynomial in $(n,p,s)$; in other words, the sampler is rapidly mixing. To provide intuition about this result and proof techniques, which will be crucial to the understanding of the theory developed in this work, we construct a detailed illustrative example.

In Example 3, $\pi$ is unimodal with respect to $(\mathcal{V}_{2},\mathcal{N}_{\mathrm{ads}})$. An important intermediate result of Yang et al. (2016) was that, with high probability, such a unimodality property still holds in the high-dimensional regime they considered, and the global mode coincides with the true model $\delta^{*}$. This implies that the add-delete-swap sampler is unlikely to get stuck at any $\delta\neq\delta^{*}$, since there always exists some neighboring state $\delta^{\prime}$ such that $\mathsf{A}(\delta,\delta^{\prime})=1$ and thus $\mathsf{P}(\delta,\delta^{\prime})=\Omega(p^{-1}s^{-1})$. This is the main intuition behind the rapid mixing proof of Yang et al. (2016). Interestingly, it was shown in Yang et al. (2016) that the unimodality is unlikely to hold on the unrestricted space $\mathcal{V}$, since local modes can easily occur on the set $\mathcal{V}\setminus\mathcal{V}_{s}$, though it has negligible posterior mass. Since these local modes can easily trap the sampler for a huge number of iterations, they had to consider the restricted space $\mathcal{V}_{s}$ in order to obtain a rapid mixing result. They also found that without using swap moves, local modes are likely to occur at the boundary of the restricted space. (In Example 3, $(0,1,1)$ is such a local mode.) This is exactly the reason why swaps are required for their rapid mixing proof. More generally, for the MCMC convergence analysis of high-dimensional model selection problems with sparsity constraints, examining the local posterior landscape near the boundary seems often a major technical challenge (Zhou and Chang, 2023). Nevertheless, the practical implications of this theoretical difficulty is largely unclear. Even if we choose $\mathcal{V}$ as the search space, in practice, the sampler is typically initialized within $\mathcal{V}_{s}$, and if $s$ is large enough, we will probably never see the chain leave $\mathcal{V}_{s}$, and thus the landscape of $\pi$ on $\mathcal{V}\setminus\mathcal{V}_{s}$ is unimportant.

The above discussion naturally leads to the following question: assuming a warm initialization, can we obtain a “conditional” rapid mixing result on the unrestricted space $\mathcal{V}$? Here, “warm” means that the posterior probability of the model is not too small, and we use $\hat{\delta}$ to denote such an estimator, which can often be obtained by a frequentist variable selection algorithm. An affirmative result was given in Pollard and Yang (2019), where a random walk MH algorithm on $\mathcal{V}$ is constructed such that it proposes models from $\mathcal{N}_{1}(\delta)$ and is allowed to immediately jump back to $\hat{\delta}$ whenever $||\delta||_{1}$ becomes too large. They assumed $\hat{\delta}$ was obtained by the thresholded lasso method (Zhou, 2010). An even stronger result was obtained recently in Atchadé (2021), who showed that by only using addition and deletion proposals, the random walk MH algorithm on $\mathcal{V}$ is still rapidly mixing given a warm start. The proof of Atchadé (2021) utilizes a novel technique based on restricted spectral gap, which we will discuss in detail later.

Now consider informed proposal schemes (Zanella, 2020), which will be the focus of our theoretical analysis. Unlike random walk proposals which draw a candidate move indiscriminately from the given neighborhood, informed schemes compare the posterior probabilities of all models in the neighborhood and then assign larger proposal probabilities to those with larger posterior probabilities. One might conjecture that such an informed MH sampler may behave similarly to a greedy search algorithm and can quickly find the global mode if the target distribution is unimodal. But a naive informed proposal scheme may lead to performance even much worse than a random walk MH sampler, as illustrated in the following example.

In order to avoid scenario similar to Example 4, Zhou et al. (2022) proposed to use informed schemes with bounded proposal weights. They constructed an informed add-delete-swap sampler and showed that, under the high-dimensional setting considered by Yang et al. (2016), the mixing time on the restricted space is $O(n)$, thus independent of the problem dimension parameter $p$. Since each iteration of informed proposal requires evaluating the posterior probabilities of all models in the current neighborhood, the actual complexity of their algorithm should be $O(psn)$, comparable to the mixing time of the random walk MH sampler. The intuition behind the proof of Zhou et al. (2022), which was based on a drift condition argument, is similar to that of Yang et al. (2016). The unimodality of $\pi$ implies that any $\delta\neq\delta^{*}$ has one or more neighboring models with larger posterior probabilities, and they showed that, for their sampler, the transition probability to such a “good” neighbor is bounded from below by a universal constant.

In the remainder of this paper, we develop general theories that extend the results discussed in this section to arbitrary discrete spaces.

Let $\pi$ denote a probability distribution on a finite space $\mathcal{X}$. Let $\mathcal{X}$ be endowed with a neighborhood relation $\mathcal{N}\colon\mathcal{X}\rightarrow 2^{\mathcal{X}}$, and we say $x^{\prime}$ is a neighbor of $x$ if and only if $x^{\prime}\in\mathcal{N}(x)$. We assume $\mathcal{N}$ satisfies three conditions: (i) $x\notin\mathcal{N}(x)$ for each $x$, (ii) $x\in\mathcal{N}(x^{\prime})$ whenever $x^{\prime}\in\mathcal{N}(x)$, and (iii) $(\mathcal{X},\mathcal{N})$ is connected.²²2“Connected” means that for any $x,x^{\prime}\in\mathcal{X}$, there is a sequence $(x_{0}=x,x_{1},\dots,x_{k}=x^{\prime})$ such that $x_{i}\in\mathcal{N}(x_{i-1})$ for $i=1,\dots,k$. For our theoretical analysis, we will assume the triple $(\mathcal{X},\mathcal{N},\pi)$ is given and analyze the convergence of MH algorithms that propose states from $\mathcal{N}$. Define

$$M=M(\mathcal{X},\mathcal{N})\coloneqq\max_{x\in\mathcal{X}}|\mathcal{N}(x)|.$$

(8)

For most high-dimensional statistical problems, $|\mathcal{X}|$ grows at a super-polynomial rate with respect to some complexity parameter $p$, while $M$ only has a polynomial growth rate. Below are some examples.

1. (i)

   Variable selection. As discussed in Section 2, we can let $\mathcal{X}=\mathcal{V}$ and $\mathcal{N}=\mathcal{N}_{\mathrm{ads}}$, in which case we have $M(\mathcal{V},\mathcal{N}_{\rm{ads}})=O(p^{2})$. If we consider the restricted space $\mathcal{V}_{s}$, we have $M(\mathcal{V}_{s},\mathcal{N}_{\rm{ads}})=O(ps)$.

2. (ii)

   Structure learning of Bayesian networks. Given $p$ variables, one can let $\mathcal{X}$ be the collection of all $p$-node Bayesian networks (i.e., labeled directed acyclic graphs). Let $\mathcal{N}(x)$ be the collection of all Bayesian networks that can be obtained from $x$ by adding, deleting or reversing an edge (Madigan et al., 1995). We have $M=O(p^{2})$ since the graph has at most $O(p^{2})$ edges, while $|\mathcal{X}|$ is super-exponential in $p$ by Robinson’s formula (Robinson, 1977).

3. (iii)

   Ordering learning of Bayesian networks. Every Bayesian network is consistent with at least one total ordering of the $p$ nodes such that node $i$ precedes node $j$ whenever there is an edge from $i$ to $j$. To learn the ordering, we can let $\mathcal{X}$ be all possible orderings of $p$, which is the symmetric group with degree $p$. Clearly, $|\mathcal{X}|=p!$ is super-exponential in $p$. We may define $\mathcal{N}(x)$ as the set of all orderings that can be obtained from $x$ by a random transposition, which interchanges any two elements of $x$ while keeping the others unchanged (Chang et al., 2023; Diaconis and Shahshahani, 1981). This yields $M=p(p-1)/2$.

4. (iv)

   Community detection. Suppose there are $p$ nodes forming $K$ communities. One can let $\mathcal{X}$ be the collection of all label assignment vectors that specify the community label for each node, and define $\mathcal{N}(x)$ as the set of all assignments that differ from $x$ by the label of only one node (Zhuo and Gao, 2021). We have $|\mathcal{X}|=K^{p}$ and $M=p(K-1)$.

5. (v)

   Dyadic CART. Consider a classification or regression tree problem where the splits are selected from $p$ pre-specified locations. Assume $p=2^{K}-1$ for some integer $K\geq 1$, and consider the dyadic CART algorithm, a special case of CART, where splits always occur at midpoints (Donoho, 1997; Castillo and Ročková, 2021; Kim and Rockova, 2023). The search space $\mathcal{X}$ is the collection of all dyadic trees with depth less than or equal to $K$ (in a dyadic tree, every non-leaf node has 2 child nodes). One can show that $|\mathcal{X}|$ is exponential in $p$. For $x\in\mathcal{X}$, we define $\mathcal{N}(x)$ as the set of dyadic trees obtained by either a “grow” or “prune” operation; “grow” means to add two child nodes to one leaf node, and “prune” means to remove two leaf nodes with a common parent node. Then $M=O(p)$.

Let $\mathsf{P}$ be the transition matrix of an irreducible, aperiodic and reversible Markov chain with stationary distribution $\pi$. For all Markov chains we will analyze, $\mathsf{P}$ moves on the graph $(\mathcal{X},\mathcal{N})$; that is, $\{x^{\prime}\colon\mathsf{P}(x,x^{\prime})>0\}=\mathcal{N}(x)$. Define the total variation distance between $\pi$ and another distribution $\zeta$ by $||\zeta-\pi||_{\mathrm{TV}}=\sup_{A\subset\mathcal{X}}|\zeta(A)-\pi(A)|$. For $\epsilon\in(0,1/2)$, let

$$\tau_{x}(\mathsf{P},\epsilon)=\min\{t\in\mathbb{N}:||\mathsf{P}^{t}(x,\cdot)-\pi||_{\mathrm{TV}}\leq\epsilon\},$$

(9)

which can be seen as the “conditional” mixing time of $\mathsf{P}$ given the initial state $x$. Define the mixing time of $\mathsf{P}$ by

$$\tau(\mathsf{P},\epsilon)=\max_{x\in\mathcal{X}}\tau_{x}(\mathsf{P},\epsilon).$$

(10)

It is often assumed in the literature that $\epsilon=1/4$, because one can show that $\tau(\mathsf{P},\epsilon)\leq\lceil\log_{2}\epsilon^{-1}\rceil\tau(\mathsf{P},1/4)$ for any $\epsilon\in(0,1/2)$ (Levin and Peres, 2017, Chap. 4.5). However, since such an inequality does not hold for $\tau_{x}$, we will treat $\epsilon$ as an arbitrary constant in $(0,1/2)$ in this work. We can quantify the complexity of an MH algorithm by the product of its mixing time and the complexity per iteration.

It is well known that mixing time can be bounded by the spectral gap (Sinclair, 1992). Our assumption on $\mathsf{P}$ implies that it has real eigenvalues $1=\lambda_{1}>\dots\geq\lambda_{|\mathcal{X}|}>-1$ (Levin and Peres, 2017, Lemma 12.1). The spectral gap is defined as $\mathrm{Gap}(\mathsf{P})=1-\max\{\lambda_{2},|\lambda_{|\mathcal{X}|}|\}$ and satisfies the inequality

$$\tau(\mathsf{P},\epsilon)\leq\mathrm{Gap}(\mathsf{P})^{-1}\log\left\{\frac{1}{\epsilon\,\pi_{\min}}\right\},$$

(11)

where $\pi_{\min}=\min_{x\in\mathcal{X}}\pi(x)$. The quantity $\mathrm{Gap}(\mathsf{P})^{-1}$ is known as the relaxation time. In our mixing time bounds, we will always work with the “lazy” version of $\mathsf{P}$, which is defined by $\mathsf{P}^{\rm{lazy}}=(\mathsf{P}+\mathsf{I})/2$. That is, for any $x\neq x^{\prime}$, we set $\mathsf{P}^{\rm{lazy}}(x,x^{\prime})=\mathsf{P}(x,x^{\prime})/2$. Since all eigenvalues of $\mathsf{P}^{\rm{lazy}}$ are non-negative, $1-\mathrm{Gap}(\mathsf{P}^{\rm{lazy}})$ always equals the second largest eigenvalue of $\mathsf{P}^{\rm{lazy}}$.

To characterize the modality and tail behavior of $\pi$, we introduce another parameter $R$:

	$\displaystyle R=R(\mathcal{X},\mathcal{N},\pi)$	$\displaystyle\coloneqq\min_{x\in\mathcal{X}\setminus\{x^{*}\}}\max_{y\in\mathcal{N}(x)}\frac{\pi(y)}{\pi(x)},$		(12)
	$\displaystyle\text{ where }x^{*}$	$\displaystyle=\operatorname*{arg\,max}_{x\in\mathcal{X}}\pi(x).$		(13)

If $\operatorname*{arg\,max}$ in the above equation is not uniquely defined, fix $x^{*}$ as one of the modes according to any pre-specified rule. If $R>1$, we say $\pi$ is unimodal with respect to $\mathcal{N}$, since for any $x\neq x^{*}$, there is some neighboring state $x^{\prime}\in\mathcal{N}(x)$ such that $\pi(x^{\prime})>\pi(x)$. In our subsequent analysis of MH algorithms, we will consider unimodal targets with $R>M$ (for random walk MH) or $R>M^{2}$ (for informed MH). These conditions ensure that the tails of $\pi$ decay sufficiently fast. To see this, for $k\geq 1$, let $\mathrm{Tail}(k)=\{x\colon\mathsf{P}^{k-1}(x,x^{*})=0,\;\mathsf{P}^{k}(x,x^{*})>0\}$ denote the set of states that need at least $k$ steps to reach $x^{*}$. Then, $\pi(\mathrm{Tail}(k))\leq(M/R)^{k}$, which decreases to $0$ exponentially fast if $R>M$.

As discussed in Section 2, for high-dimensional variable selection, the unimodal property has been rigorously established on $\mathcal{X}=\mathcal{V}_{s}$ with $\mathcal{N}=\mathcal{N}_{\mathrm{ads}}$. A careful examination of the proof of Yang et al. (2016) reveals that, since $M$ is polynomial in $p$, the proof for $R>1$ and that for $R>M^{2}$ are essentially the same; see the discussion in Section S3 of Zhou et al. (2022). The unimodal condition has also been proved for other high-dimensional statistical problems, including structure learning of Markov equivalence classes (Zhou and Chang, 2023), community detection (Zhuo and Gao, 2021), and dyadic CART (Kim and Rockova, 2023). It should be noted that how to choose a proper $\mathcal{N}$ so that $R>M$ or $R>M^{2}$ can be a very challenging question (e.g., for structure learning), which we do not elaborate on in this paper.

In order to achieve rapid mixing, such unimodal conditions are arguably necessary, since for general multimodal targets, the chain can get trapped at some local mode for an arbitrarily large amount of time. We further highlight two reasons why the unimodal analysis is important and not as restrictive as it may seem. First, mixing time bounds for unimodal targets are often the building blocks for more general results in multimodal scenarios. One strategy that will be discussed shortly is to use restricted spectral gap. Another approach is to apply state decomposition techniques (Madras and Randall, 2002; Jerrum et al., 2004; Guan and Krone, 2007), which works for general multimodal targets; see, e.g., Zhou and Smith (2022). Second, as explained in Remark 4 of Zhou and Smith (2022), the condition $R>M$ is very similar to log-concavity on Euclidean spaces, since it implies unimodality and exponentially decaying tails. Moreover, as we have illustrated in Example 3, for a variable selection problem with fixed $p$ and sufficiently large $n$, $\pi$ is typically unimodal, but $\pi(\delta)$ may not increase as $\delta$ gets closer to the mode in $L^{1}$ distance. This suggests that our unimodal condition is conceptually more general than log-concavity, since the latter also implies unimodality in a slice taken along any direction.

For the analysis beyond the unimodal setting, we consider a subset $\mathcal{X}_{0}\subset\mathcal{X}$ and only impose unimodality on $\mathcal{X}_{0}$. Let $\mathcal{N}|_{\mathcal{X}_{0}}$ denote $\mathcal{N}$ restricted to $\mathcal{X}_{0}$, that is, $\mathcal{N}|_{\mathcal{X}_{0}}(x)=\mathcal{N}(x)\cap\mathcal{X}_{0}$. Mimicking the definition of $R$, we define $R|_{\mathcal{X}_{0}}$ by restricting ourselves to $\mathcal{X}_{0}$:

	$\displaystyle R|_{\mathcal{X}_{0}}=R(\mathcal{X}_{0},\mathcal{N}|_{\mathcal{X}_{0}},\pi)$	$\displaystyle=\min_{x\in\mathcal{X}_{0}\setminus\{x^{*}_{0}\}}\max_{x^{\prime}\in\mathcal{N}|_{\mathcal{X}_{0}}(x)}\frac{\pi(x^{\prime})}{\pi(x)},$		(14)
	$\displaystyle\text{ where }x^{*}_{0}$	$\displaystyle=\operatorname*{arg\,max}_{x\in\mathcal{X}_{0}}\pi(x).$		(15)

If $R|_{\mathcal{X}_{0}}>1$, we say $\pi$ is unimodal on $\mathcal{X}_{0}$ with respect to $\mathcal{N}$. Note that $R|_{\mathcal{X}_{0}}>1$ also implies that $(\mathcal{X}_{0},\mathcal{N})$ is connected. In Section 5, we will study the mixing times of MH algorithms assuming $R|_{\mathcal{X}_{0}}>M$ (or $R|_{\mathcal{X}_{0}}>M^{2}$) for some $\mathcal{X}_{0}$ such that $\pi(\mathcal{X}_{0})$ is sufficiently large.

We first review the mixing time bound for random walk algorithms obtained in Zhou and Chang (2023) under a unimodal condition. Recall that we assume the random walk proposal scheme can be expressed as

$$\mathsf{K}(x,x^{\prime})=|\mathcal{N}(x)|^{-1}\mathbbm{1}_{\mathcal{N}(x)}(x^{\prime}).$$

The resulting transition matrix of random walk MH can be written as

$\displaystyle\mathsf{P}_{0}(x,x^{\prime})=\left\{\begin{array}[]{cc}\min\left\{\frac{1}{|\mathcal{N}(x)|},\,\frac{\pi(x^{\prime})}{\pi(x)|\mathcal{N}(x^{\prime})|}\right\},&\text{ if }x^{\prime}\in\mathcal{N}(x),\vskip 3.0pt plus 1.0pt minus 1.0pt\\ 1-\sum_{\tilde{x}\in\mathcal{N}(x)}\mathsf{P}_{0}(x,\tilde{x}),&\text{ if }x^{\prime}=x,\vskip 3.0pt plus 1.0pt minus 1.0pt\\ 0,&\text{ otherwise. }\end{array}\right.$

(19)

Yang et al. (2016) used (11) and the “canonical path ensemble” argument to bound the mixing time of the add-delete-swap sampler for high-dimensional variable selection. As observed in Zhou and Smith (2022) and Zhou and Chang (2023), the method of Yang et al. (2016) is applicable to the general setting we consider. The only assumption one needs is that the triple $(\mathcal{X},\mathcal{N},\pi)$ satisfies $R>M$, which ensures that, for any $x\neq x^{*}$, we can identify a path $(x_{0}=x,x_{1},\dots,x_{k-1},x_{k}=x^{*})$ such that $\pi(x_{j})/\pi(x_{j-1})>M$ for each $j$. A canonical path ensemble is a collection of such paths, one for each $x\neq x^{*}$, and then a spectral gap bound can be obtained by identifying the maximum length of a canonical path and the edge subject to the most congestion. It was shown in Zhou and Smith (2022) that the bound of Yang et al. (2016) can be further improved by measuring the length of each path using a metric depending on $\pi$ (instead of counting the number of edges). The following result is a direct consequence of (11) and Lemma 3 of Zhou and Smith (2022).

Informed MH algorithms generate proposal moves after evaluating the posterior probabilities of all neighboring states. Given a weighting function $h:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$, we define the informed proposal by

$\displaystyle\mathsf{K}_{h}\left(x,x^{\prime}\right)=\frac{h\left(\frac{\pi\left(x^{\prime}\right)}{\pi(x)}\right)}{Z_{h}(x)}\mathbbm{1}_{\mathcal{N}(x)}\left(x^{\prime}\right),\text{ where }Z_{h}(x)=\sum_{\tilde{x}\in\mathcal{N}(x)}h\left(\frac{\pi(\tilde{x})}{\pi(x)}\right);$

(22)

that is, $\mathsf{K}_{h}(x,\cdot)$ draws $x^{\prime}$ from the set $\mathcal{N}(x)$ with probability proportional to $h(\pi(x^{\prime})/\pi(x))$. Intuitively, one wants to let $h$ be non-decreasing so that states with larger posterior probabilities receive larger proposal probabilities. Choices such as $h(u)=1+u,\sqrt{u}$ or $1\wedge u$ were analyzed in Zanella (2020). It was observed in Zhou et al. (2022) and also illustrated by our Example 4 that for problems such as high-dimensional variable selection, such choices could be problematic and lead to even worse mixing than random walk MH algorithms. To gain a deeper insight, we write down the transition matrix of the induced MH algorithm:

$\displaystyle\mathsf{P}_{h}\left(x,x^{\prime}\right)=\left\{\begin{array}[]{cc}\mathsf{K}_{h}\left(x,x^{\prime}\right)\min\left\{1,\frac{\pi\left(x^{\prime}\right)\mathsf{K}_{h}\left(x^{\prime},x\right)}{\pi(x)\mathsf{K}_{h}\left(x,x^{\prime}\right)}\right\},&\text{ if }x^{\prime}\neq x,\\ 1-\sum_{\tilde{x}\neq x}\mathsf{P}_{h}(x,\tilde{x}),&\text{ if }x^{\prime}=x.\end{array}\right.$

(23)

We expect that $\mathsf{K}_{h}(x,\cdot)$ proposes with high probability some $x^{\prime}$ such that $\pi(x^{\prime})\gg\pi(x)$. But $\mathsf{K}_{h}(x^{\prime},x)$, which depends on the local landscape of $\pi$ on $\mathcal{N}(x^{\prime})$, can be arbitrarily small if $h$ is unbounded, causing the acceptance probability of the proposal move from $x$ to $x^{\prime}$ to be exceedingly small.

The solution proposed in Zhou et al. (2022) was to use some $h$ that is bounded both from above and from below. In this work, we consider the following choice of $h$:

$$h(u)=\mathrm{clip}(u,\ell,L)\coloneqq\left\{\begin{array}[]{cc}\ell,&\text{ if }u<\ell,\\ u,&\quad\text{ if }\ell\leq u\leq L,\\ L,&\text{ if }u>L,\end{array}\right.$$

(24)

where $\ell<L$ are some constants. Henceforth, whenever we write $\mathsf{K}_{h}$ or $\mathsf{P}_{h}$, it is understood that $h$ takes the form given in (24). We now revisit Example 3 and show that this bounded weighting scheme overcomes the issue of diminishing acceptance probabilities.

What we have observed in Example 5 is not a coincidence. The following lemma gives simple conditions under which an informed proposal is guaranteed to have acceptance probability equal to $1$.

The assumption $Z_{h}(x)\geq L$ used in Lemma 1 is weak, and it will be satisfied if $x$ has one neighboring state $z$ such that $\pi(z)/\pi(x)\geq L$, which is true if $L\leq R$.

Using Lemma 1 and Theorem 2 of Zhou and Chang (2023), we can now prove a sharp mixing time bound for informed MH algorithms under the assumption $R>M^{2}$.

The proof of Theorem 2 utilizes the multicommodity flow method (Sinclair, 1992), which generalizes the canonical path ensemble argument and allows us to select multiple likely transition routes between any two states. If one uses the canonical path ensemble, the resulting relaxation time bound for informed MH algorithms will still involve a factor of $M$ as in Theorem 1. A detailed review of the multicommodity flow method is given in Appendix A.

The mixing time analysis we have carried out so far relies on employing path methods to establish bounds on the spectral gap—an approach that is predominant in the literature on finite-state Markov chains. One exception was the recent work of Zhou et al. (2022), who used drift condition to study the mixing time of an informed add-delete-swap sampler for Bayesian variable selection. In this section, we show that their method can also be applied to general discrete spaces, and it can be used to improve the mixing time bound obtained in Theorem 2 in an asymptotic setting.

We still let $x^{*}=\operatorname*{arg\,max}_{x\in\mathcal{X}}\pi(x)$. The strategy of the proof is to establish a drift condition on $\mathcal{X}\setminus\{x^{*}\}$, which means that for some $V\colon\mathcal{X}\rightarrow[1,\infty)$ and $\alpha\in(0,1)$,

$$(\mathsf{P}_{h}V)(x)\leq\alpha V(x),\quad\text{ for any }x\in\mathcal{X}\setminus\{x^{*}\},$$

(26)

where $(\mathsf{P}_{h}V)(x)=\sum_{x^{\prime}\in\mathcal{X}}\mathsf{P}_{h}(x,x^{\prime})V(x^{\prime}).$ Then, one can apply Theorem 4.5 of Jerison (2016) to obtain a mixing time bound of the order $(1-\alpha)^{-1}\log V$. The main challenge is to find an appropriate $V$ such that $\alpha$ is as small as possible. After several trials, we find that the choice,

$$V(x)=\pi(x)^{1/\log\pi_{\mathrm{min}}}=\exp\left(\frac{\log\pi(x)}{\log\pi_{\mathrm{min}}}\right),$$

(27)

yields the desired mixing rate given in Theorem 3 below. This drift function is similar to but simpler than the one used in Zhou et al. (2022) for variable selection. By definition, $V(x)$ always decreases as $\pi(x)$ increases since $\pi_{\mathrm{min}}<1$, and $V$ is bounded on $[1,e]$.