Supplement to “Markov Chain Monte Carlo Based on Deterministic Transformations”

S-1 Proof of detailed balance for TMCMC

The detailed balance condition is proved as follows: Suppose $\mathbf{y}=T_{\mathbf{z}}(\mathbf{x},\boldsymbol{\epsilon})\in T_{\mathbf{z}}(\mathbf{x},\mathcal{Y})$ , then $\mathbf{x}=T_{\mathbf{z}^{c}}(\mathbf{y},\boldsymbol{\epsilon})$ . Hence, the kernel $K$ satisfies,

	$\displaystyle\pi(\mathbf{x})K(\mathbf{x}\to\mathbf{y})$	$\displaystyle=$	$\displaystyle\pi(\mathbf{x})~P(T_{\mathbf{z}})~g(\boldsymbol{\epsilon})\min\left\{1,\dfrac{P(T_{\mathbf{z}^{c}})\pi(\mathbf{y})}{P(T_{\mathbf{z}})\pi(\mathbf{x})}J_{\mathbf{z}}(\mathbf{x},\boldsymbol{\epsilon})\right\}$
		$\displaystyle=$	$\displaystyle g(\boldsymbol{\epsilon})~\min\left\{\pi(\mathbf{x})~P(T_{\mathbf{z}}),\pi(\mathbf{y})P(T_{\mathbf{z}^{c}})~J_{\mathbf{z}}(\mathbf{x},\boldsymbol{\epsilon})\right\}$

and

	$\displaystyle\pi(\mathbf{y})K(\mathbf{y}\to\mathbf{x})$	$\displaystyle=$	$\displaystyle\pi(\mathbf{y})~P(T_{\mathbf{z}^{c}})~g(\boldsymbol{\epsilon})J_{\mathbf{z}}(\mathbf{x},\boldsymbol{\epsilon})\min\left\{1,\dfrac{P(T_{\mathbf{z}})\pi(\mathbf{x})}{P(T_{\mathbf{z}^{c}})\pi(\mathbf{y})}J_{\mathbf{z}^{c}}(\mathbf{y},\boldsymbol{\epsilon})~\right\}$
		$\displaystyle=$	$\displaystyle g(\boldsymbol{\epsilon})~\min\left\{\pi(\mathbf{y})~P(T_{\mathbf{z}^{c}})J_{\mathbf{z}}(\mathbf{x},\boldsymbol{\epsilon}),\pi(\mathbf{x})P(T_{\mathbf{z}})\right\}$

where $J_{\mathbf{z}}(\mathbf{x},\boldsymbol{\epsilon})=\left|\partial(T_{\mathbf{z}}(\mathbf{x},\boldsymbol{\epsilon}),\boldsymbol{\epsilon})/\partial(\mathbf{x},\boldsymbol{\epsilon})\right|$ satisfies

J_{\mathbf{z}^{c}}(T_{\mathbf{z}}(\mathbf{x},\boldsymbol{\epsilon}),\boldsymbol{\epsilon})\times J_{\mathbf{z}}(\mathbf{x},\boldsymbol{\epsilon})=1\quad\textrm{ since }\quad T_{\mathbf{z}^{c}}(T_{\mathbf{z}}(\mathbf{x},\boldsymbol{\epsilon}),\boldsymbol{\epsilon})=\mathbf{x}.

S-2 General TMCMC algorithm based on a single $\epsilon$

Algorithm S-2.1

General TMCMC algorithm based on a single $\epsilon$ .

•

Input: Initial value ${\mathbf{x}}^{(0)}$ , and number of iterations $N$ .

•

For $t=0,\ldots,N-1$

1.

Generate $\epsilon\sim g(\cdot)$ and an index $i\sim\mathcal{M}(1;p_{1},\ldots,p_{3^{k}-1})$ independently. Again, actual simulation from the high-dimensional multinomial distribution is not necessary; see Section 3.1 of DB.

2.

\mathbf{x}^{\prime}=T_{\mathbf{z}_{i}}({\mathbf{x}}^{(t)},\epsilon)\quad\textrm{ and }\quad\alpha({\mathbf{x}}^{(t)},\epsilon)=\min\left(1,\dfrac{P(T_{\mathbf{z}^{c}_{i}})}{P(T_{\mathbf{z}_{i}})}~\dfrac{\pi(\mathbf{x}^{\prime})}{\pi({\mathbf{x}}^{(t)})}~\left|\frac{\partial(T_{\mathbf{z}_{i}}({\mathbf{x}}^{(t)},\epsilon),\epsilon)}{\partial({\mathbf{x}}^{(t)},\epsilon)}\right|\right)

3.

Set

