A quantum parallel Markov chain Monte Carlo

Parallel MCMC techniques use multiple proposals to obtain efficiency gains over MCMC algorithms such as Metropolis-Hastings (Metropolis et al., 1953; Hastings, 1970) and its progeny that use only a single proposal. Neal (2003) first develops efficient MCMC transitions for inferring the states of hidden Markov models by proposing a ‘pool’ of candidate states and using dynamic programming to select among them. Next, Tjelmeland (2004) considers inference in the general setting and shows how to maintain detailed balance for an arbitrary number $P$ of proposals. Consider a probability distribution $\pi(\mbox{d}\bm{\theta})$ defined on $\mathbb{R}^{D}$ that admits a probability density $\pi(\bm{\theta})$ with respect to the Lebesgue measure, i.e., $\pi(\mbox{d}\bm{\theta})=:\pi(\bm{\theta})\mbox{d}\bm{\theta}$. To generate samples from the target distribution $\pi$, we craft a kernel $P(\bm{\theta}_{0},\mbox{d}\bm{\theta})$ that satisfies

for all $A\subset\mathbb{R}^{D}$. Letting $\bm{\theta}_{0}$ denote the current state of the Markov chain, Tjelmeland (2004) proposes sampling from such a kernel $P(\bm{\theta}_{0},\cdot)$ by drawing $P$ proposals $\bm{\Theta}_{-0}=(\bm{\theta}_{1},\dots,\bm{\theta}_{P})$ from a joint distribution $Q(\bm{\theta}_{0},\mbox{d}\bm{\Theta}_{-0})=:q(\bm{\theta}_{0},\bm{\Theta}_{-0})\mbox{d}\bm{\Theta}_{-0}$ and selecting the next Markov chain state from among the current and proposed states with probabilities

Here, if $\bm{\Theta}=(\bm{\theta}_{0},\dots,\bm{\theta}_{P})$ is the $D\times(P+1)$ matrix having the current state and all $P$ proposals as columns, then $\bm{\Theta}_{-p}$ is the $D\times P$ matrix obtained by removing the $p$th column. Tjelmeland (2004) shows that the kernel $P(\bm{\theta}_{0},\cdot)$ constructed in such a manner maintains detailed balance and hence satisfies (1). Others have since built on this work, developing parallel MCMC methods that generate or recycle multiple proposals (Frenkel, 2004; Delmas and Jourdain, 2009; Calderhead, 2014; Yang et al., 2018; Luo and Tjelmeland, 2019; Schwedes and Calderhead, 2021; Holbrook, 2021). Most recently, Glatt-Holtz et al. (2022) place all these developments in their natural measure theoretic context, allowing one to trivially apply (1) and (2) to distributions over discrete-valued random variables.

Taken together, these developments demonstrate the ability of parallel MCMC methods to alleviate inferential challenges such as multimodality and to deliver performance gains over single-proposal competitors as measured by reduced autocorrelation between samples. In the following, we focus on the specification presented in Algorithm 1 but note that the techniques we present may also be effective for generalizations of this basic algorithm.

One need not use parallel computing to implement Algorithm 1, but the real promise and power of parallel MCMC comes from its natural parallelizability (Calderhead, 2014). Contemporary hardware design emphasizes architectures that enable execution of multiple mathematical operations simultaneously. Parallel MCMC techniques stand to leverage technological developments that keep modern computation on track with Moore’s Law, which predicts that processing power doubles every two years. For example, the algorithm of Tjelmeland (2004) generates $P$ conditionally independent proposals and then evaluates the probabilities of (2). One may parallelize the proposal generation step using parallel pseudorandom number generators (PRNG) such as those advanced in Salmon et al. (2011). The computational complexity of the target evaluations $\pi(\bm{\theta}_{p})$ is linear in the number of proposals. This presents a significant burden when $\pi(\cdot)$ is computationally expensive, e.g., in big data settings or for Bayesian inversion, but evaluation of the target density over the $P$ proposals is again a naturally parallelizable task. Moreover, widely available machine learning software such as TensorFlow allows users to easily parallelize both random number generation and target evaluations on general purpose graphics processing units (GPU) (Lao et al., 2020). Finally, when generating independent proposals using a proposal distribution of the form $q(\bm{\theta}_{0},\bm{\Theta}_{-0})=\prod_{p=1}^{P}q(\bm{\theta}_{0},\bm{\theta}_{p})$, the acceptance probabilities (2) require the ${\cal O}\left(P^{2}\right)$ evaluation of the $P+1\choose 2$ terms $q(\bm{\theta}_{p},\bm{\theta}_{p^{\prime}})$, but Holbrook et al. (2021a, b) demonstrate the natural parallelizability of such pairwise operations, obtaining orders-of-magnitude speedups with contemporary GPUs. The proposed method directly addresses the acceptance step of Algorithm 1, while leaving the choice of parallelizing (or not parallelizing) the proposal step to the practitioner.

