Fast Markov Chain Monte Carlo Algorithms via Lie Groups

TL;DR

This paper introduces a new MCMC algorithm derived from Lie group principles, demonstrating faster convergence on complex models like spin glasses compared to existing methods.

Contribution

It presents a novel MCMC algorithm based on Lie group theory, expanding the toolkit for efficient sampling in complex probability spaces.

Findings

The new algorithm converges faster than traditional MCMC methods on spin glass models.

Explicit numerical examples illustrate the effectiveness of the proposed algorithms.

Derivation of multiple MCMC variants from Lie group considerations.

Abstract

From basic considerations of the Lie group that preserves a target probability measure, we derive the Barker, Metropolis, and ensemble Markov chain Monte Carlo (MCMC) algorithms, as well as variants of waste-recycling Metropolis-Hastings and an altogether new MCMC algorithm. We illustrate these constructions with explicit numerical computations, and we empirically demonstrate on a spin glass that the new algorithm converges more quickly than its siblings.

Figures

Equations

$$[X,[Y,Z]]+[Y,[Z,X]]+[Z,[X,Y]]=0.$$

$$STO(n):=\{P\in GL(n,\mathbb{R}):P1=1\}$$

$$⟨p⟩:=\{P\in STO(n):pP=p\}$$

$$e_{(j,k)}:=e_j(e_k^T-e_n^T),$$

$$({\delta }_{k\ell }-{\delta }_{n\ell })e_{(j,m)}-({\delta }_{mj}-{\delta }_{nj})e_{(\ell ,k)}.$$

$$e_{(j,k)}^{(p)}:=e_{(j,k)}-r_j{e}_{(n,k)}=\left(e_j-r_j{e}_n\right)(e_k^T-e_n^T).$$