{\mathbf{x}}^{(t+1)}=\left\{\begin{array}[]{ccc}\mathbf{x}^{\prime}&\textsf{ with probability }&\alpha({\mathbf{x}}^{(t)},\epsilon)\\ {\mathbf{x}}^{(t)}&\textsf{ with probability }&1-\alpha({\mathbf{x}}^{(t)},\epsilon)\end{array}\right.

•

End for

S-3 Convergence properties of additive TMCMC

In this section we prove some convergence properties of the TMCMC in the case of the additive transformation. Before going into our main result we first borrow some definitions from the MCMC literature.

Definition 1 (Irreducibility)

A Markov transition kernel $K$ is $\varphi-$ irreducible, where $\varphi$ is a nontrivial measure, if for every $x\in\mathcal{X}$ and for every measurable set $A$ of $\mathcal{X}$ with $\varphi(A)>0$ , there exists $n\in\mathbb{N}$ , such that $K^{n}(x,A)>0.$

Definition 2 (Small set)

A measurable subset $E$ of $\mathcal{X}$ is said to be small if there is an $n\in\mathbb{N}$ , a constant $c>0$ , possibly depending on $E$ and a finite measure $\nu$ such that

K^{n}(x,A)\geq c~\nu(A),\qquad\forall~A\in\mathcal{B}(\mathcal{X}),~\forall~x\in E

Definition 3 (Aperiodicity)

A Markov kernel $K$ is said to be periodic with period $d>0$ if the state-space $\mathcal{X}$ can be partitioned into $d$ disjoint subsets $\mathcal{X}_{1},\mathcal{X}_{2},\ldots,\mathcal{X}_{d}$ with

K(x,\mathcal{X}_{i+1})=1~\forall~x\in\mathcal{X}_{i},~i=1,2,\ldots,d-1

and $K(x,\mathcal{X}_{1})=1~\forall~x\in\mathcal{X}_{d}$ .

A Markov kernel $K$ is aperiodic if for no $d\in\mathbb{N}$ it is periodic with period $d$ .

S-3.1 Additive transformation with singleton $\epsilon$

Consider now the case where $\mathcal{X}=\mathbb{R}^{k}$ , $\mathcal{D}=\mathbb{R}$ and $T_{\mathbf{z}}(\mathbf{x},\epsilon)=(x_{1}+z_{1}a_{1}\epsilon,x_{2}+z_{2}a_{2}\epsilon,\ldots,x_{k}+z_{k}a_{k}\epsilon)$ where, for $i=1,\ldots,k$ , $z_{i}=\pm 1$ , and $a_{i}>0$ . In this case $\mathcal{Y}=[0,\infty)$ . Suppose that $g$ is a density on $\mathcal{Y}$ .

Theorem 1

Suppose that $\pi$ is bounded and positive on every compact subset of $\mathbb{R}^{k}$ and that $g$ is positive on every compact subset of $(0,\infty)$ . Then the chain is $\lambda$ -irreducible, aperiodic. Moreover every nonempty compact subset of $\mathbb{R}^{k}$ is small.

Proof 1

Without loss we may assume that $a_{i}=1;$ $i=1,\ldots,k$ . For notational convenience we shall prove the theorem for $k=2$ . The general case can be seen to hold with suitably defined ‘rotational’ matrices on $\mathbb{R}^{k}$ similar to (S-3.1).

Suppose $E$ is a nonempty compact subset of $\mathbb{R}^{k}$ . Let $C$ be a compact rectangle whose sides are parallel to the diagonals $\{(x,y)~:~|y|=|x|\}$ and containing $E$ such that $\lambda(C)>0$ . We shall show that $E$ is small, i.e., $\exists~c>0$ such that

K^{2}(\mathbf{x},A)\geq c\lambda_{C}(A)\qquad\forall A\in\mathcal{B}(\mathbb{R}^{2})\textrm{ and }\forall x\in E.

It is clear that the points reachable from $\mathbf{x}$ in two steps are of the form

\begin{pmatrix}x_{1}\pm\epsilon_{1}\pm\epsilon_{2}\\ x_{2}\pm\epsilon_{1}\pm\epsilon_{2}\end{pmatrix},\qquad\epsilon_{1}\geq 0,\epsilon_{2}\geq 0

Thus, if we define the matrices

\begin{split}&M_{1}=\begin{pmatrix}1&1\\ 1&-1\end{pmatrix}\quad M_{2}=\begin{pmatrix}-1&1\\ 1&1\end{pmatrix}\quad M_{3}=\begin{pmatrix}1&-1\\ -1&-1\end{pmatrix}\quad M_{4}=\begin{pmatrix}-1&-1\\ -1&1\end{pmatrix}\\ &\tilde{M}_{1}=\begin{pmatrix}1&1\\ -1&1\end{pmatrix}\quad\tilde{M}_{2}=\begin{pmatrix}1&-1\\ 1&1\end{pmatrix}\quad\tilde{M}_{3}=\begin{pmatrix}-1&1\\ -1&-1\end{pmatrix}\quad\tilde{M}_{4}=\begin{pmatrix}-1&-1\\ 1&-1\end{pmatrix}\end{split}

(S-3.1)

then the points reachable from $\mathbf{x}$ in two steps, other than the points lying on the diagonals passing through $\mathbf{x}$ itself, are of the form

\mathbf{x}+M_{i}\left(\begin{smallmatrix}\epsilon_{1}\\ \epsilon_{2}\end{smallmatrix}\right)\quad\textrm{ and }\quad\mathbf{x}+\tilde{M}_{i}\left(\begin{smallmatrix}\epsilon_{1}\\ \epsilon_{2}\end{smallmatrix}\right),\quad\epsilon_{1}>0,\epsilon_{2}>0,~i=1,\ldots,4.

Define

m=\inf_{\mathbf{y}\in C}\pi(\mathbf{y})>0\qquad M=\sup_{\mathbf{y}\in C}\pi(\mathbf{y})<\infty\qquad a=\inf_{0<\epsilon<R}g(\epsilon)>0

where $R$ is the length of the diagonal of the rectangle $C$ ¹¹1Actually $R/\sqrt{2}$ suffices.. Fix an element $\mathbf{x}\in E$ . For any set $A\in\mathcal{B}(\mathbb{R}^{2})$ , let $A^{*}=A\cap C$ and define,

\begin{split}A_{i}&=\{\boldsymbol{\epsilon}\in(0,\infty)^{2}~:~\mathbf{x}+M_{i}\boldsymbol{\epsilon}\in A^{*}\}\\ \tilde{A}_{i}&=\{\boldsymbol{\epsilon}\in(0,\infty)^{2}~:~\mathbf{x}+\tilde{M}_{i}\boldsymbol{\epsilon}\in A^{*}\}\end{split}

(S-3.2)

The need for defining such sets illustrated in the following example: to make a transition from the state $\mathbf{x}$ to a state in $A^{*}$ in two steps, first making a forward transition in both coordinates and then a forward transition in first coordinate and a backward transition in the second coordinate is same as applying the transformation $\mathbf{x}\to\mathbf{x}+M_{1}\boldsymbol{\epsilon}$ for some $\boldsymbol{\epsilon}\in A_{1}$ in two steps, i.e. first

\mathbf{x}\to\mathbf{x}+M_{1}(\epsilon_{1},0)^{T}=\mathbf{x}+(\epsilon_{1},\epsilon_{1})^{T}\quad\textrm{ then }\quad\mathbf{x}+M_{1}(\epsilon_{1},\epsilon_{2})^{T}\to\mathbf{x}+M_{1}\boldsymbol{\epsilon}

Also note that for any $\boldsymbol{\epsilon}=(\epsilon_{1},\epsilon_{2})\in A_{i}$ , $A^{*}\subset C$ implies that the intermediate point $\mathbf{x}+M_{i}(\epsilon_{1},0)^{T}\in C$ and similarly for $\tilde{A}_{i}~(i=1,\ldots,4)$ . Now, with ${\underline{p}}$ and $\bar{p}$ as the minimum and maximum of the move probabilities.

	$\displaystyle K^{2}(\mathbf{x},A)\geq K^{2}(\mathbf{x},A^{*})$	(S-3.3)
$\displaystyle\geq$	$\displaystyle{\underline{p}}^{2}\sum_{i=1}^{4}\int\displaylimits_{A_{i}}g(\epsilon_{1})g(\epsilon_{2})\min\left\{\frac{{\underline{p}}\pi(\mathbf{x}+M_{i}(\epsilon_{1},0)^{T})}{{\bar{p}}\pi(\mathbf{x})},1\right\}\min\left\{\frac{{\underline{p}}\pi(\mathbf{x}+M_{i}(\epsilon_{1},\epsilon_{2})^{T})}{{\bar{p}}\pi(\mathbf{x}+M_{i}(\epsilon_{1},0)^{T})},1\right\}d\epsilon_{1}d\epsilon_{2}$
$\displaystyle+$	$\displaystyle{\underline{p}}^{2}\sum_{i=1}^{4}\int\displaylimits_{\tilde{A}_{i}}g(\epsilon_{1})g(\epsilon_{2})\min\left\{\frac{{\underline{p}}\pi(\mathbf{x}+\tilde{M}_{i}(\epsilon_{1},0)^{T})}{{\bar{p}}\pi(\mathbf{x})},1\right\}\min\left\{\frac{{\underline{p}}\pi(\mathbf{x}+\tilde{M}_{i}(\epsilon_{1},\epsilon_{2})^{T})}{{\bar{p}}\pi(\mathbf{x}+\tilde{M}_{i}(\epsilon_{1},0)^{T})},1\right\}d\epsilon_{1}d\epsilon_{2}$
$\displaystyle\geq$	$\displaystyle{\underline{p}}^{2}a^{2}\left(\min\left\{\dfrac{\underline{p}m}{\bar{p}M},1\right\}\right)^{2}\left(\sum_{i=1}^{4}\lambda(A_{i})+\sum_{i=1}^{4}\lambda(\tilde{A}_{i})\right)$
$\displaystyle=$	$\displaystyle{\underline{p}}^{2}a^{2}\left(\min\left\{\dfrac{\underline{p}m}{\bar{p}M},1\right\}\right)^{2}\times 2\times\sum_{i=1}^{4}\lambda(A_{i})$

Since $(\epsilon_{1},\epsilon_{2})\in A_{i}\iff(\epsilon_{2},\epsilon_{1})\in\tilde{A}_{i}$ , so that, $\lambda(A_{i})=\lambda(\tilde{A}_{i})$ . Now notice that, if we define for $i=1,\ldots,4$

f_{i}:(0,\infty)^{2}\to\mathbb{R}^{2}\ni\boldsymbol{\epsilon}\mapsto\mathbf{x}+M_{i}\boldsymbol{\epsilon}

and

A_{\mathbf{x}}=\{(\epsilon,0)^{T}~:~\epsilon>0,(x_{1}\pm\epsilon,x_{2}\pm\epsilon)\in A^{*}\}

then,

A^{*}=\bigcup_{i=1}^{4}f_{i}(A_{i}\cup A_{x})\quad\Longrightarrow\quad\lambda(A^{*})=\sum_{i=1}^{4}f_{i}(A_{i})\quad=\quad 2\times\sum_{i=1}^{4}\lambda(A_{i}),

since, $f_{i}(A_{i})$ ’s are pairwise disjoint, $\lambda(f_{i}(A_{x}))=0$ and $\lambda(f_{i}(A_{i}))=2\lambda(A_{i})$ for $1\leq i\leq 4$ . It follows from (S-3.3) that

K^{2}(\mathbf{x},A)\quad\geq\quad{\underline{p}}^{2}a^{2}\left(\min\left\{\dfrac{\underline{p}m}{\bar{p}M},1\right\}\right)^{2}\lambda(A^{*})\quad=\quad c\lambda_{C}(A)

where $c={\underline{p}}^{2}a^{2}\left(\min\left\{\dfrac{\underline{p}m}{\bar{p}M},1\right\}\right)^{2}>0$ .
This completes the proof that $E$ is small.

That the chain is irreducible, follows easily, for any $\mathbf{x}$ , the set $\{\mathbf{x}\}$ is a compact set and for a measurable set $A$ with $\lambda(A)>0$ we may choose $C$ in the first part of the proof such that $\lambda(C\cap A)>0$ . Now,

K^{2}(\mathbf{x},A)\quad\geq\quad c\lambda(C\cap A)>0

Also aperiodicity follows trivially from the observation that any set with positive $\lambda$ -measure can be accessed in at most 2 steps.

S-4 General TMCMC algorithm with single $\epsilon$ and dependent $\mathbf{z}$

Also, let Let $h_{1}(\mathbf{p})$ , $h_{2}(\mathbf{q})$ be the specified joint distributions of $\mathbf{p}$ and $\mathbf{q}$ induced by the Gaussian distributions of $\mathbf{w}_{1},\mathbf{w}_{2},\mathbf{w}_{3}$ , and let $P(\mathbf{z}|\mathbf{p},\mathbf{q})=\prod_{i=1}^{k}f_{i}(z_{i}|p_{i},q_{i})$ denote the conditional probability of $\mathbf{z}$ , given $\mathbf{p},\mathbf{q}$ , where $f_{i}(\cdot|p_{i},q_{i})$ is the conditional probability of $z_{i}$ given $p_{i}$ and $q_{i}$ . Then the general TMCMC algorithm with singleton $\epsilon$ and dependent $\mathbf{z}$ is given as follows.

Algorithm S-4.1

General TMCMC algorithm based on single $\epsilon$ and dependent $\mathbf{z}$ .

•

Input: Initial value ${\mathbf{x}}^{(0)}$ , and number of iterations $N$ .

1.
For $t=0,\ldots,N-1$
1. (a)
  
  Generate $\mathbf{w}_{1}\sim N_{k}(\mbox{\boldmath$\mu$}_{1},\mbox{\boldmath$\Sigma$}_{1})$ , $\mathbf{w}_{2}\sim N_{k}(\mbox{\boldmath$\mu$}_{2},\mbox{\boldmath$\Sigma$}_{2})$ , and $\mathbf{w}_{3}\sim N_{k}(\mu_{3},\mbox{\boldmath$\Sigma$}_{3})$ .
2. (b)
  
  For $i=1,\ldots,k$ , set $p_{i}=\exp\left(w_{1i}\right)/\sum_{j=1}^{3}\exp\left(w_{ji}\right)$ , $q_{i}=\exp\left(w_{2i}\right)/\sum_{j=1}^{3}\exp\left(w_{ji}\right)$ , and $1-p_{i}-q_{i}=\exp\left(w_{3i}\right)/\sum_{j=1}^{3}\exp\left(w_{ji}\right)$ .
3. (c)
  
  Generate $\epsilon\sim g(\cdot)$ and an index $i\sim\mathcal{M}(1;p_{1},\ldots,p_{3^{k}-1})$ independently.

2.

\mathbf{x}^{\prime}=T_{\mathbf{z}_{i}}({\mathbf{x}}^{(t)},\epsilon)\quad\textrm{ and }\quad\alpha({\mathbf{x}}^{(t)},\epsilon)=\min\left(1,\dfrac{P(\mathbf{z}^{c}_{i}|\mathbf{p},\mathbf{q})}{P(\mathbf{z}_{i}|\mathbf{p},\mathbf{q})}~\dfrac{\pi(\mathbf{x}^{\prime})}{\pi({\mathbf{x}}^{(t)})}~\left|\frac{\partial(T_{\mathbf{z}_{i}}({\mathbf{x}}^{(t)},\epsilon),\epsilon)}{\partial({\mathbf{x}}^{(t)},\epsilon)}\right|\right)

3.

Set

{\mathbf{x}}^{(t+1)}=\left\{\begin{array}[]{ccc}\mathbf{x}^{\prime}&\textsf{ with probability }&\alpha({\mathbf{x}}^{(t)},\epsilon)\\ {\mathbf{x}}^{(t)}&\textsf{ with probability }&1-\alpha({\mathbf{x}}^{(t)},\epsilon)\end{array}\right.

•

End for

S-5 Proof of detailed balance for TMCMC with dependent $\mathbf{z}$

Let $\mathbf{y}=T_{\mathbf{z}}(\mathbf{x},\epsilon)\in T_{\mathbf{z}}(\mathbf{x},\mathcal{Y})$ , then $\mathbf{x}=T_{\mathbf{z}^{c}}(\mathbf{y},\epsilon)$ . The kernel $K$ satisfies,

	$\displaystyle\pi(\mathbf{x})K(\mathbf{x}\to\mathbf{y})$	$\displaystyle=$	$\displaystyle\pi(\mathbf{x})h_{1}(\mathbf{p})h_{2}(\mathbf{q})P(\mathbf{z}\|\mathbf{p},\mathbf{q})~g(\epsilon)\min\left\{1,\dfrac{P(\mathbf{z}^{c}\|\mathbf{p},\mathbf{q})\pi(\mathbf{y})}{P(\mathbf{z}\|\mathbf{p},\mathbf{q})\pi(\mathbf{x})}J_{\mathbf{z}}(\mathbf{x},\epsilon)\right\}$
		$\displaystyle=$	$\displaystyle h_{1}(\mathbf{p})h_{2}(\mathbf{q})g(\boldsymbol{\epsilon})\min\left\{\pi(\mathbf{x})P(\mathbf{z}\|\mathbf{p},\mathbf{q}),\pi(\mathbf{y})P(\mathbf{z}^{c}\|\mathbf{p},\mathbf{q})J_{\mathbf{z}}(\mathbf{x},\epsilon)\right\}$

and

	$\displaystyle\pi(\mathbf{y})K(\mathbf{y}\to\mathbf{x})$	$\displaystyle=$	$\displaystyle\pi(\mathbf{y})h_{1}(\mathbf{p})h_{2}(\mathbf{q})P(\mathbf{z}^{c}\|\mathbf{p},\mathbf{q})g(\epsilon)J_{\mathbf{z}}(\mathbf{x},\epsilon)\min\left\{1,\dfrac{P(\mathbf{z}\|\mathbf{p},\mathbf{q})\pi(\mathbf{x})}{P(\mathbf{z}^{c}\|\mathbf{p},\mathbf{q})\pi(\mathbf{y})}J_{\mathbf{z}^{c}}(\mathbf{y},\epsilon)~\right\}$
		$\displaystyle=$	$\displaystyle h_{1}(\mathbf{p})h_{2}(\mathbf{q})g(\epsilon)\min\left\{\pi(\mathbf{y})P(\mathbf{z}^{c}\|\mathbf{p},\mathbf{q})J_{\mathbf{z}}(\mathbf{x},\epsilon),\pi(\mathbf{x})P(\mathbf{z}\|\mathbf{p},\mathbf{q})\right\}$

S-6 Improved acceptance rates of additive TMCMC with singleton $\epsilon$ compared to joint updating using RWMH

The joint RWMH algorithm generates $\boldsymbol{\epsilon}=(\epsilon_{1},\ldots,\epsilon_{k})^{\prime}$ independently from $N(0,1)$ , and then uses the transformation of the form $x^{\prime}_{i}=x_{i}+a_{i}\epsilon_{i}$ , where $a_{i}>0$ are appropriate scaling constants. For large $k$ , the so-called “curse of dimensionality” can force the acceptance rate to be close to zero. On the other hand, the additive-transformation based TMCMC also updates $(x_{1},\ldots,x_{k})$ simultaneously in a single block, but instead of using $k$ different $\epsilon_{i}$ , it uses a single $\epsilon$ for updating all the $x_{i}$ variables. In other words, for TMCMC based on additive transformation $\boldsymbol{\epsilon}$ is of the form $\boldsymbol{\epsilon}=(\pm\epsilon,\ldots,\pm\epsilon)^{\prime}$ , where $\epsilon\sim N(0,1)\mathbb{I}_{\{\epsilon>0\}}$ . Thus, relative to RWMH, the dimension in the TMCMC case is effectively reduced to 1, avoiding the curse of dimensionality. Thus, it is expected that additive TMCMC will have a much higher acceptance rate than RWMH. In this section we formalize and compare the issues related to acceptance rates of additive TMCMC and RWMH.

S-6.1 Discussion on optimal scaling and optimal acceptance rate of additive TMCMC and RWMH

A reasonable approach to compare the acceptance rates of additive TMCMC and RWMH is to develop the optimal scaling theory for additive TMCMC, obtain the optimal acceptance rate, and then compare the latter with the optimal acceptance rates for RWMH, which are already established in the MCMC literature. Indeed, optimal scaling and optimal acceptance rate of additive TMCMC and comparison with those of RWMH is the subject of Dey and Bhattacharya (2013), where it is shown that additive TMCMC has a much higher optimal acceptance rate compared to RWMH. Before we summarize the results of Dey and Bhattacharya (2013) we first provide a brief overview of optimal scaling and optimal acceptance rate of RWMH.

S-6.1.1 Brief overview of optimal scaling and optimal acceptance rate for RWMH

Roughly, the optimal random walk proposal variance, represented as an inverse function of the dimension $k$ , is the one that maximizes the speed of convergence to the stationary distribution of a relevant diffusion process to which a ‘sped-up’ version of RWMH weakly converges as the dimension $k$ increases to infinity. The optimal acceptance rate corresponds to the optimal proposal variance. Under various assumptions on the form of the target distribution $\pi$ , ranging from the $iid$ assumption (Roberts et al. (1997)), through independent but non-identical set-up (Bedard (2007)), to a more general dependent structure (Mattingly et al. (2011)), the optimal acceptance rate turns out to be 0.234.

S-6.1.2 Optimal scaling and optimal acceptance rate for additive TMCMC

In Dey and Bhattacharya (2013) it has been proved in the case of additive TMCMC, assuming $p_{i}=q_{i}=1/2$ , that the optimal acceptance rate, as $k\rightarrow\infty$ , is 0.439 under the set-ups ( $iid$ , independent but non-identical, and dependent) for which the optimal acceptance rate for RWMH has been studied and established to be 0.234. Thus, the optimal acceptance rate for additive TMCMC is much higher than that of RWMH. The optimal scalings, that is, the optimal values of the scales $a_{1},\ldots,a_{k}$ are also available using the optimal scaling theory. As shown in Dey and Bhattacharya (2013), all these results for additive TMCMC and RWMH remain true even in all the aforementioned set-ups if some of the co-ordinates of $\mathbf{x}$ are updated at random, conditioning on the remaining co-ordinates.

S-6.2 Comparison between the asymptotic forms of the acceptance rates of additive TMCMC and RWMH for strongly log-concave target densities

The results on optimal scaling and optimal acceptance rate discussed in Sections S-6.1.1 and S-6.1.2 are available only for special forms of the target distribution $\pi$ . In this section we obtain the asymptotic forms of the acceptance rates associated with RWMH and additive TMCMC assuming that the target density is strongly log-concave. In particular, under suitable conditions we show that as the dimension increases, the acceptance rate of RWMH converges to zero at a much faster rate than that of additive TMCMC.

Assuming without loss of generality that the marginal variances of the target density $\pi$ are all unity (achieved after suitable scaling perhaps), for RWMH we consider the following proposed value $\mathbf{x}^{\prime}$ given the current value $\mathbf{x}$ : $\mathbf{x}^{\prime}=\mathbf{x}+\boldsymbol{\epsilon}$ , where $\boldsymbol{\epsilon}\sim N_{k}(\mbox{\boldmath$0$},\mathbf{I}_{k})$ . On the other hand, for additive TMCMC, we consider $\mathbf{x}^{\prime}=\mathbf{x}+\epsilon\mbox{\boldmath$\delta$}$ , where $\epsilon\sim N(0,1)\mathbb{I}(\epsilon>0)$ and the components $\delta_{i}$ of $\boldsymbol{\delta}$ are $iid$ taking values $\pm 1$ with probability $1/2$ each.

To proceed we consider the following form of acceptance rate for our asymptotic framework. Letting $R(\mathbf{x}^{\prime}|\mathbf{x})$ denote the acceptance probability of $\mathbf{x}^{\prime}$ given the current value $\mathbf{x}$ , and letting $U\sim Uniform(0,1)$ , the acceptance rate is given by

$\displaystyle AR$	$\displaystyle=\int R(\mathbf{x}^{\prime}\|\mathbf{x})q(\mathbf{x}^{\prime}\|\mathbf{x})\pi(\mathbf{x})d\mathbf{x}d\mathbf{x}^{\prime}$
	$\displaystyle=\int Pr\left(U<R(\mathbf{x}^{\prime}\|\mathbf{x})\right)q(\mathbf{x}^{\prime}\|\mathbf{x})\pi(\mathbf{x})d\mathbf{x}d\mathbf{x}^{\prime}$
	$\displaystyle=\int\left[\int Pr\left(U<R(\mathbf{x}^{\prime}\|\mathbf{x})\right)q(\mathbf{x}^{\prime}\|\mathbf{x})d\mathbf{x}^{\prime}\right]\pi(\mathbf{x})d\mathbf{x}$
	$\displaystyle=\int\left[\int_{0}^{1}Pr\left(R(\mathbf{x}^{\prime}\|\mathbf{x})>u\right)du\right]\pi(\mathbf{x})d\mathbf{x}$	(S-6.1)

In the above formula for acceptance rate note that $Pr\left(R(\mathbf{x}^{\prime}|\mathbf{x})>u\right)\rightarrow 1$ as $u\rightarrow 0$ and $Pr\left(R(\mathbf{x}^{\prime}|\mathbf{x})>u\right)\rightarrow 0$ as $u\rightarrow 1$ . Hence, given any $\eta_{1}>0,\eta_{2}>0$ , we can choose $\psi_{1},\psi_{2}\in(0,1)$ sufficiently small such that $\int_{0}^{\psi_{1}}Pr\left(R(\mathbf{x}^{\prime}|\mathbf{x})>u\right)du<\eta_{1}$ and $\int_{1-\psi_{2}}^{1}Pr\left(R(\mathbf{x}^{\prime}|\mathbf{x})>u\right)du<\eta_{2}$ . Hence, re-writing (S-6.1) as

	$\displaystyle AR$	$\displaystyle=$	$\displaystyle\int\left[\int_{0}^{\psi_{1}}Pr\left(R(\mathbf{x}^{\prime}\|\mathbf{x})>u\right)du\right]\pi(\mathbf{x})d\mathbf{x}+\int\left[\int_{\psi_{1}}^{1-\psi_{2}}Pr\left(R(\mathbf{x}^{\prime}\|\mathbf{x})>u\right)du\right]\pi(\mathbf{x})d\mathbf{x}$
			$\displaystyle\quad\quad\quad\quad+\int\left[\int_{1-\psi_{2}}^{1}Pr\left(R(\mathbf{x}^{\prime}\|\mathbf{x})>u\right)du\right]\pi(\mathbf{x})d\mathbf{x},$

we find that the first and the third term on the right hand side are negligible for any algorithm. So, for the purpose of comparing algorithms with respect to their acceptance rates, we consider only the middle term; in all that follow we denote

AR=\int\left[\int_{\psi_{1}}^{1-\psi_{2}}Pr\left(R(\mathbf{x}^{\prime}|\mathbf{x})>u\right)du\right]\pi(\mathbf{x})d\mathbf{x}.

(S-6.2)

For our purpose, we consider a target density $\pi(\mathbf{x})$ of $k$ variables that is strongly log-concave, that is,

-M_{k}\mathbf{I}_{k}\leq\nabla^{2}\log\pi(\mathbf{x})\leq-m_{k}\mathbf{I}_{k},

(S-6.3)

where we assume that $M_{k}=c_{k}+m_{k}$ , with $m_{k},~c_{k}>0$ for every $k$ . We further assume that $m_{k}\rightarrow\infty$ , and the sequence $\{c_{k}\}$ is such that $c_{k}/m_{k}\rightarrow 0$ as $k\rightarrow\infty$ . Then clearly, $M_{k}\asymp m_{k}$ , meaning $M_{k}/m_{k}\rightarrow 1$ as $k\rightarrow\infty$ . In fact, we assume that $M_{k}/m_{k}$ approaches 1 at a sufficiently fast rate, so that $k\left|\frac{M_{k}}{m_{k}}-1\right|\rightarrow 0$ . For our purpose we assume that $c_{k}=O(k^{s});s\geq 1$ and $m_{k}=O(k^{t});t>s+1\geq 2$ , so that $M_{k}=O(k^{t})$ . It is easy to verify that these choices satisfy the above conditions.

It is important to note that our assumption $m_{k},M_{k}\rightarrow\infty$ need not hold for all strongly log-concave distributions. For instance, when $\pi$ is the $iid$ product of standard normals, that is, when $\mathbf{x}\sim N_{k}\left(\mbox{\boldmath$0$},\mathbf{I}_{k}\right)$ under $\pi$ , $\nabla^{2}\log\pi(\mathbf{x})=\mathbf{I}_{k}$ . In this case, $m_{k}=M_{k}=1$ for every $k\geq 1$ . In general, even if $m_{k}$ and $M_{k}$ remains finite as $k$ grows to infinity, our proofs remain valid provided that $M_{k}\asymp m_{k}$ and $k\left|\frac{M_{k}}{m_{k}}-1\right|\rightarrow 0$ . The case of $\pi$ being an $iid$ product of standard normals clearly satisfies the above conditions.

S-6.2.1 Asymptotic form of the acceptance rate for RWMH

Let $\mathbf{x}^{*}$ denote the mode of the target density $\pi(\cdot)$ . Then for every $r\in(0,1),$

			$\displaystyle Pr~(R(\mathbf{x}^{\prime}\|\mathbf{x})<r)=Pr~(\pi(\mathbf{x}^{\prime})/\pi(\mathbf{x})<r)=Pr~(\log\pi(\mathbf{x}^{\prime})-\log\pi(\mathbf{x})<\log r)$
		$\displaystyle=$	$\displaystyle Pr~(\left[\log\pi(\mathbf{x}^{\prime})-\log\pi(\mathbf{x^{}})\right]-\left[\log\pi(\mathbf{x}^{\prime})-\log\pi(\mathbf{x^{}})\right]<\log r)$
		$\displaystyle=$	$\displaystyle Pr~\left(\left[\nabla\log\pi(\mathbf{x^{}})^{T}(\mathbf{x}^{\prime}-\mathbf{x^{}})+(1/2)(\mathbf{x}^{\prime}-\mathbf{x^{}})^{T}\nabla^{2}\log\pi(\boldsymbol{\xi_{1}(\mathbf{x}^{\prime},\mathbf{x}^{})})(\mathbf{x}^{\prime}-\mathbf{x^{*}})\right]\right.$
			$\displaystyle\quad\quad\left.-\left[\nabla\log\pi(\mathbf{x^{}})^{T}(\mathbf{x}-\mathbf{x^{}})+(1/2)(\mathbf{x}-\mathbf{x^{}})^{T}\nabla^{2}\log\pi(\boldsymbol{\xi_{2}(\mathbf{x},\mathbf{x}^{})})(\mathbf{x}-\mathbf{x^{*}})\right]<\log r\right),$
			$\displaystyle\quad\textrm{~for some }\boldsymbol{\xi}_{1}(\mathbf{x}^{\prime},\mathbf{x}^{}),\boldsymbol{\xi}_{2}(\mathbf{x},\mathbf{x}^{})\ \ \mbox{depending upon}\ \ (\mathbf{x}^{\prime},\mathbf{x}^{})\ \ \mbox{and}(\mathbf{x},\mathbf{x}^{})\ \ \mbox{respectively};$
		$\displaystyle=$	$\displaystyle Pr~\left(\left[(1/2)(\mathbf{x}^{\prime}-\mathbf{x^{}})^{T}\nabla^{2}\log\pi(\boldsymbol{\xi_{1}(\mathbf{x}^{\prime},\mathbf{x}^{})})(\mathbf{x}^{\prime}-\mathbf{x^{*}})\right]\right.$
			$\displaystyle\quad\quad\left.-\left[(1/2)(\mathbf{x}-\mathbf{x^{}})^{T}\nabla^{2}\log\pi(\boldsymbol{\xi_{2}(\mathbf{x}^{\prime},\mathbf{x}^{})})(\mathbf{x}-\mathbf{x^{*}})\right]<\log r\right)$
			$\displaystyle\quad\quad\quad\quad\quad\mbox{since}\ \ \nabla\log\pi(\mathbf{x^{*}})=\mbox{\boldmath$0$}.$

Thus from the assumption in (S-6.3), and noting that $(\mathbf{x}^{\prime}-\mathbf{x^{*}})^{T}(\mathbf{x}^{\prime}-\mathbf{x^{*}})=(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})+2(\mathbf{x}-\mathbf{x^{*}})^{T}\boldsymbol{\epsilon}+\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}$ it follows that