While parallel MCMC algorithms are increasingly well-suited for developing many-core computational architectures, there are trade-offs that need to be taken into account when choosing how to allocate computational resources. On one end of the spectrum, Gelman and Rubin (1992a, b) demonstrate the usefulness of generating, combining and comparing multiple independent Markov chains that target the same distribution, and one may perform this embarrassingly parallel task by assigning the operations for each individual chain to a separate central processing unit (CPU) core or GPU work-group. In this multi-chain context, simultaneously running multiple parallel MCMC chains could limit resources available for the within-chain parallelization described in the previous paragraph. On the other end of the spectrum, one may find it useful to allocate resources to overcome computational bottlenecks within a single Markov chain that uses a traditional accept-reject step. In big data contexts, Holbrook et al. (2021a, b, 2022a, 2022b) use multi- and many-core processing to accelerate log-likelihood and log-likelihood gradient calculations within single Metropolis-Hastings and Hamiltonian Monte Carlo (Neal, 2011) generated chains. This big data strategy might again limit resources available for parallelization across proposals and target evaluations described above.

In the presence of these trade-offs, it is worthwhile to develop additional parallel computing tools for parallel MCMC so that future scientists may be better able to flexibly assign sub-routines to different computing resources in a way that is tailored to their specific inference task. In particular, we show that quantum parallelization can be a useful tool for parallel MCMC when evaluation of the target density represents the computational bottleneck.

Quantum computers use quantum mechanics to store information and manipulate data. By leveraging the idiosyncracies of quantum physics, these computers are able to deliver speedups over conventional computing for a relatively small number of computational problems. Some of these speedups are very large. Quantum computers can achieve exponential speedups over conventional computers for some computational tasks.. Shor’s quantum algorithm for factoring an integer $N$ (Shor, 1994) is polynomial in $\log N$ compared to a super-polynomial classical optimum. The HHL algorithm for solving sparse $N$-dimensional linear systems (Harrow et al., 2009) is ${\cal O}\left(\log N\right)$ compared to a classical ${\cal O}\left(N\right)$. Other algorithms deliver a still impressive polynomial speedup. For example, the algorithms considered in the following (Section 2.2) achieve quadratic speedups over conventional computing, turning ${\cal O}\left(N\right)$ runtimes into $\mathcal{O}(\sqrt{N})$. Both the magnitude of quantum computing’s successes and the relative rarity of those successes mean that there is a general interest in leveraging quantum algorithms for previously unconsidered computational problems. We will show that—with the help of the Gumbel-max trick (Section B)—established quantum optimization techniques are actually useful for sampling from computationally expensive discrete distributions and apply this insight to achieving quadratic speedups for parallel MCMC.

We provide a short introduction to the most basic elements of quantum computing in Appendix A. The interested reader may look to Nielsen and Chuang (2002) for a canonical introduction to quantum algorithms or Lopatnikova and Tran (2021); Wang and Liu (2022) for surveys geared toward statisticians. In Section 2.1, we review and compare three different quantum search algorithms that enjoy quadratic speedups over conventional algorithms. In Section 2.2 we review the quantum minimization algorithm of Durr and Hoyer (1996) and investigate the use of warm-starting therein. We then combine the Gumbel-max trick (Section B) and recent insights in parallel MCMC with the quantum minimization algorithm to create the quantum parallel MCMC algorithm for both continuously-valued (Section 3.2) and discrete-valued (Section 3.3) targets.

Grover (1996) demonstrates how use a quantum computer to find a single marked element within a finite set of $N$ items with only $\mathcal{O}(\sqrt{N})$ queries, and a result from Bennett et al. (1997) shows that Grover’s algorithm is optimal up to a multiplicative constant. Boyer et al. (1998) extend Grover’s algorithm to multiple marked items or solutions, further extend it to the case when the number of solutions is unknown and provide a rigorous bounds on the algorithms’ error rates. Finally, Durr and Hoyer (1996) use the results of Boyer et al. (1998) to obtain the minimum value within a discrete ordered set. Here, we briefly review these advances and provide illustrative simulations. In Section B, the Gumbel-max trick allows us to extend these results to the problem of sampling from a discrete distribution.