$${\left(e_{(j,k)}^{(p)}\right)}^i=\{\begin{array}{ll}I,& i=0;\\ {\left({\delta }_{jk}+r_j\right)}^{i-1}{e}_{(j,k)}^{(p)},& i>0.\end{array}$$

$$\left[e_{(j,k)}^{(p)},e_{(\ell ,m)}^{(p)}\right]=\left({\delta }_{k\ell }+r_{\ell }\right)e_{(j,m)}^{(p)}-\left({\delta }_{mj}+r_j\right)e_{(\ell ,k)}^{(p)}.$$

$$f_{(j,k)}^{(p)}(t):=\frac{e^{-t({\delta }_{jk}+r_j)}-1}{{\delta }_{jk}+r_j}.$$

$$ST{O}^{+}(n):=\{P\in M(n,\mathbb{R}):P1=1\text{ and }P\ge 0\},$$

$$⟨p{⟩}^{+}:=\{P\in ST{O}^{+}(n):pP=p\}.$$

$$\exp \left(-t{e}_{(j,j)}^{(p)}\right)=I+f_{(j,j)}^{(p)}(t)\cdot e_{(j,j)}^{(p)}\in ⟨p{⟩}^{+}.$$

$$\frac{P_{jn}}{P_{nj}}=\frac{1}{r_j}=\frac{p_n}{p_j}.$$

$${\left(I-\min (1,r_j^{-1})\cdot e_{(j,j)}^{(p)}\right)}_{nj}=\min (1,r_j)$$

$${\gamma }_{(\mathcal{J})}={\alpha }_{(\mathcal{J})}(I+1r_{(\mathcal{J})}){\beta }_{(\mathcal{J})}.$$

$$A_{(\mathcal{J})}^{(p;\omega )}:=\omega \sum _{u,v}\left({\delta }_{j_u{j}_v}-\frac{1}{1+r_{(\mathcal{J})}1}{r}_{j_v}\right)e_{(j_u,j_v)}^{(p)}.$$

$$\exp t{A}_{(\mathcal{J})}^{(p;\omega )}=I+\frac{e^{\omega t}-1}{\omega }{A}_{(\mathcal{J})}^{(p;\omega )}.$$

$${\mathcal{B}}_{(\mathcal{J})}^{(p)}:=I-{\omega }^{-1}{A}_{(\mathcal{J})}^{(p;\omega )}$$

$${\mathcal{M}}_{(\mathcal{J})}^{(p)}:=I-\frac{1}{\max \Delta \left(A_{(\mathcal{J})}^{(p;\omega )}\right)}{A}_{(\mathcal{J})}^{(p;\omega )}$$

$$A_{(\mathcal{J})}^{(p;\omega )}=\frac{\omega }{16}\left(\begin{array}{ccccc}15& -2& -3& 0& -10\\ -1& 14& -3& 0& -10\\ -1& -2& 13& 0& -10\\ 0& 0& 0& 0& 0\\ -1& -2& -3& 0& 6\end{array}\right).$$

$$\exp \left(\log 2\cdot A_{(\mathcal{J})}^{(p;1)}\right)=\frac{1}{32}\left(\begin{array}{ccccc}17& 2& 3& 0& 10\\ 1& 18& 3& 0& 10\\ 1& 2& 19& 0& 10\\ 0& 0& 0& 32& 0\\ 1& 2& 3& 0& 26\end{array}\right).$$

$${\mathcal{B}}_{(\mathcal{J})}^{(p)}=\frac{1}{16}\left(\begin{array}{ccccc}1& 2& 3& 0& 10\\ 1& 2& 3& 0& 10\\ 1& 2& 3& 0& 10\\ 0& 0& 0& 16& 0\\ 1& 2& 3& 0& 10\end{array}\right);{\mathcal{M}}_{(\mathcal{J})}^{(p)}=\frac{1}{15}\left(\begin{array}{ccccc}0& 2& 3& 0& 10\\ 1& 1& 3& 0& 10\\ 1& 2& 2& 0& 10\\ 0& 0& 0& 15& 0\\ 1& 2& 3& 0& 9\end{array}\right).$$

$$\frac{1}{\omega }{\left(A_{(\mathcal{J})}^{(p;\omega )}\right)}_{j_u{j}_u}$$

$$\frac{1}{\omega }{\left(A_{(\mathcal{J})}^{(p;\omega )}\right)}_{n{j}_u}$$

$$\frac{1}{\omega }{\left(A_{(\mathcal{J})}^{(p;\omega )}\right)}_{nn}$$

$${\left({\mathcal{B}}_{(\mathcal{J})}^{(p)}\right)}_{n{j}_u}$$

$$\frac{1}{\omega }\max \Delta \left(A_{(\mathcal{J})}^{(p;\omega )}\right)=\frac{1+r_{(\mathcal{J})}1-\min \{1,\min r_{(\mathcal{J})}\}}{1+r_{(\mathcal{J})}1}$$

$${\left({\mathcal{M}}_{(\mathcal{J})}^{(p)}\right)}_{n{j}_u}$$

$${\left({\mathcal{M}}_{(\mathcal{J})}^{(p)}\right)}_{nn}$$

$$p(s):=Z^{-1}\exp \left(-\frac{\beta }{\sqrt{N}}{\sum }_{jk}{J}_{jk}{s}_j{s}_k\right)$$

$$(1_{\mathcal{J}}{)}_j:=\{\begin{array}{ll}1& \text{if }j\in \mathcal{J}\cup \{n\}\\ 0& \text{otherwise},\end{array}$$

$${\tau }_{(\mathcal{J})}^{(p)}=\left(\begin{array}{c}{I}_{n-1}\\ -r_{\mathcal{J}}^{-}\end{array}\right)\tau \left(\begin{array}{cc}{I}_{n-1}& -1_{\mathcal{J}}^{-}\end{array}\right),$$

Full text

1 Introduction

The basic problem that Markov chain Monte Carlo (MCMC) algorithms solve is to estimate expected values using a Markov chain that has the desired probability measure as its invariant measure. Originally developed to solve problems in computational physics, MCMC algorithms have since become omnipresent in statistical inference and machine learning. Indeed, it is reasonable to suggest in line with Richey (2010) and Brooks et al. (2011) that MCMC algorithms comprise the most ubiquitous and important class of high-level numerical algorithms discovered to date.

Consequently, the literature on MCMC algorithms is vast. However, the mostly unexplored interface of MCMC algorithms and the theory of Lie groups and Lie algebras holds a surprise. As we shall see, the space of transition matrices with a given invariant measure is a monoid that is closely related to a Lie group. Searching for elements of this monoid with closed form expressions naturally leads to the classical Barker and Metropolis MCMC samplers. Generalizing this search leads to higher-order versions of these samplers which respectively correspond to the ensemble MCMC algorithm of Neal (2011) and an algorithm of Delmas and Jourdain (2009). Further generalizing this search leads to an algorithm which we call the higher-order programming solver and whose convergence appears to improve on the state of the art. For each of these algorithms (only treated for finite state spaces), the acceptance mechanism is specified, but the proposal mechanism is not (though one can always extend any proposal mechanism for single states to multiple states by repeated sampling).

In this paper, we first review the basics of MCMC, Lie theory, and related work in §2. Next, we introduce the Lie group generated by a probability measure in §3. Here, Lemma 1 exhibits a convenient basis of the stochastic Lie algebra. Lemma 2 and Theorem 1 next yield a convenient basis of the Lie subalgebra that annihilates a target probability measure $p$. Critically, this basis only requires knowledge of $p$ up to a multiplicative factor. In §4, we consider a closely related monoid, and Lemma 3 shows how we can analytically produce nonnegative transition matrices that leave $p$ invariant. We then exhibit the construction of the Barker and Metropolis samplers from Lie-theoretic considerations in §5. In §6, Lemma 4 extends Lemma 2, and Theorem 2 extends Theorem 1 in such a way as to yield generalizations of the preceding samplers that entertain multiple proposals at once. These higher-order Barker and Metropolis samplers are explicitly constructed in §7. We then demonstrate their behavior on a simple example of a spin glass in §8. In §9, Theorem 3 yields multiple-proposal transition matrices that are closest in Frobenius norm to the “ideal” transition matrix $1p$: in this section we introduce and demonstrate the resulting higher-order programming solver. Finally, we close with remarks in §10. Proofs, though brief, are relegated to Appendix A.

2 Background

2.1 Markov chain Monte Carlo

As mentioned in §1 and Brémaud (1999), the basic problem of MCMC is to estimate expected values of functions with respect to a probability measure $p$ that is infeasible to construct. A common instance is where $p_{j}\equiv\mathcal{L}_{j}/Z$, where it is easy to compute $\mathcal{L}$ but hard to compute the normalizing constant $Z$ due to the scale of the problem. The approach of MCMC is to construct an irreducible, ergodic Markov chain that has $p$ as its invariant measure without using global information.

If now $X_{t}$ is the state of such a chain at time $t$, then in the limit we have $X_{t}\sim p$ for any initial condition. For $f$ suitable, $\mathbb{E}_{p}f(X)=\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{j=1}^{t}f(X_{j})$ even though the $X_{j}$ are correlated.

A MCMC algorithm typically depends on respective proposal and acceptance probabilities $q_{jk}:=\mathbb{P}(X^{\prime}=k|X_{t}=j)$ and $\alpha_{jk}:=\mathbb{P}(X_{t+1}=k|X^{\prime}=k,X_{t}=j)$, which yield $P_{jk}:=\mathbb{P}(X_{t+1}=k|X_{t}=j)=q_{jk}\alpha_{jk}$ for the elements of the chain transition matrix.

The Hastings algorithm uses a symmetric matrix $s$ to accept a proposal with probability $\alpha_{jk}=\frac{s_{jk}}{1+t_{jk}}$, where $t_{jk}:=\frac{p_{j}q_{jk}}{p_{k}q_{kj}}$. This requires $s_{jk}\leq 1+\min(t_{jk},t_{kj})$. Taking $s_{jk}=1$ yields the Barker sampler; Peskun (1973) shows that the optimal choice $s_{jk}=1+\min(t_{jk},t_{kj})$ yields the Metropolis-Hastings sampler.

2.2 Lie groups and Lie algebras

For general background on Lie groups and algebras, we refer to Onishchik and Vinberg (1990) and Kirillov (2008). Here, we briefly restate the basic concepts in the real and finite-dimensional setting.

A Lie group is a group that is also a manifold, and for which the group operations are smooth. The tangent space of a Lie group $G$ at the identity is a Lie algebra that we denote by $\mathfrak{lie}(G)$. Besides its vector space structure, this Lie algebra inherits a version of the Lie group structure through a bilinear antisymmetric bracket $[\cdot,\cdot]$ that satisfies the Jacobi identity

[TABLE]

In particular, Ado’s theorem implies that a real finite-dimensional Lie group is isomorphic to a subgroup of the group $GL(n,\mathbb{R})$ of invertible $n\times n$ matrices over $\mathbb{R}$. Meanwhile, the corresponding Lie algebra is isomorphic to a Lie subalgebra of real $n\times n$ matrices, for which the bracket is the usual matrix commutator: $[X,Y]:=XY-YX$. In the other direction, the usual matrix exponential gives a map from a matrix Lie algebra to the corresponding Lie group that respects both the algebra and group structures.

2.3 Related work

Our higher-order Barker and Metropolis samplers respectively correspond to constructions in Neal (2011) and Delmas and Jourdain (2009). Besides ensemble algorithms, Robert et al. (2018) details a large body of work on accelerating MCMC algorithms by techniques such as multiple try algorithms as in Liu et al. (2000), Martino (2018), and Martino et al. (2018); or by parallelization, as in Calderhead (2014).

Work by Niepert (2012a), Niepert (2012b), Bui et al. (2013), Shariff et al. (2015), Van den Broeck and Niepert (2015) and Anand et al. (2016) dealt with accelerating MCMC algorithms by exploiting discrete symmetries that preserve the (exact or approximate) level sets of a target measure. In a related vein, “group moves” for MCMC algorithms were considered by Liu and Wu (1999) and Liu and Sabatti (2000). There is also a long tradition of learning and exploiting symmetries in data representations for machine learning (see, e.g. Lüdtke et al. (2018), and Anselmi et al. (2019)): including for neural networks, for which see Cohen and Welling (2016) and Cohen et al. (2018). However, to our knowledge, the present paper is the first attempt to consider continuous symmetries that preserve a target measure in the context of MCMC.

The study of Markov models on groups has been studied in considerable depth, as in Saloff-Coste (2001) and Ceccherini-Silberstein et al. (2008). However, although the idea of applying Lie theory to Markov models motivates work on the stochastic group, actual applications themselves are few and far between, with Sumner et al. (2012) serving as an exemplar.

If we sacrifice analytical tractability and/or computational convenience, it is possible to consider generic MCMC algorithms that optimize some criterion over the relevant monoid. Optimal control considerations lead to algorithms such as those of Suwa and Todo (2010), Chen and Hwang (2013), Bierkens (2016), and Takahashi and Ohzeki (2016) that optimize convergence at the cost of reversibility/detailed balance. Alternatively, Frigessi et al. (1992), Pollet et al. (2004), Chen et al. (2012), Wu and Chu (2015), and Huang et al. (2012) try to optimize the asymptotic variance.

3 The Lie group generated by a measure

For $1<n\in\mathbb{N}$, let $p$ be a probability measure on $[n]:=\{1,\dots,n\}$. Relying on context to resolve ambiguity, write $1=(1,\dots,1)^{T}\in\mathbb{R}^{n}$. Following Johnson (1985), Poole (1995), Boukas et al. (2015), and Guerra and Sarychev (2018), define the stochastic group

[TABLE]

as the stabilizer fixing $1$ on the left in $GL(n,\mathbb{R})$, and

[TABLE]

as the stabilizer fixing $p$ on the right: we call $\langle p\rangle$ the group generated by $p$. $STO(n)$ and $\langle p\rangle$ are Lie groups of respective dimension $n(n-1)$ and $(n-1)^{2}$.

If $P\in STO(n)$ is irreducible and ergodic, then it has a unique invariant measure that we write as $\langle P\rangle:=1^{T}(P-I+11^{T})^{-1}$, so that $pP=p$ iff $\langle p\rangle=\langle\langle P\rangle\rangle$. We have that $\langle p\rangle-I\subset\mathfrak{lie}(\langle p\rangle)\subset\mathfrak{lie}(STO(n))$.

For $(j,k)\in[n]\times[n-1]$, define

[TABLE]

where $\{e_{j}\}_{j\in[n]}$ is the standard basis of $\mathbb{R}^{n}$.

Lemma 1**.**

The matrices $\{e_{(j,k)}\}_{(j,k)\in[n]\times[n-1]}$ form a basis of $\mathfrak{lie}(STO(n))$ and $\left[e_{(j,k)},e_{(\ell,m)}\right]$ equals

[TABLE]

This basis has the obvious advantage of computationally trivial decompositions.

For $j,k\in[n-1]$, define $r_{j}:=p_{j}/p_{n}$ and

[TABLE]

Observe that if $p_{j}\equiv\mathcal{L}_{j}/Z$, then $r_{j}=\mathcal{L}_{j}/\mathcal{L}_{n}$ does not depend on $Z$. This is why MCMC methods allow us to avoid computing such normalization factors, which in turn is why MCMC methods are useful.

For future reference, we define $r:=(r_{1},\dots,r_{n-1},1)$ and $r^{-}:=(r_{1},\dots,r_{n-1})$.

Lemma 2**.**

For $i\in\mathbb{N}$,

[TABLE]

Theorem 1**.**

The $e_{(j,k)}^{(p)}$ form a basis for $\mathfrak{lie}(\langle p\rangle)$ and

[TABLE]

For later convenience, we write

[TABLE]

4 The positive monoid of a measure

Most of the elements of $STO(n)$ are not bona fide stochastic matrices because they have negative entries; meanwhile, stochastic matrices need not be invertible. We therefore consider the monoids (i.e., semigroups with identity; cf. Hilgert and Neeb (1993))

[TABLE]

where $P\geq 0$ is interpreted per entry, and

[TABLE]

Note that $STO^{+}(n)\not\subset STO(n)$ and $\langle p\rangle^{+}\not\subset\langle p\rangle$ owing to the noninvertible elements on the LHSs. Also, $STO^{+}(n)$ and $\langle p\rangle^{+}$ are bounded convex polytopes.

Lemma 3**.**

If $t_{j}\geq 0$, then $\exp\left(-\sum_{j}t_{j}e_{(j,j)}^{(p)}\right)\in\langle p\rangle^{+}$.

In particular, for $t\geq 0$ we have that

[TABLE]

Unfortunately, aside from (10), Lemma 3 does not give a way to construct explicit elements of $\langle p\rangle^{+}$ in closed form, or even algorithmically. This situation is an analogue of the quantum compilation problem (see Dawson and Nielsen (2006)), which is by no means trivial.

Indeed, even if the sum in the lemma’s statement has only two terms, we are immediately confronted with the formidable Zassenhaus formula (see Casas et al. (2012)). While $\exp\left(-t_{(j,k)}e_{(j,k)}^{(p)}-t_{(\ell,m)}e_{(\ell,m)}^{(p)}\right)$ can be evaluated in closed form with a computer algebra package, the results involve many pages of arithmetic for the case corresponding to Lemma 3, and the other possibilities all yield some negative entries.

5 Barker and Metropolis samplers

Despite the weak foothold that Lemma 3 affords for explicit analytical constructions, we can still use (10) to produce a MCMC algorithm parametrized by $t$. We use a simple trick of relabeling the current state as $n$ and then reversing the relabeling, so that the transition $n\rightarrow j$ becomes generic. For $P=\exp\left(-te_{(j,j)}^{(p)}\right)$, we have $P_{jj}=1+f_{(j,j)}^{(p)}(t)$, $P_{jn}=-f_{(j,j)}^{(p)}(t)$, $P_{nj}=-f_{(j,j)}^{(p)}(t)r_{j}$, and $P_{nn}=1+f_{(j,j)}^{(p)}(t)r_{j}$. In particular,

[TABLE]

That is, detailed balance is automatic.

From the point of view of convergence, the optimal value for $t$ is the one that maximizes the off-diagonal terms, i.e., $t=\infty$. Here we get $P_{jj}=\frac{r_{j}}{1+r_{j}}$, $P_{jn}=\frac{1}{1+r_{j}}$, $P_{nj}=\frac{r_{j}}{1+r_{j}}$, and $P_{nn}=\frac{1}{1+r_{j}}$. The corresponding MCMC algorithm is a Barker sampler.

In light of (10), we can improve on the Barker sampler almost trivially. We have that $I-\tau e_{(j,j)}^{(p)}\in\langle p\rangle^{+}$ iff $0\leq\tau\leq\min(1,r_{j}^{-1})$. But

[TABLE]

is just the Metropolis acceptance ratio. That is, we have derived the Barker and Metropolis samplers from basic considerations of symmetry and (in the latter case) optimality.

Note that the proposal mechanism that selects the state $j$ is unspecified and unconstrained by our construction (i.e., $\mathcal{J}$ can be drawn from an arbitrary joint distribution on the subset of $[n-1]^{d}$ without duplicate entries). That is, our approach separates concerns between the proposal and acceptance mechanisms, and only focuses on the latter. In later sections of this paper, the proposal mechanism that selects a set of states will similarly be unspecified and unconstrained. However, in §8 we select the elements of proposal sets uniformly at random without replacement for illustrative purposes. That said, a good proposal mechanism is of paramount importance for MCMC algorithms.

6 Some algebra

The Barker and Metropolis samplers can be regarded as among the very “simplest” MCMC methods in the sense that (10) is among the very sparsest possible nontrivial matrices in $\langle p\rangle^{+}$. This suggests the question: what happens if we are willing to sacrifice some sparsity? In other words, what if we consider possible transitions to more than one state? It is natural to expect both better convergence and increased complexity. The (utterly impractical and degenerate) limiting case is the matrix $1p$, and the practical starting case is the Barker and Metropolis samplers. Meanwhile, it is also natural to wonder how (or if) we can analytically construct more general elements of $\langle p\rangle^{+}$ than (10).

The following generalization of Lemma 2 is the first step toward an answer to the preceding questions. For $\mathcal{J}:=\{j_{1},\dots,j_{d}\}\subseteq[n-1]$ and a matrix $\alpha\in M_{n-1,n-1}$, define $\alpha_{(\mathcal{J})}\in M_{d,d}$ by $(\alpha_{(\mathcal{J})})_{uv}:=\alpha_{j_{u}j_{v}}$, $\alpha_{(\mathcal{J})}^{(p)}:=\sum_{u,v=1}^{d}\alpha_{j_{u}j_{v}}e_{(j_{u},j_{v})}^{(p)}\in\mathfrak{lie}(\langle p\rangle)$, and $r_{(\mathcal{J})}:=(r_{j_{1}},\dots,r_{j_{d}})$.

Lemma 4**.**

Let $\mathcal{J}:=\{j_{1},\dots,j_{d}\}\subseteq[n-1]$. If $\gamma_{(\mathcal{J})}^{(p)}=\alpha_{(\mathcal{J})}^{(p)}\beta_{(\mathcal{J})}^{(p)}$, then

[TABLE]

We remark that introducing heavy notation for Lemma 4 is worth it: the proof takes just three lines, whereas the case $d=2$ takes about a page of algebra to check otherwise. Using Lemma 4, we can readily construct an analytically convenient matrix in $\mathfrak{lie}(\langle p\rangle)$.

Theorem 2**.**

Let $\mathcal{J}:=\{j_{1},\dots,j_{d}\}\subseteq[n-1]$, $\omega\in\mathbb{R}$ and

[TABLE]

Then

[TABLE]

Moreover, $\exp\left(-tA_{(\mathcal{J})}^{(p;\omega)}\right)\in\langle p\rangle^{+}\cap GL(n,\mathbb{R})$ if $t\geq 0$. In particular, the Barker matrix

[TABLE]

is in $\langle p\rangle^{+}$, and does not depend on $\omega$.

Let $\Delta$ denote the map that takes a matrix to the vector of its diagonal entries, and indicate the boundary of a set using $\partial$.

Lemma 5**.**

The Metropolis matrix

[TABLE]

is in $\partial\langle p\rangle^{+}$ and does not depend on $\omega$.

6.1 Example

As an example, consider $p=(1,2,3,4,10)/20$ and $\mathcal{J}=\{1,2,3\}$. Now (13) is given by

[TABLE]

For $\omega=1$ and $t=-\log 2$, (14) is given by

[TABLE]

Finally, (15) and (16) are respectively given by

[TABLE]

7 Higher-order samplers

The idea now is to let $n\rightarrow j\in\mathcal{J}$ correspond to a generic transition as in §5. (Again, we do not specify or constrain a proposal that produces $\mathcal{J}$.) This yields novel MCMC algorithms using (15) and (16) which we respectively call higher-order Barker and Metropolis samplers and abbreviate as HOBS and HOMS.

The corresponding matrix elements are readily obtained with a bit of arithmetic: we have that

[TABLE]

which yields the HOBS:

[TABLE]

Meanwhile,

[TABLE]

which yields the HOMS:

[TABLE]

The HOBS turns out to be equivalent to the ensemble MCMC algorithm of Neal (2011) as described in Martino (2018) and Martino et al. (2018). While in §8 the proposal mechanism we use for the HOBS essentially (apart from non-replacement, which technically induces jointness) amounts to the independent ensemble MCMC sampler, in general this is not the case. A more sophisticated proposal mechanism with joint structure will be more powerful. That said, we reiterate that our approach is agnostic with respect to the details of proposals.

On the other hand, the HOMS is different than a multiple-try Metropolis sampler (MTMS), including the independent MTMS described in Martino (2018). In the HOMS, we sample from $\mathcal{J}\cup\{n\}$ to perform a state transition in a single step according to (7), whereas a MTMS first samples from $\mathcal{J}$ and then accepts or rejects the result. The HOMS (and for that matter, also the HOBS) turns out to be a slightly special case of a construction in §2.3 of Delmas and Jourdain (2009). This work uses a “proposition kernel” that is defined by assigning a probability distribution on the power set of the state space to each element of the state space. Roughly speaking, the HOMS and HOBS result if this distribution is independent of the individual element (i.e., varying only with the subset).

8 Behavior

As $d=|\mathcal{J}|$ increases and/or $p$ becomes less uniform (e.g., in a low-temperature limit), the difference between the HOBS and HOMS decreases, since in either limit we have $\min\{1,\min r_{(\mathcal{J})}\}\ll 1+r_{(\mathcal{J})}1$. Although these limits are where one might hope to gain the most utility from improved MCMC algorithms, the HOMS can still provide an advantage in, e.g. the high-temperature part of a parallel tempering scheme (see Earl and Deem (2005)) or for $d>1$ but small, with elements chosen in complementary ways (uniformly at random, near current/previous states, etc.).

We exhibit the the behavior of the HOBS and HOMS on a Sherrington-Kirkpatrick (SK) spin glass in Figure 1. As Bolthausen and Bovier (2007) and Panchenko (2012) remark, the SK spin glass is the distribution

[TABLE]

over spins $s\in\{\pm 1\}^{N}$, where $J$ is a symmetric $N\times N$ matrix with IID standard Gaussian entries and $\beta$ is the inverse temperature.

The disordered landscape of the SK model suits a straightforward evaluation of higher-order samplers: more detailed benchmarks seem to require specific assumptions (e.g., exploiting the particular form of a spin Hamiltonian for Swendsen-Wang updates) and/or parameters (e.g., the choice of additional temperatures for parallel tempering or of a vorticity matrix for non-reversible Metropolis-Hastings). In particular, we do not consider sophisticated or diverse ways to generate elements of proposal sets $\mathcal{J}$: instead, we merely select the elements of $\mathcal{J}$ uniformly at random without replacement. We also use the same pseudorandom number generator initial state for the HOBS and HOMS simulations in order to highlight their relative behavior, and pick $\beta$ low enough ($1/4$ and $1$) so that the behavior of a single run is sufficiently representative to make qualitative judgments.

The inset figure shows that although higher-order samplers converge more quickly, this requires more evaluations of probability ratios. Parallelism is therefore necessary to make higher-order samplers worthwhile.

As noted above, the HOMS gives results very close to the HOBS except for small values of $d$ or a more uniform target distribution $p$. Increasing the number of spins and/or considering an Edwards-Anderson spin glass also gives qualitatively similar results.

9 A linear program

We can push these ideas further by using a linear program to (implicitly) construct transition matrices with the desired invariant measure and that are optimal in some sense, though the regime of utility then narrows to situations where computing likelihoods is hard enough and parallel resources are sufficient to justify the added computational costs. For example, the particular objective function $-1_{\mathcal{J}}^{T}\tau_{(\mathcal{J})}^{(p)}r_{\mathcal{J}}^{T}$ considered immediately after (25) yields an optimal sparse approximation of the “ultimate” transition matrix $1p$. (To the best of our knowledge, this construction has not been considered elsewhere.)

Toward this end, define $1_{\mathcal{J}}\in\mathbb{R}^{n}$ by

[TABLE]

$1_{\mathcal{J}}^{-}:=((1_{\mathcal{J}})_{1},\dots,(1_{\mathcal{J}})_{n-1})^{T}$, $r_{\mathcal{J}}:=r\odot 1_{\mathcal{J}}^{T}$, and $r_{\mathcal{J}}^{-}:=r^{-}\odot(1_{\mathcal{J}}^{-})^{T}$, where $\odot$ is the entrywise or Hadamard product (note that $r_{\mathcal{J}}\in\mathbb{R}^{n}$, while $r_{(\mathcal{J})}\in\mathbb{R}^{|\mathcal{J}|}$ has been defined previously).

Writing $\Delta$ for the matrix diagonal map, using the notation of Lemma 4, and noting that

[TABLE]

we have that $I-\tau_{(\mathcal{J})}^{(p)}\in\langle p\rangle^{+}$ iff

[TABLE]

Constraints (22b)-(22d) respectively force the first $n-1$ entries of the last column, the first $n-1$ entries of the last row, and the bottom right matrix entry of $\tau_{(\mathcal{J})}^{(p)}$ to be in the unit interval; (22a) forces the relevant entries of the “coefficient matrix” $\tau$ (as an upper left submatrix of $\tau_{(\mathcal{J})}^{(p)}$) to be in the unit interval.

Furthermore, it is convenient to set to zero the irrelevant/unspecified rows and columns of $\tau$ that do not contribute to $\tau_{(\mathcal{J})}^{(p)}$ via the constraints

[TABLE]

If we impose (23), then (22a) can be replaced with

[TABLE]

The “diagonal” case corresponding to Lemma 3 shows that the constraints (22) and (23) jointly have nontrivial solutions. It is therefore natural to consider suitable objectives and the corresponding linear programs for optimizing the MCMC transition matrix $I-\tau_{(\mathcal{J})}^{(p)}$. Toward this end, we introduce the vectorization map vec that sends a matrix to the vector obtained by stacking the matrix columns in order, and which obeys the useful identity $\text{vec}(XYZ^{T})=(Z\otimes X)\text{vec}(Y)$, where $\otimes$ denotes the tensor product.

A reasonably generic objective to maximize is

[TABLE]

for suitable vectors $x$ and $y$. In practice, we shall take $x=1_{\mathcal{J}}$ and $y=-r_{\mathcal{J}}^{T}$, so that our objective maximizes the Frobenius inner product of $I-\tau_{(\mathcal{J})}^{(p)}$ and $1_{\mathcal{J}}r_{\mathcal{J}}$ as a consequence of the equality

[TABLE]

Alternatives such as $x=e_{n},y=e_{n}$ (which discourages self-transitions) can result in convergence that slows catastrophically as $d=|\mathcal{J}|$ increases, because high-probability states are less likely to remain occupied. More surprisingly, the same sort of slowing down occurs for $x=e_{n},y=-r_{\mathcal{J}}^{T}$ and even for variations upon the $n$th component of $y$: we suspect that the cause is the same, though mediated indirectly through an objective that “overfits” the proposed transition probabilities to the detriment of remaining in place (or in some cases “underfits” by yielding the identity matrix). In general, it appears nontrivial to select better choices for $x$ and $y$ than our defaults.

By (21) we get

[TABLE]

and in turn $(y^{T}\otimes x^{T})\text{vec}\left(\tau_{(\mathcal{J})}^{(p)}\right)$ equals

[TABLE]

At this point both the constraints and the objective of the linear program are explicitly specified in terms of the “coefficient” matrix $\tau$. We can rephrase the constraints into a more computationally convenient form, respectively rephrasing (22b)-(22d), (23), and (24) as

[TABLE]

[TABLE]

[TABLE]

Therefore, writing

[TABLE]

and

[TABLE]

we can at last write the desired linear program (noting the inclusion of a minus sign in $w_{(\mathcal{J})}^{(p)}$ and a minimization versus a maximization as a result) in the MATLAB-ready form

[TABLE]

As a result of the preceding discussion, we have

Theorem 3**.**

For any $x,y\in\mathbb{R}^{n}$, the linear program (32) has a solution in $\langle p\rangle^{+}$. ∎

9.1 Example

As in §6.1, consider $p=(1,2,3,4,10)/20$ and $\mathcal{J}=\{1,2,3\}$. The solution of the linear program with $x=1_{\mathcal{J}}$ and $y=-r_{\mathcal{J}}^{T}$ yields the following element of $\langle p\rangle^{+}$:

[TABLE]

For comparison, we recall that the last row of $\mathcal{M}_{(\mathcal{J})}^{(p)}$ equals $(0.0\bar{6},0.1\bar{3},0.2,0,0.6)$.

9.2 The higher-order programming sampler

Call the sampler obtained from (25) and (32) with $x=-1_{\mathcal{J}}$ and $y=r_{\mathcal{J}}^{T}$ the higher-order programming sampler (HOPS). In figures 2 and 3, we compare the HOMS and HOPS (cf. Figure 1). It is clear from the figures that the HOPS improves upon the HOMS, which in turn improves upon the HOBS.

10 Remarks

Aside from providing a framework that unifies several different MCMC algorithms, our perspective has uncovered the apparently new HOPS algorithm of §9, which may enhance existing MCMC techniques specifically tailored for parallel computation, as in Conrad et al. (2018). In particular, the Bayesian approach to inverse problems detailed in Dashti and Stuart (2015) offers fertile ground for useful applications.

As we have already indicated, the present paper is agnostic with respect to proposals, and focuses on acceptance mechanisms. However, the proposal arguably plays a more important role in practice than the acceptance mechanism, particularly for differentiable distributions. In any practical application, a stateful and/or problem-specific proposal mechanism with joint structure would likely confer significant additional power to our approach, though we leave investigations along these lines open for now (one possibility is suggested by particle MTMS algorithms as in Martino (2014) and exploiting tensor product structure in transition matrices and $\langle p\rangle$). It is also tempting to try to incorporate some limited proposal mechanism into the objective of (32), but it is not clear how to usefully do this in general. We used the SK spin glass to illustrate our ideas precisely because its highly disordered structure (and discrete state space) are suited for separating concerns about proposals and acceptance.

It would obviously be interesting to extend the considerations of this paper to continuous variables. However, this seems to require a much more technical treatment, as infinite-dimensional Lie theory, distributions, etc. would inevitably arise at least in principle. We leave this for future work. In a complementary vein, it would be interesting to see if the full construction of Delmas and Jourdain (2009) can be recovered from considerations of symmetry/Lie theory alone.

While Barker and Metropolis samplers are reversible, it is not clear if the HOPS is, though Bierkens (2016) points out ways to transform reversible kernels into irreversible ones and vice versa.

We note that is possible to produce transiton matrices (as it turns out, even in closed form) in which the $n$th row is nonnegative but other rows have negative entries. It is not immediately clear if using such a matrix inevitably poisons a MCMC algorithm. Though our experiments in this direction were not encouraging, we have not found a compelling argument that rules out the use of such matrices.

It is tempting to try to sample from the vertices of the polytope $\langle p\rangle^{+}$. However, (even approximately) uniformly sampling vertices of a polytope is $\mathbf{NP}$-hard by Theorem 1 of Khachiyan (2001); see also Khachiyan et al. (2008).

Acknowledgements

We thank BAE Systems FAST Labs for its support; Carlo Beenakker, Jun Liu, Lorenzo Najt, Allyson O’Brien, and Daniel Zwillinger for helpful comments; and reviewers for their careful efforts, especially for bringing our attention to Delmas and Jourdain (2009).

Appendix A Proofs

Proof of Lemma 1..

$e_{(j,k)}e_{(\ell,m)}=(\delta_{k\ell}-\delta_{n\ell})e_{(j,m)}$. Considering $j\leftrightarrow\ell,k\leftrightarrow m$, we are done. ∎

Proof of Lemma 2..

Using the rightmost expression in (5) and using $j,k,\ell,m\neq n$ to simplify the product of the innermost two factors, we have that

[TABLE]

Taking $j=\ell$ and $k=m$ establishes the result for $i\leq 2$. The general case follows by induction on $i$. ∎

Proof of Theorem 1..

Note that

[TABLE]

Furthermore, linear independence and the commutation relations are obvious, so it suffices to show that $\exp te_{(j,k)}^{(p)}\in\langle p\rangle$ for all $t\in\mathbb{R}$. By Lemma 2,

[TABLE]

∎

Proof of Lemma 3..

By hypothesis and (5), $-\sum_{j}t_{j}e_{(j,j)}^{(p)}$ has nonpositive diagonal entries and nonnegative off-diagonal entries (i.e., it is a generator matrix for a continuous-time Markov process); the result follows. ∎

Proof of Lemma 4..

[TABLE]

where the second equality follows from (A) and the third from bookkeeping. ∎

Proof of Theorem 2..

The Sherman-Morrison formula (see Horn and Johnson (2013)) gives that

[TABLE]

and the elements of this matrix are precisely the coefficients in (13). Using the notation of Lemma 4, we can therefore rewrite (13) as

[TABLE]

whereupon invoking the lemma itself yields $\left(A_{(\mathcal{J})}^{(p;\omega)}\right)^{i+1}=\omega^{i}A_{(\mathcal{J})}^{(p;\omega)}$ for $i\in\mathbb{N}$. The result now follows similarly to Theorem 1. ∎

Proof of Lemma 5..

Writing $A\equiv A_{(\mathcal{J})}^{(p;\omega)}$ here for clarity, the result follows from three elementary observations: $\Delta(A)\geq 0$, $\max\Delta\left(A\right)>0$, and $A-\Delta\left(\Delta\left(A\right)\right)\leq 0$. ∎