\begin{split}&Pr~\left(\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-m_{k}(\mathbf{x}-\mathbf{x^{*}})^{T}\boldsymbol{\epsilon}-\frac{m_{k}}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}<\log r\right)\\ &\leq Pr~(R(\mathbf{x}^{\prime}|\mathbf{x})<r)\\ &\leq Pr~\left(-\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-M_{k}(\mathbf{x}-\mathbf{x^{*}})^{T}\boldsymbol{\epsilon}-\frac{M_{k}}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}<\log r\right);\end{split}

(S-6.4)

so that

\begin{split}&Pr~\left(\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-m_{k}(\mathbf{x}-\mathbf{x^{*}})^{T}\boldsymbol{\epsilon}-\frac{m_{k}}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}>\log r\right)\\ &\geq Pr~(R(\mathbf{x}^{\prime}|\mathbf{x})>r)\\ &\geq Pr~\left(-\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-M_{k}(\mathbf{x}-\mathbf{x^{*}})^{T}\boldsymbol{\epsilon}-\frac{M_{k}}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}>\log r\right),\end{split}

(S-6.5)

Hence, using (S-6.2) it can be seen that the acceptance rate is bounded above and below as follows

	$\displaystyle\int\left[\int_{\psi_{1}}^{1-\psi_{2}}\left\{\int_{A^{k}_{2,\boldsymbol{\epsilon},u}}\frac{1}{(2\pi)^{k/2}}\exp\left\{-\frac{1}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}\right\}d\boldsymbol{\epsilon}\right\}du\right]\pi(\mathbf{x})d\mathbf{x}$
	$\displaystyle\geq AR^{(RWMH)}$		(S-6.6)
	$\displaystyle\geq\int\left[\int_{\psi_{1}}^{1-\psi_{2}}\left\{\int_{A^{k}_{1,\boldsymbol{\epsilon},u}}\frac{1}{(2\pi)^{k/2}}\exp\left\{-\frac{1}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}\right\}d\boldsymbol{\epsilon}\right\}du\right]\pi(\mathbf{x})d\mathbf{x},$