With the rudiments of quantum minimization in hand, we present a quantum parallel MCMC (QPMCMC). The general structure of the algorithm the same as Algorithm 1: it constructs a Markov chain by iterating between generating many proposals and randomly selecting new Markov chain states among these proposals. Unlike classical single-proposal MCMC, parallel MCMC proposes many states and chooses one according to the generalized acceptance probabilities of (2). In general MCMC, evaluation of the target density function $\pi(\cdot)$ is the rate-limiting computational step, becoming particularly onerous in big data scenarios (Holbrook et al., 2021a). While parallel MCMC benefits from improved mixing as a function of MCMC iterations, the increased computational burden of $P$ target evaluations at each step can lead to less favorable comparisons when accounting for elapsed wall-clock time.

Having successfully generated proposals $\bm{\theta}_{p},$ $p=1,\dots,P$ and evaluated the corresponding proposal densities $q(\cdot,\cdot)$, we would like to use quantum parallelism and an appropriate oracle gate to efficiently compute the $\pi(\bm{\theta}_{p})$s but immediately encounter a Catch-22 when we seek to draw a sample $\bm{\theta}^{(s+1)}\sim$ Discrete$(\bm{\pi})$ for $\bm{\pi}$ the vector of probabilities $\pi_{0},\pi_{1},\dots,\pi_{P}$ defined in (2). We can use neither a quantum nor a classical circuit to draw the sample! On the one hand, drawing the sample within a quantum circuit would somehow require that all superposed states have access to all the $\pi(\bm{\theta}_{p})$s at once to perform the required normalization. On the other hand, drawing the sample within the classical circuit would require access to all the $\pi(\bm{\theta}_{p})$s, but measurement in the computational basis can only return one.

In light of this dilemma, we propose to use the Gumbel-max trick (Papandreou and Yuille, 2011) to transform the generalized accept-reject step into a discrete optimization procedure (Section B). From there, we use quantum minimization to efficiently sample from $Discrete(\bm{\pi})$. Crucially, we get away with quantum parallel evaluation of the target $\pi(\cdot)$ over all superposed states because the Gumbel-max trick does not rely on a normalization step: each superposed state requires no knowledge of any target evaluation other than its own.

Once one has effectively parallelized the target evaluations $\pi(\bm{\theta}_{p})$, the new computational bottleneck is the calculations required to obtain the $P+1$ proposal density terms $q(\bm{\theta}_{p},\bm{\Theta}_{-p})$, for $p=0,\dots,P$. Under an independent proposal mechanism with a proposal distribution of the form $q(\bm{\theta}_{0},\bm{\Theta}_{-0})=\prod_{p=1}^{P}q(\bm{\theta}_{0},\bm{\theta}_{p})$, the acceptance probabilities (2) require the ${\cal O}\left(P^{2}\right)$ evaluation of the $P+1\choose 2$ terms $q(\bm{\theta}_{p},\bm{\theta}_{p^{\prime}})$. To avoid such calculations, we follow Holbrook (2021) and use symmetric proposal distributions for which $q(\bm{\theta}_{p},\bm{\Theta}_{-p})=q(\bm{\theta}_{p^{\prime}},\bm{\Theta}_{-p^{\prime}})$ for $p,p^{\prime}=0,\dots,P$.

We have shown that parallel MCMC enjoys quadratic speedups by combining quantum minimization with the Gumbel-max trick. Within a QPMCMC iteration, the current state represents a warm-start for the sampling-turned-optimization problem, leading to increased efficiencies for the quantum minimization algorithm. Moreover, combining this approach with the Tjelmeland correction (Holbrook, 2021) results in a total complexity reduction from $\mathcal{O}(P^{2})$ to $\mathcal{O}(\sqrt{P})$. Preliminary evidence suggests that our strategy may find use when target distributions exhibit extreme multimodality. While the algorithm still must construct long Markov chains to reach equilibrium, generating each individual Markov chain state requires significantly fewer target evaluations.