where

A^{k}_{1,\boldsymbol{\epsilon},u}=\left\{\mathbf{x}:-\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-M_{k}(\mathbf{x}-\mathbf{x^{*}})^{T}\boldsymbol{\epsilon}-\frac{M_{k}}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}>\log u\right\}

and

A^{k}_{2,\boldsymbol{\epsilon},u}=\left\{\mathbf{x}:\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-m_{k}(\mathbf{x}-\mathbf{x^{*}})^{T}\boldsymbol{\epsilon}-\frac{m_{k}}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}>\log u\right\}.

Now, note that for some $\mbox{\boldmath$\xi$}(\mathbf{x},\mathbf{x}^{*})$ depending upon $\mathbf{x}$ and $\mathbf{x}^{*}$ ,

	$\displaystyle\pi(\mathbf{x})$	$\displaystyle=$	$\displaystyle\exp\left\{\log\pi(\mathbf{x})\right\}d\mathbf{x}$
		$\displaystyle=$	$\displaystyle\exp\left\{\log\pi(\mathbf{x^{}})+\frac{1}{2}(\mathbf{x}-\mathbf{x^{}})^{T}\nabla^{2}\log\pi(\boldsymbol{\xi}(\mathbf{x},\mathbf{x}^{}))(\mathbf{x}-\mathbf{x^{}})\right\},$

so that the inequalities related to strong convexity, given by (S-6.3) yield

\frac{(2\pi)^{k/2}}{m^{k}_{k}}\pi(\mathbf{x}^{*})\frac{m^{k}_{k}}{(2\pi)^{k/2}}\exp\left\{-\frac{m_{k}}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})\right\}\geq\pi(\mathbf{x})

\geq\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\frac{M^{k}_{k}}{(2\pi)^{k/2}}\exp\left\{-\frac{M_{k}}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})\right\}

(S-6.8)

Using the lower bound of $\pi(\mathbf{x})$ given by (S-6.8) and Fubini’s theorem, the lower bound of the acceptance rate given by (S-6.6) can be further bounded below as

$\displaystyle AR^{(RWMH)}$	$\displaystyle\geq\int\int_{\psi_{1}}^{1-\psi_{2}}\int_{A^{k}_{1,\boldsymbol{\epsilon},u}}\frac{1}{(2\pi)^{k/2}}\exp\left\{-\frac{1}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}\right\}\pi(\mathbf{x})~d\mathbf{x}~du~d\boldsymbol{\epsilon}$
	$\displaystyle\geq\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\int\int_{\psi_{1}}^{1-\psi_{2}}\int_{A^{k}_{1,\boldsymbol{\epsilon},u}}\frac{1}{(2\pi)^{k/2}}\exp\left\{-\frac{1}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}\right\}$
	$\displaystyle\quad\quad\quad\quad\times\frac{M^{k}_{k}}{(2\pi)^{k/2}}\exp\left\{-\frac{M_{k}}{2}(\mathbf{x}-\mathbf{x^{}})^{T}(\mathbf{x}-\mathbf{x^{}})\right\}~d\mathbf{x}~du~d\boldsymbol{\epsilon}$
	$\displaystyle\geq\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\inf_{u\in(\psi_{1},1-\psi_{2})}Pr~(A^{k}_{1,\boldsymbol{\epsilon},u}),$	(S-6.9)

where $Pr~(A^{k}_{1,\boldsymbol{\epsilon},u})$ must be calculated with respect to $\boldsymbol{\epsilon}\sim N_{k}(\mbox{\boldmath$0$},\mathbf{I}_{k})$ , and independently, $\mathbf{x}-\mathbf{x}^{*}\sim N_{k}(\mbox{\boldmath$0$},M^{-1}_{k}\mathbf{I}_{k})$ .

Similarly, using the upper bound of $\pi(\mathbf{x})$ given by (S-6.8) the upper bound of the acceptance rate given by (S-6.6) can be further bounded above as

$\displaystyle AR^{(RWMH)}$	$\displaystyle\leq\int\int_{\psi_{1}}^{1-\psi_{2}}\int_{A^{k}_{2,\boldsymbol{\epsilon},u}}\frac{1}{(2\pi)^{k/2}}\exp\left\{-\frac{1}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}\right\}\pi(\mathbf{x})~d\mathbf{x}~du~d\boldsymbol{\epsilon}$
	$\displaystyle=\frac{(2\pi)^{k/2}}{m^{k}_{k}}\pi(\mathbf{x}^{*})\int_{\psi_{1}}^{1-\psi_{2}}Pr~(A^{k}_{2,\boldsymbol{\epsilon},u})du$
	$\displaystyle\leq\frac{(2\pi)^{k/2}}{m^{k}_{k}}\pi(\mathbf{x}^{*})\sup_{u\in(\psi_{1},1-\psi_{2})}Pr~(A^{k}_{2,\boldsymbol{\epsilon},u})$	(S-6.10)

The probability $Pr~(A^{k}_{2,\boldsymbol{\epsilon},u_{k}})$ must be calculated with respect to $\boldsymbol{\epsilon}\sim N_{k}(\mbox{\boldmath$0$},\mathbf{I}_{k})$ , and independently, $\mathbf{x}-\mathbf{x}^{*}\sim N_{k}(\mbox{\boldmath$0$},m^{-1}_{k}\mathbf{I}_{k})$ . Thus, we have

\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\inf_{u\in(\psi_{1},1-\psi_{2})}Pr~(A^{k}_{1,\boldsymbol{\epsilon},u})\leq AR^{(RWMH)}\leq\frac{(2\pi)^{k/2}}{m^{k}_{k}}\pi(\mathbf{x}^{*})\sup_{u\in(\psi_{1},1-\psi_{2})}Pr~(A^{k}_{2,\boldsymbol{\epsilon},u}).

(S-6.11)

We first focus on the lower bound in (S-6.11). As $k\rightarrow\infty$ ,
$-\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-M_{k}(\mathbf{x}-\mathbf{x^{*}})^{T}\boldsymbol{\epsilon}-\frac{M_{k}}{2}\boldsymbol{\epsilon}^{T}\boldsymbol{\epsilon}$

\displaystyle\sim

\displaystyle AN\left(-\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)+M_{k}\right],\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)^{2}+2M_{k}+M^{2}_{k}\right]\right),

(S-6.12)

where $AN(\mu,\sigma^{2})$ denotes asymptotic normal with mean $\mu$ and variance $\sigma^{2}$ . From (S-6.12) it follows that

	$\displaystyle\inf_{u\in(\psi_{1},1-\psi_{2})}Pr~(A^{k}_{1,\boldsymbol{\epsilon},u})$	$\displaystyle\asymp$	$\displaystyle 1-\sup_{u\in(\psi_{1},1-\psi_{2})}\Phi\left(\frac{\log u+\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)+M_{k}\right]}{\sqrt{\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)^{2}+2M_{k}+M^{2}_{k}\right]}}\right)$		(S-6.13)
		$\displaystyle=$	$\displaystyle 1-\Phi\left(\frac{\log(1-\psi_{2})+\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)+M_{k}\right]}{\sqrt{\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)^{2}+2M_{k}+M^{2}_{k}\right]}}\right).$		(S-6.13)

Combining (S-6.9) and (S-6.13) we obtain

AR^{(RWMH)}\geq\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\left\{1-\Phi\left(\frac{\log(1-\psi_{2})+\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)+M_{k}\right]}{\sqrt{\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)^{2}+2M_{k}+M^{2}_{k}\right]}}\right)\right\}.

(S-6.14)

Now focusing our attention on the upper bound of $AR^{(RWMH)}$ we similarly obtain

\displaystyle AR^{(RWMH)}

\displaystyle\leq

\displaystyle\frac{(2\pi)^{k/2}}{m^{k}_{k}}\pi(\mathbf{x}^{*})\left\{1-\Phi\left(\frac{\log\psi_{1}-\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{m_{k}}\right)-m_{k}\right]}{\sqrt{\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{m_{k}}\right)^{2}+2m_{k}+m^{2}_{k}\right]}}\right)\right\}.

In other words,

	$\displaystyle\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})$	$\displaystyle\left\{1-\Phi\left(\frac{\log(1-\psi_{2})+\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)+M_{k}\right]}{\sqrt{\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)^{2}+2M_{k}+M^{2}_{k}\right]}}\right)\right\}\leq AR^{(RWMH)}$
		$\displaystyle\leq\frac{(2\pi)^{k/2}}{m^{k}_{k}}\pi(\mathbf{x}^{*})\left\{1-\Phi\left(\frac{\log\psi_{1}-\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{m_{k}}\right)-m_{k}\right]}{\sqrt{\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{m_{k}}\right)^{2}+2m_{k}+m^{2}_{k}\right]}}\right)\right\}.$

Since $m_{k}\asymp M_{k}$ , it is easy to see that

		$\displaystyle\frac{\log(1-\psi_{2})+\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)+M_{k}\right]}{\sqrt{\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)^{2}+2M_{k}+M^{2}_{k}\right]}}\asymp\sqrt{\frac{k}{2}},\quad\mbox{and}$
		$\displaystyle\frac{\log\psi_{1}-\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{m_{k}}\right)-m_{k}\right]}{\sqrt{\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{m_{k}}\right)^{2}+2m_{k}+m^{2}_{k}\right]}}\asymp\sqrt{\frac{k}{2}}.$

Hence, it follows that

AR^{(RWMH)}\asymp\frac{(2\pi)^{k/2}}{M^{k}_{k}}\left\{1-\Phi\left(\sqrt{\frac{k}{2}}\right)\right\}.

(S-6.17)

S-6.2.2 Asymptotic bounds of the acceptance rate for additive TMCMC

Next let us obtain lower and upper bounds of $AR^{(TMCMC)}$ associated with TMCMC with additive transformation. Recall that in this case, $\mathbf{x}^{\prime}=\mathbf{x}+\epsilon\boldsymbol{\delta}$ where $\epsilon\sim N(0,1)\mathbb{I}(\epsilon>0)$ and the components $\delta_{i}$ of $\boldsymbol{\delta}$ are $iid$ taking values $\pm 1$ with probability $1/2$ each. In this set up (S-6.5) becomes

\begin{split}&Pr~\left(\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-m_{k}\epsilon(\mathbf{x}-\mathbf{x^{*}})^{T}\mbox{\boldmath$\delta$}-\frac{m_{k}}{2}k\epsilon^{2}>\log r\right)\\ &\geq Pr~(R(\mathbf{x}^{\prime}|\mathbf{x})>r)\\ &\geq Pr~\left(-\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-M_{k}\epsilon(\mathbf{x}-\mathbf{x^{*}})^{T}\mbox{\boldmath$\delta$}-\frac{M_{k}}{2}k\epsilon^{2}>\log r\right),\end{split}

(S-6.18)

Now notice that, under the lower bound of $\pi(\mathbf{x})$ provided in (S-6.8), as $k\rightarrow\infty$ ,

\frac{M_{k}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})}{k}\stackrel{{\scriptstyle\alpha.s.}}{{\longrightarrow}}1,

and

\frac{\sqrt{M_{k}}(\mathbf{x}-\mathbf{x^{*}})^{T}\mbox{\boldmath$\delta$}}{k}\stackrel{{\scriptstyle\alpha.s.}}{{\longrightarrow}}0.

Similarly, under the upper bound of $\pi(\mathbf{x})$ in (S-6.8), the above hold with $M_{k}$ replaced with $m_{k}$ . From these it follow that the asymptotic forms of the lower and the upper bounds of (S-6.18) are given by

Pr~\left(-\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-M_{k}\epsilon(\mathbf{x}-\mathbf{x^{*}})^{T}\mbox{\boldmath$\delta$}-\frac{M_{k}}{2}k\epsilon^{2}>\log r\right)

\asymp\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\left\{2\Phi\left(\sqrt{-\frac{2}{kM_{k}}\log r-\left(\frac{M_{k}-m_{k}}{M^{2}_{k}}\right)}\right)-1\right\}

and

Pr~\left(\frac{(M_{k}-m_{k})}{2}(\mathbf{x}-\mathbf{x^{*}})^{T}(\mathbf{x}-\mathbf{x^{*}})-m_{k}\epsilon(\mathbf{x}-\mathbf{x^{*}})^{T}\mbox{\boldmath$\delta$}-\frac{m_{k}}{2}k\epsilon^{2}>\log r\right)