There are major technical barriers to the practical implementation of QPMCMC. The framework, like other quantum machine learning (QML) schemes, requires on-loading and off-loading classical data to and from a quantum machine. In the seminal review of QML, Biamonte et al. (2017) discuss what they call the ‘input problem’. Cortese and Braje (2018) present a general purpose quantum circuit for loading $B$ classical bits into a quantum data structure comprising $\log_{2}(B)$ qubits with circuit depth ${\cal O}\left(\log(B)\right)$. For continuous targets, QPMCMC requires reading ${\cal O}\left(P\right)$ real vectors $\bm{\theta}_{p}\in\mathbb{R}^{D}$ onto the quantum machine at each MCMC iteration. If $b$ is the number of bits used to represent a single real value, then reading $B={\cal O}\left(PDb\right)$ classical bits into the quantum machine requires time ${\cal O}\left(\log(PDb)\right)$. This is true whether one opts for fixed-point (Jordan, 2005) or floating-point (Haener et al., 2018) representations for real values. The outlook can be better for discrete distributions. For example, the QPMCMC scheme presented in Section 4 only requires loading the entire configuration state once prior to sampling. A $D$-spin Ising model requires $D$ classical bits to encode, and these bits load onto a quantum machine in time ${\cal O}\left(\log(D)\right)$. Once one has encoded the entire system state, each QPMCMC iteration only requires loading the addresses of the ${\cal O}\left(P\right)$ proposal states. If one uses $b$ bits to encode each address, then the total time required to load data onto the quantum machine is ${\cal O}\left(\log(Pb)\right)$ for each QPMCMC iteration. On the other hand, the speedup for discrete targets assumes the ability to hold an entire configuration in QRAM. Conveniently, the ‘output problem’ is less of an issue for QPMCMC, as only a single integer $\hat{p}\in\{0,\dots,P\}$ need be extracted within Algorithm 4.

This work leads to three interesting questions. First, what is the status of inference obtained by QPMCMC? The QPMCMC selection step relies on quantum minimization, an algorithm that only achieves success with probability $1-\epsilon$. While empirical studies suggest that this error induces little bias, it would be helpful to use this $\epsilon$ to bound the distance between the target distribution and the distribution generated by QPMCMC. Such theoretical efforts would need to extend recent formalizations of the multiproposal based parallel MCMC paradigm (Glatt-Holtz et al., 2022) to account for biased MCMC kernels.

Second, can QPMCMC be combined with established quantum algorithms that make use of ‘quantum walks’ on graphs to sample from discrete target distributions? Szegedy (2004) presents a quantum analog to classical ergodic reversible Markov chains and shows that such quantum walks sometimes provide quadratic speedups over classical Markov chains. Szegedy (2004) also points out that Grover’s search, a key component within QPMCMC, may be interpreted as performing just such a quantum walk on a complete graph. Wocjan and Abeyesinghe (2008) accelerate the quantum walk by using ancillary Markov chains to improve mixing and apply their method to simulated annealing. Given a quantum algorithm $\mathbb{A}$ for producing a discrete random sample with variance $\sigma^{2}$, Montanaro (2015) develops a quantum algorithm for estimating the mean of algorithm $\mathbb{A}$’s output with error $\epsilon$ by running algorithm $\mathbb{A}$ a mere $\widetilde{\mathcal{O}}(\sigma/\epsilon)$ times, where $\widetilde{\mathcal{O}}$ hides polylogarithmic terms. Importantly, this quadratic quantum speedup over classical Monte Carlo applies for quantum algorithms $\mathbb{A}$ that only feature a single measurement such as certain simple quantum walk algorithms. Unfortunately, this assumption fails for quantizations of Metropolis-Hastings, in general, and QPMCMC, in particular. More promising for QPMCMC, Lemieux et al. (2020) develop a quantum circuit that applies Metropolis-Hastings to the Ising model without the need for oracle calls. An interesting question is whether similar non-oracular quantum circuits exist for the basic parallel MCMC backbone to QPMCMC. In general, however, comparison between QPMCMC and other quantum Monte Carlo techniques is challenging because the foregoing literature (a) largely focuses on MCMC as a tool for discrete optimization, with algorithms that only ever return a single Monte Carlo sample or function thereof, and (b) restricts itself to a handfull of simple, stylized and discrete target distributions. On the other hand, QPMCMC is a general inferential framework for sampling from general discrete and continuous distributions alike.

Third, there are other models and algorithms in statistics, machine learning and statistical mechanics that require sampling from potentially costly discrete distributions. Can one adapt our quantum Gumbel-max trick to these? Approximate Bayesian computation (Csilléry et al., 2010) extends probabilistic inference to contexts within which one only has access to a complex, perhaps computationally intensive, forward model. Having sampled model parameters from the prior, one accepts or rejects these parameter values based on the discrepancy between the observed and simulated data. If one succeeds in embedding forward model dynamics within a quantum circuit, then one may plausibly select from many parameter values using our quantum Gumbel-max trick. The trick may also find use within sequential Monte Carlo (Doucet et al., 2001). For example, Berzuini and Gilks (2001) present a sequential importance resampling algorithm that uses MCMC-type moves to encourage particle diversity and avoid the need for bootstrap resampling. Multiproposals accelerated by the quantum Gumbel-max trick could add speed and robustness to such MCMC-resampling.

Large-scale, practical quantum computing is still a long way off, but quantum algorithms are ripe for mainstream computational statistics.

This work was supported by grants NIH K25 AI153816, NSF DMS 2152774 and NSF DMS 2236854.