\asymp\frac{(2\pi)^{k/2}}{m^{k}_{k}}\pi(\mathbf{x}^{*})\left\{2\Phi\left(\sqrt{-\frac{2}{km_{k}}\log r+\left(\frac{M_{k}-m_{k}}{m^{2}_{k}}\right)}\right)-1\right\}.

Using the above results, it follows as in the case of $AR^{(RWMH)}$ that

	$\displaystyle\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\left\{2\inf_{u\in(\psi_{1},1-\psi_{2})}\Phi\left(\sqrt{-\frac{2}{kM_{k}}\log u-\left(\frac{M_{k}-m_{k}}{M^{2}_{k}}\right)}\right)-1\right\}$
	$\displaystyle\leq AR^{(TMCMC)}\leq\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\left\{2\sup_{u\in(\psi_{1},1-\psi_{2})}\Phi\left(\sqrt{-\frac{2}{km_{k}}\log u-\left(\frac{m_{k}-M_{k}}{m^{2}_{k}}\right)}\right)-1\right\}.$

Substituting the infimum and supremum over $(\psi_{1},1-\psi_{2})$ we obtain

	$\displaystyle\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\left\{2\Phi\left(\sqrt{-\frac{2}{kM_{k}}\log(1-\psi_{2})-\left(\frac{M_{k}-m_{k}}{M^{2}_{k}}\right)}\right)-1\right\}$
	$\displaystyle\leq AR^{(TMCMC)}\leq\frac{(2\pi)^{k/2}}{m^{k}_{k}}\pi(\mathbf{x}^{*})\left\{2\Phi\left(\sqrt{-\frac{2}{km_{k}}\log\psi_{1}-\left(\frac{m_{k}-M_{k}}{m^{2}_{k}}\right)}\right)-1\right\}.$

Since $k\left|\frac{M_{k}}{m_{k}}-1\right|\rightarrow 0$ and $m_{k}\asymp M_{k}$ , it follows that

	$\displaystyle-\frac{2}{kM_{k}}\log(1-\psi_{2})-\left(\frac{M_{k}-m_{k}}{M^{2}_{k}}\right)\asymp-\frac{2}{kM_{k}}\log(1-\psi_{2})\quad\mbox{and}$
	$\displaystyle-\frac{2}{km_{k}}\log\psi_{1}-\left(\frac{m_{k}-M_{k}}{m^{2}_{k}}\right)\asymp-\frac{2}{km_{k}}\log\psi_{1}\asymp-\frac{2}{kM_{k}}\log\psi_{1}.$

Hence,

	$\displaystyle\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\left\{2\Phi\left(\sqrt{-\frac{2}{kM_{k}}\log(1-\psi_{2})}\right)-1\right\}$	$\displaystyle\leq$	$\displaystyle AR^{(TMCMC)}$
		$\displaystyle\leq$	$\displaystyle\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})\left\{2\Phi\left(\sqrt{-\frac{2}{kM_{k}}\log\psi_{1}}\right)-1\right\}.$

For comparing (LABEL:eq:AR_TMCMC_bounds) with (S-6.17) where $M_{k}=O\left(k^{t}\right);t>2$ , it can be easily verified using L’Hospital’s rule that for any $\zeta_{1}>0$ , $\zeta_{2}>0$ ,

\frac{2\Phi\left(\frac{\zeta_{1}}{\sqrt{kM_{k}}}\right)-1}{1-\Phi\left(\zeta_{2}\sqrt{k}\right)}\rightarrow\infty.

(S-6.20)

The above result will continue to hold if instead of $M_{k}=O\left(k^{t}\right);t>2$ , $M_{k}\rightarrow a$ , where $a>0$ is some constant. Hence, $AR^{(TMCMC)}$ converges to zero at a much slower rate compared to $AR^{(RWMH)}$ .

S-7 Comparison of TMCMC with HMC

Motivated by Hamiltonian dynamics, Duane et al. (1987) introduced HMC, an MCMC algorithm with deterministic proposals based on approximations of the Hamiltonian equations. We will show that this algorithm is a special case of TMCMC, but first we provide a brief overview of HMC. More details can be found in Liu (2001), Cheung and Beck (2009) and the references therein.

S-7.1 Overview of HMC

If $\pi(\mathbf{x})$ is the target distribution, a fictitious dynamical system may be considered, where $\mathbf{x}(t)\in\mathbb{R}^{d}$ can be thought of as the $d$ -dimensional position vector of a body of particles at time $t$ . If $\mathbf{v}(t)=\dot{\mathbf{x}}(t)=\frac{d{\mathbf{x}}}{dt}$ is the speed vector of the particles, $\dot{\mathbf{v}}(t)=\frac{d{\mathbf{v}}}{dt}$ is its acceleration vector, and $\vec{F}$ is the force exerted on the particle; then, by Newton’s law of motion $\vec{F}=\mathbf{m}\dot{\mathbf{v}}(t)=(m_{1}\dot{v_{1}},\ldots,m_{d}\dot{v_{d}})(t)$ , where $\mathbf{m}\in\mathbb{R}^{d}$ is a mass vector. The momentum vector, $\mathbf{p}=\mathbf{m}\mathbf{v}$ , often used in classical mechanics, can be thought of as a vector of auxiliary variables brought in to facilitate simulation from $\pi(\mathbf{x})$ . The kinetic energy of the system is defined as $W(\mathbf{p})=\mathbf{p}^{\prime}\mathbf{M}^{-1}\mathbf{p}$ , $\mathbf{M}$ being the mass matrix. Usually, $\mathbf{M}$ is taken as $\mathbf{M}=diag\{m_{1},\ldots,m_{d}\}$ .

The target density $\pi(\mathbf{x})$ is linked to the dynamical system via the potential energy field of the system, defined as $U(\mathbf{x})=-\log\pi(\mathbf{x})$ . The total energy (Hamiltonian function), is given by $H(\mathbf{x},\mathbf{p})=U(\mathbf{x})+W(\mathbf{p})$ . A joint distribution over the phase-space $(\mathbf{x},\mathbf{p})$ is then considered, given by

f(\mathbf{x},\mathbf{p})\propto\exp\left\{-H(\mathbf{x},\mathbf{p})\right\}=\pi(\mathbf{x})\exp\left(-\mathbf{p}^{\prime}\mathbf{M}^{-1}\mathbf{p}/2\right)

(S-7.1)

Since the marginal density of $f(\mathbf{x},\mathbf{p})$ is $\pi(\mathbf{x})$ , it now remains to provide a joint proposal mechanism for simulating $(\mathbf{x},\mathbf{p})$ jointly; ignoring $\mathbf{p}$ yields $\mathbf{x}$ marginally from $\pi(\cdot)$ .

For the joint proposal mechanism, HMC makes use of Newton’s law of motion, derived from the law of conservation of energy, and often written in the form of Hamiltonian equations, given by

	$\displaystyle\dot{\mathbf{x}}(t)$	$\displaystyle=$	$\displaystyle\frac{\partial H(\mathbf{x},\mathbf{p})}{\partial\mathbf{p}}=\mathbf{M}^{-1}\mathbf{p},$
	$\displaystyle\dot{\mathbf{p}}(t)$	$\displaystyle=$	$\displaystyle-\frac{\partial H(\mathbf{x},\mathbf{p})}{\partial\mathbf{x}}=-\nabla U(\mathbf{x}),$

where $\nabla U(\mathbf{x})=\frac{\partial U(\mathbf{x})}{\partial\mathbf{x}}$ . The Hamiltonian equations can be approximated by the commonly used leap-frog algorithm (Hockney (1970)), given by,

	$\displaystyle\mathbf{x}(t+\delta t)$	$\displaystyle=\mathbf{x}(t)+\delta t\mathbf{M}^{-1}\left\{\mathbf{p}(t)-\frac{\delta t}{2}\nabla U\left(\mathbf{x}(t)\right)\right\}$		(S-7.2)
	$\displaystyle\mathbf{p}(t+\delta t)$	$\displaystyle=\mathbf{p}(t)-\frac{\delta t}{2}\left\{\nabla U\left(\mathbf{x}(t)\right)+\nabla U\left(\mathbf{x}(t+\delta t)\right)\right\}$		(S-7.3)

Given choices of $\mathbf{M}$ , $\delta t$ , and $L$ , the HMC is then given by the following algorthm:

Algorithm S-7.1

HMC

•

Initialise $\mathbf{x}$ and draw $\mathbf{p}\sim N(\mbox{\boldmath$0$},\mathbf{M})$ .

•

Assuming the current state to be $(\mathbf{x},\mathbf{p})$ , do the following:

1.

Generate $\mathbf{p}^{\prime}\sim N\left(\mbox{\boldmath$0$},\mathbf{M}\right)$ ;
2.

Letting $(\mathbf{x}(0),\mathbf{p}(0))=(\mathbf{x},\mathbf{p}^{\prime})$ , run the leap-frog algorithm for $L$ time steps, to yield $(\mathbf{x}^{\prime\prime},\mathbf{p}^{\prime\prime})=\left(\mathbf{x}(t+L\delta t),\mathbf{p}(t+L\delta t)\right)$ ;

3.

Accept $(\mathbf{x}^{\prime\prime},\mathbf{p}^{\prime\prime})$ with probability

\min\left\{1,\exp\left\{-H(\mathbf{x}^{\prime\prime},\mathbf{p}^{\prime\prime})+H(\mathbf{x},\mathbf{p}^{\prime})\right\}\right\},

(S-7.4)

and accept $(\mathbf{x},\mathbf{p}^{\prime})$ with the remaining probability.

In the above algorithm, it is not required to store simulations of $\mathbf{p}$ . Next we show that HMC is a special case of TMCMC.

S-7.2 HMC is a special case of TMCMC

To see that HMC is a special case of TMCMC, note that the leap-frog step of the HMC algorithm (Algorithm S-7.1) is actually a deterministic transformation of the form $g^{L}:(\mathbf{x}(0),\mathbf{p}(0))\to(\mathbf{x}(L),\mathbf{p}(L))$ (see Liu (2001)). This transformation satisfies the following: if $(\mathbf{x}^{\prime},\mathbf{p}^{\prime})=g^{L}(\mathbf{x},\mathbf{p})$ , then $(\mathbf{x},-\mathbf{p})=g^{L}(\mathbf{x}^{\prime},-\mathbf{p}^{\prime})$ .

The Jacobian of this transformation is 1 because of the volume preservation property, which says that if $V(0)$ is a subset of the phase space, and if $V(t)=\left\{(\mathbf{x}(t),\mathbf{p}(t)):(\mathbf{x}(0),\mathbf{p}(0))\in V(0)\right\}$ , then the volume $|V(t)|={\int\int}_{V(t)}d\mathbf{x}d\mathbf{p}={\int\int}_{V(0)}d\mathbf{x}d\mathbf{p}=|V(0)|$ . As a result, the Jacobian does not feature in the HMC acceptance probability (S-7.4).

For any dimension, there is only one move type defined for HMC, which is the forward transformation $g^{L}$ . Hence, this move type has probability one of selection, and all other move types which we defined in general terms in connection with TMCMC, have zero probability of selection. As a result, the corresponding TMCMC acceptance ratio needs slight modification—it must be made free of the move-type probabilities, which is exactly the case in (S-7.4).

The momentum vector $\mathbf{p}$ can be likened to $\boldsymbol{\epsilon}$ of TMCMC, but note that $\mathbf{p}$ must always be of the same dimensionality as $\mathbf{x}$ ; this is of course, permitted by TMCMC as a special case.

S-7.3 Comparison of acceptance rate for $L=1$ with RWMH and TMCMC

For $L=1$ , the proposal corresponding to HMC is given by (see Cheung and Beck (2009))

q(\mathbf{x}^{\prime}\mid\mathbf{x}(t))=N\left(\mathbf{x}^{\prime}:\mbox{\boldmath$\mu$}(t),\mbox{\boldmath$\Sigma$}(t)\right),

(S-7.5)

where (S-7.5) is a normal distribution with mean and variance given, respectively, by the following:

	$\displaystyle\mbox{\boldmath$\mu$}(t)$	$\displaystyle=\mathbf{x}(t)+\frac{1}{2}\mathbf{M}^{-1}\delta t\nabla\log\left(\pi(\mathbf{x}(t))\right)$		(S-7.6)
	$\displaystyle\mbox{\boldmath$\Sigma$}(t)$	$\displaystyle=\delta t\mathbf{M}^{-1}$		(S-7.7)

Assuming diagonal $\mathbf{M}$ with $m_{i}$ being the $i$ -th diagonal element, the proposal can be re-written in the following more convenient manner: for $i=1,\ldots,k$ ,

x^{\prime}_{i}=x_{i}(t)+\epsilon_{i},

(S-7.8)

where $s_{i}(t)$ denotes the $i$ -th component of $\nabla\log\left(\pi(\mathbf{x}(t))\right)$ , and $\epsilon_{i}\sim N\left(\frac{1}{2}\frac{\delta ts_{i}(t)}{m_{i}},\frac{\delta t}{m_{i}}\right)$ . Assuming, as is usual, that $m_{i}=1$ for each $i$ , it follows that

\frac{\parallel\mathbf{x}^{\prime}-\mathbf{x}\parallel^{2}}{\delta t^{2}}=\sum_{i=1}^{k}\left(\frac{\epsilon_{i}}{\delta t}\right)^{2}=\sum_{i=1}^{k}\epsilon^{\prime 2}_{i}\sim\chi^{2}_{k}(\lambda),

(S-7.9)

where $\chi^{2}_{k}(\lambda)$ is a non-central $\chi^{2}$ distribution with $k$ degrees of freedom and non-centrality parameter $\lambda=\frac{\delta t^{2}}{4}\sum_{i=1}^{k}s^{2}_{i}(t)$ . Since, as either $k\rightarrow\infty$ or $\lambda\rightarrow\infty$ ,

\frac{\chi^{2}_{k}(\lambda)-(k+\lambda)}{\sqrt{2(k+2\lambda)}}\stackrel{{\scriptstyle\mathcal{L}}}{{\rightarrow}}N(0,1),

(S-7.10)

assuming the same strong log-concavity conditions on the target density $\pi$ as provided in Section S-6.2 it follows as in (S-6.2.1) that,

	$\displaystyle\frac{(2\pi)^{k/2}}{M^{k}_{k}}\pi(\mathbf{x}^{*})$	$\displaystyle\left\{1-\Phi\left(\frac{\log(1-\psi_{2})+\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)+M_{k}\delta t^{2}\left(1+\frac{\lambda}{k}\right)\right]}{\sqrt{\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{M_{k}}\right)^{2}+2M_{k}\delta t^{2}\left(1+\frac{\lambda}{k}\right)+M^{2}_{k}\delta t^{4}\left(1+\frac{2\lambda}{k}\right)\right]}}\right)\right\}\leq AR^{(HMC)}$
		$\displaystyle\leq\frac{(2\pi)^{k/2}}{m^{k}_{k}}\pi(\mathbf{x}^{*})\left\{1-\Phi\left(\frac{\log\psi_{1}-\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{m_{k}}\right)-m_{k}\delta t^{2}\left(1+\frac{\lambda}{k}\right)\right]}{\sqrt{\frac{k}{2}\left[\left(\frac{M_{k}-m_{k}}{m_{k}}\right)^{2}+2m_{k}\delta t^{2}\left(1+\frac{\lambda}{k}\right)+m^{2}_{k}\delta t^{4}\left(1+\frac{2\lambda}{k}\right)\right]}}\right)\right\}$

If $\lambda/k\rightarrow 0$ as $k\rightarrow\infty$ , it follows as in Section S-6.2.1 that

AR^{(HMC)}\asymp\frac{(2\pi)^{k/2}}{M^{k}_{k}}\left\{1-\Phi\left(\sqrt{\frac{k}{2}}\right)\right\},

(S-7.12)

which is of the same asymptotic form as (S-6.17), corresponding to the RWMH acceptance rate. On the other hand, if $\lambda/k\rightarrow\infty$ as $k\rightarrow\infty$ , then it follows that

AR^{(HMC)}\asymp\frac{(2\pi)^{k/2}}{M^{k}_{k}}\left\{1-\Phi\left(\frac{\sqrt{\frac{k}{2}\left(1+\frac{\lambda}{k}\right)}}{\sqrt{2}\sqrt{\frac{1}{M_{k}\delta t^{2}}+1}}\right)\right\},

(S-7.13)

which clearly tends to zero at a much faster rate compared to (S-7.12).

To summarize, if $\lambda/k\rightarrow 0$ as $k\rightarrow\infty$ , then both HMC and RWMH have the same asymptotic acceptance rate, tending to zero much faster than that of additive TMCMC. On the other hand, if $\lambda/k\rightarrow\infty$ as $k\rightarrow\infty$ , the acceptance rate of HMC tends to zero much faster than that of RWMH, while that of additive TMCMC maintains its slowest convergence rate to zero. Also observe that the above conclusions will continue to hold if $m_{k}$ and $M_{k}$ tend to finite positive constants satisfying $M_{k}\asymp m_{k}$ and $k\left|\frac{m_{k}}{M_{k}}-1\right|\rightarrow 0$ as $k\rightarrow\infty$ .

S-8 Generalized Gibbs/Metropolis approaches and comparisons with TMCMC

It is important to make it clear at the outset of this discussion that the goals of TMCMC and generalized Gibbs/Metropolis methods are different, even though both use moves based on transformations. While the strength of the latter lies in improving mixing of the standard Gibbs/MH algorithms by adding transformation-based steps to the underlying collection of usual Gibbs/MH steps, TMCMC is an altogether general method of simulating from the target distribution which does not require any underlying step of Gibbs or MH.

The generalized Gibbs/MH methods work in the following manner. Suppose that an underlying Gibbs or MH algorithm for exploring a target distribution has poor mixing properties. Then in order to improve mixing, one may consider some suitable transformation of the random variables being updated such that mixing is improved under the transformation. Such a transformation needs to chosen carefully since it is important to ensure that invariance of the Markov chain is preserved under the transformation. It is convenient to begin with an overview of the generalized Gibbs method with a sequential updating scheme and then proceed to the discussion on the issues and the importance of the block updating idea in the context of improving mixing of standard Gibbs/MH methods.

Liu and Sabatti (2000) (see also Liu and Yu (1999)) propose simulation of a transformation from some appropriate probability distribution, and then applying the transformation to the co-ordinate to be updated. For example, in a $d$ -dimensional target distribution, for updating $\mathbf{x}=(x_{1},x_{2},\ldots,x_{d})$ to $\mathbf{x}^{\prime}=(x^{\prime}_{1},x_{2},\ldots,x_{d})$ , using an additive transformation, one can select $\epsilon$ from some appropriate distribution and set $x^{\prime}_{1}=x_{1}+\epsilon$ . Similarly, if a scale transformation is desired, then one can set $x^{\prime}_{1}=\gamma x_{1}$ , where $\gamma$ must be sampled from some suitable distribution. The suitable distributions of $\epsilon$ and $\gamma$ are chosen such that the target distribution is invariant with respect to the move $\mathbf{x}^{\prime}$ , the forms of which are provided in Liu and Sabatti (2000). For instance, if $\pi(\cdot)$ denotes the target distribution, then for the additive transformation, $\epsilon$ may be sampled from $\pi(x_{1}+\epsilon,x_{2},\ldots,x_{d})$ , and for the multiplicative transformation, one may sample $\gamma$ from $|\gamma|\pi(\gamma x_{1},x_{2},\ldots,x_{d})$ . Since direct sampling from such distributions may be impossible, Liu and Sabatti (2000) suggest a Metropolis-type move with respect to a transformation-invariant transition kernel.

Thus, in the generalized Gibbs method, sequentially all the variables must be updated, unlike TMCMC, where all the variables can be updated simultaneously in a single block. Here we note that for irreducibility issues the generalized Gibbs approach is not suitable for updating the variables blockwise using some transformation that acts on all the variables in a given block. To consider a simple example, with say, $d=2$ and a single block consisting of both the variables, if one considers the additive transformation, then starting with $\mathbf{x}=(x_{1},x_{2})$ , where $x_{1}<x_{2}$ , one can not ever reach $\mathbf{x}^{\prime}=(x^{\prime}_{1},x^{\prime}_{2})$ , where $x^{\prime}_{1}>0,x^{\prime}_{2}<0$ . This is because $x^{\prime}_{1}=x_{1}+z$ and $x^{\prime}_{2}=x_{2}+z$ , for some $z$ , and $x^{\prime}_{1}>0,x^{\prime}_{2}<0$ implies $z>-x_{1}$ and $z<-x_{2}$ , which is a contradiction. The scale transformation implies the move $\mathbf{x}=(x_{1},\ldots,x_{d})\rightarrow(\gamma x_{1},\ldots,\gamma x_{d})=\mathbf{x}^{\prime}$ . If one initializes the Markov chain with all components positive, for instance, then in every iteration, all the variables will have the same sign. The spaces where some variables are positive and some negative will never be visited, even if those spaces have positive (in fact, high) probabilities under the target distribution. This shows that the Markov chain is not irreducible. In fact, with the aforementioned approach, no transformation, whatever distribution they are generated from, can guarantee irreducibility in general if blockwise updates using the transformation strategy of generalized Gibbs is used.

Although blockwise transformations are proposed in Liu and Sabatti (2000) (see also Kou et al. (2005) who propose a MH-based rule for blockwise transformation), they are meant for a different purpose than that discussed above. The strength of such blockwise transformations lies in improving the mixing behaviour of standard Gibbs or MH algorithms. Suppose that an underlying Gibbs or MH algorithm for exploring a target distribution has poor mixing properties. Then in order to improve mixing, one may consider some suitable transformation of the set of random variables being updated such that mixing is improved under the transformation. This additional step involving transformation of the block of random variables can be obtained by selecting a transformation from the appropriate probability distribution provided in Liu and Sabatti (2000). This “appropriate” probability distribution guarantees that stationarity of the transformed block of random variables is preserved. Examples reported in Liu and Sabatti (2000), Müller and Czado (2005), Kou et al. (2005), etc. demonstrate that this transformation also improves the mixing behaviour of the chain, as desired.

Thus, to improve mixing using the methods of Liu and Sabatti (2000) or Kou et al. (2005) one needs to run the usual Gibbs/MH steps, with an additional step involving transformations as discussed above. This additional step induces more computational burden compared to the standard Gibbs/MH steps, but improved mixing may compensate for the extra computational labour. In very high dimensions, of course, this need not be a convenient approach since computational complexity usually makes standard Gibbs/MH approaches infeasible. Since the additional transformation-based step works on the samples generated by standard Gibbs/MH, impracticality of the latter implies that the extra transformation-based step of Liu and Sabatti (2000) for improving mixing is of little value in such cases.

It is important to point out that the generalized Gibbs/MH methods can be usefully employed by even TMCMC to further improve its mixing properties. In other words, a step of generalized Gibbs/MH can be added to the computational fast TMCMC. This additional step can significantly improve the mixing properties of TMCMC. That TMCMC is much faster computationally than standard Gibbs/MH methods imply that even in very high-dimensional situations the generalized Gibbs/MH step can ve very much successful while working in conjunction with TMCMC.

Flight no.	Failure	Temp	Flight no.	Failure	Temp
14	1	53	2	1	70
9	1	57	11	1	70
23	1	58	6	0	72
10	1	63	7	0	73
1	0	66	16	0	75
5	0	67	21	1	75
13	0	67	19	0	76
15	0	67	22	0	76
4	0	68	12	0	78
3	0	69	20	0	79
8	0	70	18	0	81
17	0	70

Table S-1: Challenger data. Temperature at flight time (degrees F) and failure of O-rings (1 stands for failure, 0 for success).

S-9 Examples of TMCMC for discrete state spaces

The ideas developed in this paper are not confined to continuous target distributions, but also to discrete cases. For the sake of illustration, we consider two examples below.

(i)

Consider an Ising model, where, for $i=1,\ldots,k$ $(k\geq 1)$ , the discrete random variable $x_{i}$ takes the value $+1$ or $-1$ with positive probabilities. We then have $\mathcal{X}=\{-1,1\}$ . To implement TMCMC, consider the forward transformation $T(x_{i},\epsilon)=sgn(x_{i}+\epsilon)$ with probability $p_{i}$ , and choose the backward transformation as $T^{b}(x_{i},\epsilon)=sgn(x_{i}-\epsilon)$ with probability $1-p_{i}$ . Here $sgn(a)=\pm 1$ accordingly as $a>0$ or $a<0$ , and $\mathcal{Y}=(1,\infty)$ . Note the difference with the continuous cases. Here even though neither of the transformations is 1-to-1 or onto, TMCMC works because of discreteness; the algorithm can easily be seen to satisfy detailed balance, irreducibility and aperiodicity. However, if $k=1$ with $x_{1}$ being the only variable, then, if $x_{1}=1$ , it is possible to choose, with probability one, the backward move-type, yielding $T^{b}(x_{1},\epsilon)=-1$ . On the other hand, if $x_{1}=-1$ , with probability one, we can choose the forward move-type, yielding $T(x_{1},\epsilon)=1$ . Only $2^{k}$ move-types are necessary for the $k$ -dimensional case for one-step irreducibility. In discrete cases, however, there will be no Jacobian of transformation, thereby simplifying the acceptance ratio.
(ii)

For discrete state spaces like $\mathbb{Z}^{k}$ , ( $\mathbb{Z}=\{0,\pm 1,\pm 2,\ldots\}$ ) the additive transformation with single epsilon does not work. For example, with $k=2$ , if the starting state is $(1,2)$ then the chain will never reach any states $(x,y)$ where $x$ and $y$ have same parity (i.e. both even or both odd) resulting a reducible Markov chain. Thus in this case we need to have more move-types than $2^{k}$ . For example, with some positive probability (say $r$ ) we may select a random coordinate and update it leaving other states unchanged. With the remaining probability (i.e. $1-r$ ) we may do the analogous version of the additive transformation:
Let $\mathcal{Y}=[1,\infty)$ . Then, can choose the forward transformation for each coordinate as $T_{i}(x_{i},\epsilon)=x_{i}+[\epsilon]$ and the backward transformation as $T^{b}_{i}(x_{i},\epsilon)=x_{i}-[\epsilon]$ , where $[a]$ denotes the largest integer not exceeding $a$ .

This chain is clearly ergodic and we still need only one epsilon to update the states.

However, in discrete cases, TMCMC reduces to Metropolis-Hastings with a mixture proposal. But it is important to note that the implementation is much efficient and computationally cheap when TMCMC-based methodologies developed in this paper, are used.

$\displaystyle AR$	$\displaystyle=\int R(\mathbf{x}^{\prime}\|\mathbf{x})q(\mathbf{x}^{\prime}\|\mathbf{x})\pi(\mathbf{x})d\mathbf{x}d\mathbf{x}^{\prime}$
	$\displaystyle=\int Pr\left(U<R(\mathbf{x}^{\prime}\|\mathbf{x})\right)q(\mathbf{x}^{\prime}\|\mathbf{x})\pi(\mathbf{x})d\mathbf{x}d\mathbf{x}^{\prime}$
	$\displaystyle=\int\left[\int Pr\left(U<R(\mathbf{x}^{\prime}\|\mathbf{x})\right)q(\mathbf{x}^{\prime}\|\mathbf{x})d\mathbf{x}^{\prime}\right]\pi(\mathbf{x})d\mathbf{x}$
	$\displaystyle=\int\left[\int_{0}^{1}Pr\left(R(\mathbf{x}^{\prime}\|\mathbf{x})>u\right)du\right]\pi(\mathbf{x})d\mathbf{x}$	(S-6.1)

			$\displaystyle Pr~(R(\mathbf{x}^{\prime}\|\mathbf{x})<r)=Pr~(\pi(\mathbf{x}^{\prime})/\pi(\mathbf{x})<r)=Pr~(\log\pi(\mathbf{x}^{\prime})-\log\pi(\mathbf{x})<\log r)$
		$\displaystyle=$	$\displaystyle Pr~(\left[\log\pi(\mathbf{x}^{\prime})-\log\pi(\mathbf{x^{}})\right]-\left[\log\pi(\mathbf{x}^{\prime})-\log\pi(\mathbf{x^{}})\right]<\log r)$
		$\displaystyle=$	$\displaystyle Pr~\left(\left[\nabla\log\pi(\mathbf{x^{}})^{T}(\mathbf{x}^{\prime}-\mathbf{x^{}})+(1/2)(\mathbf{x}^{\prime}-\mathbf{x^{}})^{T}\nabla^{2}\log\pi(\boldsymbol{\xi_{1}(\mathbf{x}^{\prime},\mathbf{x}^{})})(\mathbf{x}^{\prime}-\mathbf{x^{*}})\right]\right.$
			$\displaystyle\quad\quad\left.-\left[\nabla\log\pi(\mathbf{x^{}})^{T}(\mathbf{x}-\mathbf{x^{}})+(1/2)(\mathbf{x}-\mathbf{x^{}})^{T}\nabla^{2}\log\pi(\boldsymbol{\xi_{2}(\mathbf{x},\mathbf{x}^{})})(\mathbf{x}-\mathbf{x^{*}})\right]<\log r\right),$
			$\displaystyle\quad\textrm{~for some }\boldsymbol{\xi}_{1}(\mathbf{x}^{\prime},\mathbf{x}^{}),\boldsymbol{\xi}_{2}(\mathbf{x},\mathbf{x}^{})\ \ \mbox{depending upon}\ \ (\mathbf{x}^{\prime},\mathbf{x}^{})\ \ \mbox{and}(\mathbf{x},\mathbf{x}^{})\ \ \mbox{respectively};$
		$\displaystyle=$	$\displaystyle Pr~\left(\left[(1/2)(\mathbf{x}^{\prime}-\mathbf{x^{}})^{T}\nabla^{2}\log\pi(\boldsymbol{\xi_{1}(\mathbf{x}^{\prime},\mathbf{x}^{})})(\mathbf{x}^{\prime}-\mathbf{x^{*}})\right]\right.$
			$\displaystyle\quad\quad\left.-\left[(1/2)(\mathbf{x}-\mathbf{x^{}})^{T}\nabla^{2}\log\pi(\boldsymbol{\xi_{2}(\mathbf{x}^{\prime},\mathbf{x}^{})})(\mathbf{x}-\mathbf{x^{*}})\right]<\log r\right)$
			$\displaystyle\quad\quad\quad\quad\quad\mbox{since}\ \ \nabla\log\pi(\mathbf{x^{*}})=\mbox{\boldmath$0$}.$