A Quantum Parallel Markov Chain Monte Carlo

Abstract

We propose a novel hybrid quantum computing strategy for parallel MCMC algorithms that generate multiple proposals at each step. This strategy makes the rate-limiting step within parallel MCMC amenable to quantum parallelization by using the Gumbel-max trick to turn the generalized accept-reject step into a discrete optimization problem. When combined with new insights from the parallel MCMC literature, such an approach allows us to embed target density evaluations within a well-known extension of Grover’s quantum search algorithm. Letting P denote the number of proposals in a single MCMC iteration, the combined strategy reduces the number of target evaluations required from $\mathcal{O}(P)$ to $\mathcal{O}\left(\sqrt{P}\right)$. In the following, we review the rudiments of quantum computing, quantum search and the Gumbel-max trick in order to elucidate their combination for as wide a readership as possible.

Keywords:

• Algorithms
• Markov chain Monte Carlo (MCMC)
• Parallel computing
• Statistical computing

1 Introduction

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 (\mathrm{d}\mathit{\theta })$ defined on ${\mathbb{R}}^D$ that admits a probability density $\pi (\mathit{\theta })$ with respect to the Lebesgue measure, that is, $\pi (\mathrm{d}\mathit{\theta })=:\pi (\mathit{\theta })\mathrm{d}\mathit{\theta }$. To generate samples from the target distribution π, we craft a kernel $P({\mathit{\theta }}_0,\mathrm{d}\mathit{\theta })$ that satisfies(1) $$\pi (A)=\int \pi (\mathrm{d}{\mathit{\theta }}_0)P({\mathit{\theta }}_0,A)$$(1) for all $A\subset {\mathbb{R}}^D$. Letting ${\mathit{\theta }}_0$ denote the current state of the Markov chain, Tjelmeland (2004) proposes sampling from such a kernel $P({\mathit{\theta }}_0,·)$ by drawing P proposals ${\mathbf{\Theta }}_{-0}=({\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_P)$ from a joint distribution $Q({\mathit{\theta }}_0,\mathrm{d}{\mathbf{\Theta }}_{-0})=:q({\mathit{\theta }}_0,{\mathbf{\Theta }}_{-0})\mathrm{d}{\mathbf{\Theta }}_{-0}$ and selecting the next Markov chain state from among the current and proposed states with probabilities(2) $${\pi }_p=\frac{\pi ({\mathit{\theta }}_p)q({\mathit{\theta }}_p,{\mathbf{\Theta }}_{-p})}{\sum _{p\prime =0}^P\pi ({\mathit{\theta }}_{p\prime })q({\mathit{\theta }}_{p\prime },{\mathbf{\Theta }}_{-p\prime })} ,\hspace{1em}p=0,\dots ,P .$$(2)

Here, if $\mathbf{\Theta }=({\mathit{\theta }}_0,\dots ,{\mathit{\theta }}_P)$ is the $D\times (P+1)$ matrix having the current state and all P proposals as columns, then ${\mathbf{\Theta }}_{-p}$ is the D × P matrix obtained by removing the pth column. Tjelmeland (2004) shows that the kernel $P({\mathit{\theta }}_0,·)$ 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 theoretical 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 ({\mathit{\theta }}_p)$ is linear in the number of proposals. This presents a significant burden when $\pi (·)$ is computationally expensive, for example, 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({\mathit{\theta }}_0,{\mathbf{\Theta }}_{-0})={\prod }_{p=1}^{P}q({\mathit{\theta }}_0,{\mathit{\theta }}_p)$, the acceptance probabilities (2) require the $\mathcal{O}(P^2)$ evaluation of the $(\begin{array}{c}P+1\\ 2\end{array})$ terms $q({\mathit{\theta }}_p,{\mathit{\theta }}_{p\prime })$, but Holbrook et al. (2021a, 2021b) 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, 1992b) 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, 2021b) and Holbrook, Ji, and Suchard (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.

Algorithm 1:

Parallel (multiproposal) MCMC (Tjelmeland 2004)

Data: Initial Markov chain state ${\mathit{\theta }}^{(0)}$; total length of Markov chain S; total number of proposals per iteration P; routine for evaluating target density $\pi (·)$; routines for drawing random samples from the proposal distribution $Q({\mathit{\theta }}^{(s)},·)$ and from a P + 1 discrete distribution $\text{Discrete}(\mathit{\pi })$ given some probability vector $\mathit{\pi }=({\pi }_0,\dots ,{\pi }_P)$.

Result: A Markov chain ${\mathit{\theta }}^{(1)},\dots ,{\mathit{\theta }}^{(S)}$.for $s\in \{1,\dots ,S\}$ do

${\mathit{\theta }}_0\leftarrow {\mathit{\theta }}^{(s-1)}$; ${\mathbf{\Theta }}_{-0}\leftarrow Q({\mathit{\theta }}_0,·)$; for $p\in \{0,\dots ,P\}$ do $${\pi }_p\leftarrow \pi ({\mathit{\theta }}_p) q({\mathit{\theta }}_p,{\mathbf{\Theta }}_{-p});$$end $$\widehat{p}\leftarrow \text{Discrete}(\mathit{\pi }/{\mathit{\pi }}^{T}1);{\mathit{\theta }}^{(s)}\leftarrow {\mathit{\theta }}_{\widehat{p}};$$

end return ${\mathit{\theta }}^{(1)},\dots ,{\mathit{\theta }}^{(S)}$.

In the presence of these tradeoffs, 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 $\hspace{0.17em}\log \hspace{0.17em}N$ compared to a super-polynomial classical optimum. The HHL algorithm for solving sparse N-dimensional linear systems (Harrow, Hassidim, and Lloyd 2009) is $\mathcal{O}(\hspace{0.17em}\log \hspace{0.17em}N)$ compared to a classical $\mathcal{O}(N)$. 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 $\mathcal{O}(N)$ runtimes into $\mathcal{O}\left(\sqrt{N}\right)$. 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 (2010) 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.

2 Quantum Search and Quantum Minimization

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}\left(\sqrt{N}\right)$ 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.

2.1 Quantum Search

Algorithm 2:

Quantum search algorithm (Grover 1996)

Data: An oracle gate ${\mathbf{U}}_f$ taking $|x〉|y〉$ to $|x〉|y\oplus f(x)〉$ for a function $f(x):\{0,\dots ,N-1\}\to \{0,1\}$ that satisfies $f(x_0)=1$ for a single x₀; $n+1={\hspace{0.17em}\log \hspace{0.17em}}_2(N)+1$ quantum states initialized to $\left|0{〉}^{\otimes n}\right|1〉$; an integer $R=\lceil \pi \sqrt{N}/4\rceil$ denoting the number of Grover iterations.

Result: An n-bit binary string x₀ satisfying $f(x_0)=1$ with error less than $1/N$.

$\left|0{〉}^{\otimes n}\right|1〉\to {\mathbf{H}}^{\otimes n+1}\left|0{〉}^{\otimes n}\right|1〉=|h〉|-〉$;/* $|h〉=\frac{1}{\sqrt{N}}{\sum }_{x=0}^{N-1}|x〉$. */

$|h〉|-〉\to {\mathbf{G}}^R|h〉|-〉\approx |x_0〉|-〉$;/* $\mathbf{G}=(2|h〉〈h|-\mathbf{I})(\mathbf{I}-2|x_0〉〈x_0|)$ */

$|x_0〉\to x_0$; return x₀ .

Given a set of N items and a function $f:\{0,1,\dots ,N-1\}\to \{0,1\}$ that evaluates to 1 for a single element, Grover (1996) develops an algorithm that uses quantum parallelism to score quadratic speedups over its classical counterparts. After only $\mathcal{O}\left(\sqrt{N}\right)$ evaluations of $f(·)$, Grover’s algorithm returns the $x\in \{0,1,\dots ,N-1\}$ satisfying f(x) = 1 with high probability. Compare this to $\mathcal{O}(N)$ requirement for the randomized classical algorithm that must evaluate $f(·)$ over at least $N/2$ items to obtain the same probability of detecting x. The algorithm takes the state $\left|0{〉}^{\otimes N}\right|1〉$ as input and applies Hadamard gates to each of the individual N + 1 input qubits. The resulting state is$$|h〉|-〉=\left(\frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}|x〉\right)\frac{1}{\sqrt{2}}\left(|0〉-|1〉\right)=\frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}|x〉|-〉 .$$

Next, we apply the oracle gate ${\mathbf{U}}_f:|x〉|y〉\to |x〉|y\oplus f(x)〉$ and note that$$\begin{array}{c}{\mathbf{U}}_f|x〉|-〉={\mathbf{U}}_f|x〉\frac{1}{\sqrt{2}}\left(|0〉-|1〉\right)\\ =\frac{1}{\sqrt{2}}(|x〉|0\oplus f(x)〉-|x〉|1\oplus f(x)〉)\\ =-1^{f(x)}|x〉|-〉 .\end{array}$$

Thus, ${\mathbf{U}}_f$ flips the sign for the state $|x_0〉$ for which $f(x_0)=1$ but leaves the other states unchanged. If we suppress the ancillary qubit $|-〉$, then ${\mathbf{U}}_f$ is equivalent to the gate ${\mathbf{U}}_{x_0}$ defined as ${\mathbf{U}}_{x_0}|x〉=-1^{{\delta }_{x_0}}|x〉$. We may succinctly write this gate as the Householder matrix that reflects vectors about the unique hyperplane through the origin that has $|x_0〉$ for a normal vector:$${\mathbf{U}}_{x_0}=\mathbf{I}-2|x_0〉〈x_0| .$$

The action of this gate on the state $|h〉$ takes the form$$\frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}|x〉\to \frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}-1^{{\delta }_{x_0}(x)}|x〉 .$$

Next, the algorithm reflects the current state about the hyperplane that has $|h〉$ as a normal vector and negates the resulting state:$$\begin{array}{c}(2|h〉〈h|-\mathbf{I})\left(\frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}-1^{{\delta }_{x_0}(x)}|x〉\right)\\ \hspace{1em}=\frac{1}{\sqrt{N}}\sum _x\left({(-1)}^{1-{\delta }_{x_0}(x)}+2\frac{(N-2)}{N}\right)|x〉\\ \hspace{1em}=\frac{(3N-4)}{N^{3/2}}|x_0〉+\sum _{x\ne x_0}\frac{(N-4)}{N^{3/2}}|x〉 .\end{array}$$

The scientist who measures the state at this moment would obtain the desired state $|x_0〉$ with a slightly higher probability of ${(3N-4)}^2/N^3$ than the individual probabilities of ${(N-4)}^2/N^3$ for the other states. Each additional application of the Grover iteration (3) $$\mathbf{G}:=-{\mathbf{U}}_h{\mathbf{U}}_{x_0}:=(2|h〉〈h|-\mathbf{I})(\mathbf{I}-2|x_0〉〈x_0|)$$(3) increases the probability of obtaining $|x_0〉$ at the time of measurement in the computational basis and $R=\lceil \pi \sqrt{N}/4\rceil$ iterations guarantees a probability of success that is greater than $1-1/N$.

In general, we may use Grover’s search algorithm to find a solution when the number of solutions M is greater than 1. While the algorithmic details change little, the number of required Grover iterations(4) $$R=\lceil \frac{\pi }{4}\sqrt{\frac{N}{M}}\rceil$$(4) and the probability of success after those R iterations do change (Boyer et al. 1998). When M is much smaller than N, the success rate is greater than $1-M/N$, and even for large M the success rate is more than 1/2. Lower-bounds are useful for establishing mathematical guarantees, but it is also helpful to understand the quality of algorithmic performance as a function of M and R. Figure 1 shows success curves as a function of the number of Grover iterations (or oracle calls) applied. The search field contains $2^{14}\approx 16$ k elements, and the number of solutions M varies. Each curve represents an individual search task. The algorithm requires more iterations for smaller numbers and fewer iterations for larger numbers of solutions. The upper bound on error M/N is only an upper bound: the final probability of success for, for example, M = 256 is 0.997 compared to the bound of 0.984.

Fig. 1 Curves depict the probability that Grover’s quantum search algorithm returns a solution as a function of the number of oracle evaluations. For the ranges depicted, the number of iterations required and the number of solutions inversely relate, but increasing the number of iterations past R (Equation (4)(4) $$R=\lceil \frac{\pi }{4}\sqrt{\frac{N}{M}}\rceil$$(4) ) can backfire and lead to extremely small probabilities of success.

Algorithm 3:

Quantum exponential searching algorithm (Boyer et al. 1998)

Data: An oracle gate ${\mathbf{U}}_f$ taking $|x〉|y〉$ to $|x〉|y\oplus f(x)〉$ for a function $f(x):\{0,\dots ,N-1\}\to \{0,1\}$ with unknown number of solutions; $n={\hspace{0.17em}\log \hspace{0.17em}}_2(N)$.

Result: If a solution exists, an n-bit binary string x₀ satisfying $f(x_0)=1$; if no solution exists, the algorithm runs forever.

$m\leftarrow 1$; $\gamma \leftarrow 6/5$; success $\leftarrow$ FALSE; while success $\ne$TRUE do

$j\leftarrow \mathit{Uniform}\{0,\dots ,m-1\}$; $\left|0{〉}^{\otimes n}\right|1〉\to {\mathbf{H}}^{\otimes n+1}\left|0{〉}^{\otimes n}\right|1〉=|h〉|-〉$; $|h〉|-〉\to {\mathbf{G}}^j|h〉|-〉=|x〉|-〉$;/* j Grover iterations from Algorithm 2 . */

$|x〉\to x$;/* Measure and check. */if f(x) = 1 then

$x_0\leftarrow x$; success $\leftarrow$ TRUE;else

$m\leftarrow \min \left(\gamma m,\sqrt{N}\right)$;/* Increase m in case of failure. */end end return x₀ .

While Grover’s algorithm delivers impressive speedups over classical search algorithms, it has a major weakness. Figure 1 hides the fact that the probability of success for Grover’s algorithm is not monatonic in the number of iterations. Running the algorithm for more than R iterations can backfire. For example, running the algorithm for $\sqrt{2N}$ iterations when M = 1 results in a probability of success less than 0.095 (Boyer et al. 1998). The non-monotinicity of Grover’s algorithm becomes particularly problematic when we do not know the number of solutions M. Taking for example $N=2^{20}$, Boyer et al. (1998) point out that 804 iterations provide an extremely high probability of success when M = 1 but a one-in-a-million chance of success when M = 4. To solve this problem and develop an effective search algorithm when M is unknown, those authors adopt the strategy of focusing on the expected number of iterations before success. In particular, they propose the quantum exponential search algorithm (Algorithm 3). When a solution exists, the algorithm returns a solution with expected total number of Grover iterations bounded above by $\frac{9}{2}\sqrt{N/M}$. Still better, this upper bound reduces to $\frac{9}{4}\sqrt{N/M}$ for the special case $M\ll N$, and simulations presented in Figure 2 show that even this bound is large. Such results come in handy when deciding whether to halt the algorithm’s progress if one believes it possible that no solutions exist. Indeed, this turns out to be useful in the context of quantum minimization (Section 2.2).

Fig. 2 [Left] Total number of oracle evaluations required by quantum exponential searching algorithm (Algorithm 3) for different numbers of solutions M and search set sizes N from 500 independent simulations each. Horizontal lines at $\frac{9}{4}\sqrt{N}$ represent upper bounds on expected total number of evaluations to obtain a solution for the M = 1 problem. [Right] Probabilities of success for the fixed-point quantum search algorithm (Yoder, Low, and Chuang 2014) for different proportions of solutions $\lambda =M/N$ and selecting different lower bounds w on M/N with error tolerance ${\delta }^2$.

Algorithm 3 is not the only quantum search algorithm that is useful when the number of solutions M is unknown. Yoder, Low, and Chuang (2014) use generalized Grover iterations that extend (3) to include general phases $({\alpha }_j,{\beta }_j)$:(5) $$\begin{array}{c}\mathbf{G}({\alpha }_j,{\beta }_j):=-{\mathbf{U}}_h({\alpha }_j){\mathbf{U}}_{x_0}({\beta }_j)\\ :=((1-e^{-i{\alpha }_j})|h〉〈h|-\mathbf{I})(\mathbf{I}-(1-e^{i{\beta }_j})|x_0〉〈x_0|).\end{array}$$(5)

This form reduces to the original Grover iteration (3) when $\alpha =\pm \pi$ and $\beta =\pm \pi$. To select general phases, the method first fixes an error tolerance ${\delta }^2>0$ and a lower bound $w\le \lambda =M/N$ and then selects an odd upper bound L on the total number of oracle queries (L – 1) such that$$L\ge \frac{\hspace{0.17em}\log \hspace{0.17em}(2/\delta )}{\sqrt{w}} .$$

Finally, letting $l=(L-1)/2$, one obtains phases $({\alpha }_j,{\beta }_j)$ for $j=1,\dots ,l$ that satisfy$${\alpha }_j=-{\beta }_{l-j+1}=2{\cot }^{-1}\left(\hspace{0.17em}\tan \hspace{0.17em}(2\pi j/L)\sqrt{1-{\gamma }^2}\right) ,$$for ${\gamma }^{-1}=\hspace{0.17em}\cos \hspace{0.17em}(\frac{1}{L}{\hspace{0.17em}\cos \hspace{0.17em}}^{-1}(\frac{1}{\delta }))$ is the reciprocal of the inverse Lth Chebyshev polynomial of the first kind. Such phases mark the only algorithmic difference between the fixed-point quantum search and the original Grover search (Algorithm 2). The upshot is a search procedure that obtains a guaranteed probability of success $1-{\delta }^2$ for any M such that there exists an $M_0<M$ that also obtains the same success threshold (Figure 2). On the one hand, this algorithm avoids the exponential quantum search algorithm’s need to perform multiple measurements of the quantum state. On the other hand, the exponential quantum search algorithm requires significantly fewer oracle evaluations when M is small (Figure 3), and it turns out that this is precisely the scenario that interests us.

Fig. 3 Total number of oracle evaluations required by fixed-point quantum search (FPQS) (Yoder, Low, and Chuang 2014) minus total required for quantum exponential searching algorithm (QESA) (Algorithm 3) for different numbers of solutions M and search set sizes N from 500 independent simulations each. FPQS underperforms partially due to a miniscule lower bound $w=1/N$ on $\lambda =M/N$, a tuning decision motivated by the potential application within quantum minimization with warm-starting (Section 2.2).

2.2 Quantum Minimization

Given a function $f(·)$ that maps the discrete set $\{0,\dots ,N-1\}$ to the integers, we wish to find the minimizer$$x_0=\text{arg} {\min }_{x\in \{0,\dots ,N-1\}}f(x) .$$

Durr and Hoyer (1996) propose a quantum algorithm for finding x₀ that iterates between updating a running minimum $F\in f(\{0,\dots ,N-1\})$ and applying the quantum exponential search algorithm (Algorithm 3) to find an element x such that f(x) < F. Having run these iterations a set number of times, the algorithm returns the minimizer x₀ with high probability. To this end, Durr and Hoyer (1996) show that their algorithm takes an expected total time less than or equal to $m_0:=\frac{45}{4}\sqrt{N}+\frac{7}{10}{\hspace{0.17em}\log \hspace{0.17em}}_2(N)$ to find the minimizer, where marking items with values less than the threshold value (Algorithm 4) requires ${\hspace{0.17em}\log \hspace{0.17em}}_2(N)$ time steps and each Grover iteration within the quantum exponential search algorithm requires one time step. From there, Markov’s inequality says$$\begin{array}{c}\text{Pr}\left(\mathrm{ }\text{total}\mathrm{ }\text{time}\mathrm{ }\text{to}\mathrm{ }\text{success}\ge \frac{m_0}{ϵ}\right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\le \frac{\mathrm{ }\mathrm{E}\left(\mathrm{ }\text{total}\mathrm{ }\text{time}\mathrm{ }\text{to}\mathrm{ }\text{success}\right)}{m_0/ϵ}\le ϵ,\end{array}$$or that we must scale the minimization procedure’s time steps by a factor of $1/ϵ$ to reach a failure rate less than ϵ. Due to the iterative nature of the algorithm, one might suppose that it is beneficial to start at a value x₀ for which $f(x_0)$ is lower relative to other values. This turns out to be the case in theory and in practice, and the benefit of warm-starting is particularly useful in the context of parallel MCMC.

Proposition 1

(Warm-starting). Suppose that the quantum minimization algorithm begins with a threshold F₀ such that $f(x)<F_0$ for only $K-1>0$ items. Then the expected total number of time steps to find the minimizer is bounded above by$$m_0^K=\left(\frac{5}{4}-\frac{1}{\sqrt{K-1}}\right)9\sqrt{N}+\frac{7}{10}{\hspace{0.17em}\log \hspace{0.17em}}_2(K){\hspace{0.17em}\log \hspace{0.17em}}_2(N) ,$$and so the following rule relates the warm-started upper bound to the generic upper bound:$$m_0^K=m_0-9\sqrt{\frac{N}{K-1}}+\frac{7}{10}{\hspace{0.17em}\log \hspace{0.17em}}_2\left(\frac{K}{N}\right){\hspace{0.17em}\log \hspace{0.17em}}_2(N) .$$

Proof.

The proof follows the exact same form as Lemma 2 of Durr and Hoyer (1996). It relies on a theoretical algorithm called the infinite algorithm that runs until the minimum is found. In this case, Lemma 1 of that paper says that the probability that the rth lowest value is ever selected when searching among K items is $p(K,r)=1/r$ for $r\le K$ and 0 otherwise. For a warm-start at element K, the expected total time spent in the exponential search algorithm is$$\begin{array}{c}\sum _{r=2}^{N}p(K,r)\frac{9}{2}\sqrt{\frac{N}{r-1}}=\sum _{r=2}^{K}p(K,r)\frac{9}{2}\sqrt{\frac{N}{r-1}}\\ =\frac{9}{2}\sqrt{N}\sum _{r=1}^{K-1}\frac{1}{r+1}\frac{1}{\sqrt{r}}\\ \le \frac{9}{2}\sqrt{N}\left(\frac{1}{2}+\sum _{r=2}^{K-1}{r}^{-3/2}\right)\\ \le \frac{9}{2}\sqrt{N}\left(\frac{1}{2}+{\int }_{r=1}^{K-1}{r}^{-3/2}\mathrm{d}r\right)\\ =\left(\frac{5}{4}-\frac{1}{\sqrt{K-1}}\right)9\sqrt{N}.\end{array}$$

An upper bound for the expected total number of time steps preparing the initial state and marking items $f(x)<f(x_0)$ follows in a similar manner. □

Proposition 1 shows that, for example, if Algorithm 4 begins at the item with second lowest value, then the expected total time to success is bounded above by $m_0^2=\frac{9}{4}\sqrt{N}+\frac{7}{10}{\hspace{0.17em}\log \hspace{0.17em}}_2(N)$, reducing the generic upper bound by $9\sqrt{N}-\frac{7}{10}{\hspace{0.17em}\log \hspace{0.17em}}_2\left(\frac{2}{N}\right){\hspace{0.17em}\log \hspace{0.17em}}_2(N)$ time steps. When N equals 1000, say, the expected total time steps is $m_0^2<78.2$. Keeping N = 1000 but letting the algorithm begin at the third lowest value (K = 3), the expected total time steps before success is $m_0^3<165.6$. Raising N to 10,000, the two numbers increase to $m_0^2<234.4$ and $m_0^3<503.4$.

Algorithm 4:

Quantum minimum searching algorithm (Durr and Hoyer 1996)

Data: A quantum sub-routine capable of evaluating a function $f(·)$ over $\{0,\dots ,N-1\}$ with with unique integer values; a maximum error tolerance $ϵ\in (0,1)$; expected total time to success $m_0=\frac{45}{4}\sqrt{N}+\frac{7}{10}{\hspace{0.17em}\log \hspace{0.17em}}_2(N)$.

Result: A ${\hspace{0.17em}\log \hspace{0.17em}}_2(N)$-bit binary string x₀ satisfying $f(x_0)=\min f$ with probability greater than $1-ϵ$. $s\leftarrow 0$; $x_0\leftarrow \text{Uniform}\{0,\dots ,N-1\}$; while $s<m_0/ϵ$ do

Prepare initial state $\frac{1}{\sqrt{N}}\sum _x|x〉|x_0〉$; Mark all items x satisfying $f(x)<f(x_0)$; $s\leftarrow s+{\hspace{0.17em}\log \hspace{0.17em}}_2(N)$; Apply quantum exponential searching algorithm (Algorithm 3);/* I time steps */

$s\leftarrow s+I$; Obtain $x\prime$ by measuring first register; if $f(x\prime )<f(x_0)$ then $$x_0\leftarrow x\prime$$end end return x₀ .

In practice, the number of time steps before finding the minimum is surprisingly small. Figure 4 corroborates the intuition that it is beneficial to start at a lower ranked element. To generate these results, we use an early stopping rule for the quantum minimization algorithm’s exponential searching sub-routine, halting it after $\frac{9}{4}\sqrt{N}$ iterations. Even with this stopping rule in place, the error rate of the quantum minimization algorithm is less than 1%. We can increase the early stopping threshold to obtain even lower error rates, but this strategy would seem to be unnecessary in the context of parallel MCMC. The benefits of warm-starting are useful in the same context, when the current state ${\mathit{\theta }}^{(s)}$ usually inhabits the high-density region of the target distribution but the majority of proposals do not.

Fig. 4 Total number of oracle evaluations required by quantum minimization algorithm (Durr and Hoyer 1996) for different starting ranks and search set sizes N from 500 independent simulations each. Here, less than 1% of the 12,500 instances fail to recover the true minimum on account of early stopping.

3 Quantum Parallel MCMC

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 (·)$ 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 ${\mathit{\theta }}_p,$ $p=1,\dots ,P$ and evaluated the corresponding proposal densities $q(·,·)$, we would like to use quantum parallelism and an appropriate oracle gate to efficiently compute the $\pi ({\mathit{\theta }}_p)$ s but immediately encounter a Catch-22 when we seek to draw a sample ${\mathit{\theta }}^{(s+1)}\sim$ Discrete $(\mathit{\pi })$ for $\mathit{\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 ({\mathit{\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 ({\mathit{\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 $\text{Discrete}(\mathit{\pi })$. Crucially, we get away with quantum parallel evaluation of the target $\pi (·)$ 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 ({\mathit{\theta }}_p)$, the new computational bottleneck is the calculations required to obtain the P + 1 proposal density terms $q({\mathit{\theta }}_p,{\mathbf{\Theta }}_{-p})$, for $p=0,\dots ,P$. Under an independent proposal mechanism with a proposal distribution of the form $q({\mathit{\theta }}_0,{\mathbf{\Theta }}_{-0})={\prod }_{p=1}^{P}q({\mathit{\theta }}_0,{\mathit{\theta }}_p)$, the acceptance probabilities (2) require the $\mathcal{O}(P^2)$ evaluation of the $(\begin{array}{c}P+1\\ 2\end{array})$ terms $q({\mathit{\theta }}_p,{\mathit{\theta }}_{p\prime })$. To avoid such calculations, we follow Holbrook (2021) and use symmetric proposal distributions for which $q({\mathit{\theta }}_p,{\mathbf{\Theta }}_{-p})=q({\mathit{\theta }}_{p\prime },{\mathbf{\Theta }}_{-p\prime })$ for $p,p\prime =0,\dots ,P$.

Algorithm 5:

The Gumbel-max trick

Data: A vector of unnormalized log-probabilities ${\mathit{\lambda }}^{*}=\hspace{0.17em}\log \hspace{0.17em}\mathit{\pi }+\hspace{0.17em}\log \hspace{0.17em}c$, for $\mathit{\pi }$ a discrete probability vector with P + 1 elements.

Result: A single sample $\widehat{p}\sim$ Discrete $(\mathit{\pi })$ satisfying $\widehat{p}\in \{0,1,\dots ,P\}$.for $p\in \{0,1,\dots ,P\}$ do $$z_p\leftarrow \text{Gumbel}(0,1);{\alpha }_p^{*}\leftarrow {\lambda }_p^{*}+z_p;$$ end

$\widehat{p}\leftarrow \text{arg} {\max }_{p=0,\dots ,P}{\alpha }_p^{*}$; return $\widehat{p}$.

3.1 Simplified Acceptance Probabilities

Evaluation of the P + 1 proposal density terms $q({\mathit{\theta }}_p,{\mathbf{\Theta }}_{-p})$ can be computationally burdensome as P grows large. Holbrook (2021) proves that two specific proposal mechanisms enforce equality across all P + 1 terms and thus enable the simplified acceptance probabilities(6) $${\pi }_p=\frac{\pi ({\mathit{\theta }}_p)}{\sum _{p\prime =0}^P\pi ({\mathit{\theta }}_{p\prime })} ,\hspace{1em}p=0,\dots ,P .$$(6)

In Section 3.2, we opt for one of these proposal mechanisms, the centered Gaussian proposal of Tjelmeland (2004). This strategy first generates a Gaussian random variable $\overline{\mathit{\theta }}$ centered at the current state ${\mathit{\theta }}_0$ and then generates P proposals centered at $\overline{\mathit{\theta }}$:(7) $${\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_P\stackrel{\perp }{\sim }{\text{Normal}}_D(\overline{\mathit{\theta }},\mathbf{\Sigma }) ,\hspace{1em}\overline{\mathit{\theta }}\sim {\text{Normal}}_D({\mathit{\theta }}_0,\mathbf{\Sigma }) .$$(7)

Holbrook (2021) shows that this strategy leads to the higher-order proposal symmetry(8) $$q({\mathit{\theta }}_0,{\mathbf{\Theta }}_{-0})=q({\mathit{\theta }}_1,{\mathbf{\Theta }}_{-1})=\dots =q({\mathit{\theta }}_P,{\mathbf{\Theta }}_{-P})$$(8) and that the simplified acceptance probabilities (6) maintain detailed balance. In fact, other location scale families also suffice in the manner, but here we focus on Gaussian proposals without loss of generality. Finally, the simplicial sampler (Holbrook 2021) accomplishes the same simplified acceptance probabilities but incurs an $\mathcal{O}(D^3)$ computational cost.

One may also use a more limited strategy to enforce (8). Glatt-Holtz et al. (2022) show that the sampler using independence proposal $q({\mathit{\theta }}_p,{\mathbf{\Theta }}_{-p})=\prod _{p\prime \ne p}q({\mathit{\theta }}_{p\prime })$, where $q(·)$ is not a function of ${\mathit{\theta }}_p$, results in acceptance probabilities(9) $${\pi }_p=\frac{\pi ({\mathit{\theta }}_p)/q({\mathit{\theta }}_p)}{\sum _{p\prime =0}^P\pi ({\mathit{\theta }}_{p\prime })/q({\mathit{\theta }}_{p\prime })} ,\hspace{1em}p=0,\dots ,P .$$(9)

When $\pi (·)$ and $q(·)$ take continuous values on a topologically compact domain or discrete values on a finite domain, one may specify $q(·)$ to be uniform and let (9) simplify to (6). This strategy proves useful in Section 3.3, in which we consider discrete-valued targets over Ising-type lattice models.

3.2 Continuously Valued Targets

Algorithm 6:

Quantum parallel MCMC for a continuously-valued target

Data: Initial Markov chain state ${\mathit{\theta }}^{(0)}$; total length of Markov chain S; total number of proposals per iteration P; routine for evaluating target density $\pi (·)$; routines for drawing random samples from a D-dimensional multivariate Gaussian ${\text{Normal}}_D(\mathit{\mu },\mathbf{\Sigma })$ and the standard Gumbel distribution $\text{Gumbel}(0,1)$.

Result: A Markov chain ${\mathit{\theta }}^{(1)},\dots ,{\mathit{\theta }}^{(S)}$.for $s\in \{1,\dots ,S\}$ do ${\mathit{\theta }}_0\leftarrow {\mathit{\theta }}^{(s-1)}$; $z_0\leftarrow \text{Gumbel}(0,1)$; $\overline{\mathit{\theta }}\leftarrow {\text{Normal}}_D({\mathit{\theta }}_0,\mathbf{\Sigma })$;/* Proposal (7) */

for $p\in \{1,\dots ,P\}$ do ${\mathit{\theta }}_p\leftarrow {\text{Normal}}_D(\overline{\mathit{\theta }},\mathbf{\Sigma })$; $z_p\leftarrow \text{Gumbel}(0,1)$; $|{\theta }_p〉\leftarrow {\mathit{\theta }}_p$;/* Load proposal onto quantum computer. */

$\widehat{p}\leftarrow \text{arg} {\min }_{p=0,\dots ,P}(f(p):=-(z_p+\hspace{0.17em}\log \hspace{0.17em}\pi ({\theta }_p)))$;/* $\mathcal{O}\left(\sqrt{P}\right)$ Algorithm 4 with warm-start at p = 0. */

${\mathit{\theta }}^{(s)}\leftarrow {\mathit{\theta }}_{\widehat{p}}$; return ${\mathit{\theta }}^{(1)},\dots ,{\mathit{\theta }}^{(S)}$ .

Algorithm 6 presents the details of QPMCMC for a continuously-valued target. The algorithm uses conventional (perhaps parallel) computing almost entirely, excluding two lines that are highlighted. The first of these loads proposals onto the quantum machine, and the second line that calls the quantum minimization algorithm presented in Algorithm 4. This sub-routine seeks to find the minimizer of the function(10) $$f(p)=-(z_p+\hspace{0.17em}\log \hspace{0.17em}\pi ({\mathit{\theta }}_p))\hspace{1em}\mathrm{ }\text{over}\hspace{1em}p=0,1,\dots ,P,$$(10) which combines the numerators of (2), the Gumbel-max trick and a trivial negation. As discussed in Section 2.2, the sub-routine requires only $\mathcal{O}\left(\sqrt{P}\right)$ oracle evaluations and, therefore, only $\mathcal{O}\left(\sqrt{P}\right)$ evaluations of the target density $\pi (·)$ using a quantum circuit. In theory, a quantum computer can perform the same computations as a classical computer, but the efficiency of the target evaluations would depend on a number of factors related to the structure of the density itself, the availability of fast quantum random access memory (Giovannetti, Lloyd, and Maccone 2008) and the ability of large numbers of quantum gates to act in concert with negligible noise (Kielpinski, Monroe, and Wineland 2002; Erhard et al. 2019).

In general MCMC, one often calibrates the scaling of a proposal kernel to balance between exploring the target’s domain and remaining within high-density regions. Optimal scaling strategies may lead to a large number of rejected proposals (Roberts and Rosenthal 2001). Indeed, Holbrook (2021) shows that parallel MCMC algorithms are no exception. When using the Gumbel-max trick to sample from proposals, this means that the current state is often quite near optimality, representing a warm-start. Figure 5 shows how this warm-starting effects the number of oracle calls (and hence target evaluations) within the quantum minimization sub-routine over the course of an MCMC trajectory. We target multi-dimensional standard normal distributions of differing dimensions using the vector $(100,\dots ,100)$ as starting position. We fix the number of Gaussian proposals to be 2000 and use standard adaptation procedures (Rosenthal 2011) to target a 50% acceptance rate while guaranteeing convergence. Although no theoretical results validate this target acceptance rate, the rate is close to empirically optimal rates from Holbrook (2021). Across iterations, QPMCMC requires relatively few oracle calls compared to the 2000 target evaluations of parallel MCMC. We also witness the expected drop in the number of oracle evaluations as the chain approaches the target’s high-density region. Remarkably, the algorithm succeeds in sampling from the true discrete distribution in 99.5% of the MCMC iterations. An independent simulation obtains similar results, requiring roughly 7% of the usual number of target evaluations. Figure 6 shows a quantile-quantile plot for all 100 dimensions of a multi-dimensional standard normal distribution. We see no appreciable impact from the rare failure of the quantum minimization sub-routine.

Fig. 5 Total number of oracle evaluations required for each of 2000 quantum parallel MCMC (QPMCMC) iterations for standard multivariate normal targets of five different dimensionalities. Regardless of target dimension, the individual QPMCMC runs require roughly 7% of the usual 4 million target evaluations. Over 99.4% of the 10,000 MCMC iterations across all dimensions successfully sample from the discrete distribution with probabilities of (2). Burn-in iterations require moderately more evaluations because the current state occupies a lower density region and represents a “less good” warm-start.

Fig. 6 Empirical accuracy of quantum parallel MCMC (QPMCMC) for a 100 dimensional spherical Gaussian target. We generate 100,000 samples using 2000 proposals each iteration and remove the first 2000 samples as burn-in. The QQ (quantile-quantile) plot shows sample quantiles adhering closely to the theoretical quantiles. Similar to the independent simulation shown in Figure 5, here QPMCMC requires less than 7.2% of the usual number of target evaluations.

Finally, we compare convergence speed for QPMCMC using 1000, 5000, and 10,000 proposals to target a massively multimodal distribution: a mixture of 1000 two-dimensional, isotropic Gaussians with standard deviations equal to 1 and means equal to $10\times (\mathbf{0},\mathbf{1},\dots ,\mathit{999})$. This target is particularly challenging because the distance between modes is significantly larger than standard deviations. In five independent simulations, we run each algorithm until it achieves and effective sample size of 100. Table 1 contains MCMC iterations and target evaluations required to achieve this effective sample size as well as relative speedups over sequential implementations and relative improvements over the 1000 proposal sampler.

Table 1 Racing to a minimum (between the two dimensions) of 100 effective samples for a target with 1000 disjoint modes.

Proposals	MCMC iterations	Target evaluations	Speedup	Efficiency gain
1000	249,398 (200,998, 311,998)	24,988,963 (20,149,132, 31,265,011)	9.98 (9.98, 9.98)	1
5000	14,398 (12,998, 16,998)	3,314,560 (2,989,418, 3,916,281)	21.72 (21.70, 21.74)	7.58 (6.25, 9.71)
10,000	5998 (4998, 6998)	1,993,484 (1,662,592, 2,330,842)	30 (29.96, 30.26)	12.87 (8.64, 18.80)

NOTE: “Speedup” is ratio between target evaluations required for sequential and quantum implementations. “Efficiency gain” is ratio between target evaluations required for 1000 proposal and 5000/10,000 proposal implementations. We report means (minima, maxima) across five independent runs.

3.3 Discrete-Valued Targets

We wish to sample from a target distribution defined over a discrete set ${\left\{{\mathit{\theta }}_{\alpha }\right\}}_{\alpha \in \mathcal{A}}$, for $\mathcal{A}$ some finite or countably infinite set of indices. Glatt-Holtz et al. (2022) establish the broader measure theoretical foundations for parallel MCMC procedures that use selection probabilities (2) to maintain detailed balance. In particular, when the target distribution has probability mass function $\pi (·)$ defined with respect to the counting measure on the power set $2^{\mathcal{A}}$, then detailed balance results in the kernel $P(·,·)$ satisfying$$\pi (\alpha )=\sum _{\alpha \prime }\pi (\alpha \prime )P(\alpha \prime ,\alpha ) ,\hspace{1em}\forall \alpha ,\alpha \prime \in \mathcal{A} .$$

Sections 3.3.1 and 3.3.2 consider targets defined over finite sets $\mathcal{A}$ and make use of the uniform independence proposal scheme described in Section 3.1, thereby facilitating simplified acceptance probabilities (6). This scheme proposes single-flip updates to the current Markov chain state ${\mathit{\theta }}_0$, but it is worth noting that other multiple-flip schemes would also be amenable to QPMCMC.

3.3.1 Ising Model on a Two-Dimensional Lattice

Here, we are interested in an Ising-type lattice model over configurations $\mathit{\theta }=({\theta }_1,\dots ,{\theta }_D)$ consisting of D individual spins ${\theta }_d\in \{-1,1\}$. In terms of the preceding section, we have $\mathcal{A}={\{-1,1\}}^D$ and $\left|\mathcal{A}\right|=2^D$. In particular, we consider targets of the form(11) $$\pi (\mathit{\theta }|\rho )\propto \hspace{0.17em}\exp \hspace{0.17em}\left(\rho \sum _{(d,d\prime )\in \mathcal{E}}{\theta }_d{\theta }_{d\prime }\right) ,$$(11) where $\rho >0$ is the interaction term and $\mathcal{E}$ is the lattice edge set. We specify a two-dimensional 500-by-500 lattice with nearest neighbors connections and fix ρ = 1. Since this latter setting encourages neighboring spins to equal each other, we begin our QPMCMC trajectories at the lowest-probability “checkerboard” state for an initial configuration. At each iteration, we generate collections of proposal states by uniformly flipping P individual spins θ_p , $p\in \{1,\dots ,P\}$, and obtaining each proposal state ${\mathit{\theta }}_p$ by updating the current state ${\mathit{\theta }}_0$ at the single location corresponding to θ_p . Figure 7 shows results from 10 independent runs using $P\in \{4,8,16,\dots ,2048\}$ proposals. Trace plots of unnormalized log-probabilities indicate that higher proposal counts enable discovery of higher probability configurations. Interestingly, QPMCMC appears to be particularly beneficial in this large P context, requiring less than one-tenth the target evaluations required using conventional computing when P = 2048.

Fig. 7 Convergence of log-posterior for a parallel MCMC sampler targeting an Ising model with uniform interaction ρ = 1 and no external field. All chains begin at minimal probability “checkerboard” state. Larger proposal counts allow discovery of higher probability configurations.

3.3.2 Bayesian Image Analysis

We apply a Bayesian image classification model to an intensity map (Figure 8) of the newly imaged supermassive black hole, Sagittarius A*, located at the Galactic Center of the Milky Way (Akiyama et al. 2022). Whereas one cannot see the black hole itself, one may see the shadow of the black hole cast by the hot, swirling cloud of gas surrounding it. We extract the black hole from the surrounding cloud by modeling the intensity at each of the $D=$4076${}^2=$16,613,776 pixels as belonging to a mixture of two truncated normals with values y_d restricted to the intensity range $[0,255]$. Namely, we follow Hurn (1997) and specify the latent variable model $$\begin{array}{l}{y}_d|({\mu }_{\mathcal{l}},{\sigma }^2,{\theta }_d)\stackrel{\text{ind}}{~}\text{ Normal}({\mu }_{\mathcal{l}},{\sigma }^2),y_d\in [0,255],\\ {\theta }_d=\mathcal{l},d\in \{1,\dots ,D\},\\ {\mu }_{\mathcal{l}}\stackrel{\text{iid}}{~}\text{ Uniform}(0,255),\mathcal{l}\in \{-1,1\},\\ \frac{1}{{\sigma }^2}~\text{ Gamma}\left(\frac{1}{2},\frac{1}{2}\right),\end{array}$$where $\mathit{\theta }=({\theta }_1,\dots ,{\theta }_D)$ share for a prior the Ising model (11) with edges joining neighboring pixels and interaction $\rho =1.2$.

Fig. 8 On the left is a 4076-by-4076 intensity map of the shadow of supermassive black hole Sagittarius A* (Akiyama et al. 2022). On the right is the pixelwise posterior mode of a Bayesian image classification model fit to intensity data. Within a Metropolis-in-Gibbs routine, quantum parallel MCMC using 1024 proposals requires less that one-tenth the posterior evaluations required by conventional parallel MCMC.

We use a QPMCMC-within-Gibbs scheme to infer the join posterior of $\mathit{\theta }$ and the three mixture model parameters. For the former, we use the same QPMCMC scheme as in Section 3.3.1 with 1024 proposals at each iteration. For the latter, we use the closed-form updates presented in Hurn (1997). We run this scheme for 20 million iterations, discarding the first 10 million as burnin. We thin the remaining sample at a ratio of 1 to 40,000 for the latent memberships $\mathit{\theta }$ and 1 to 4000 for the three parameters ${\mu }_{-1}$, μ₁ and ${\sigma }^2$. Using the R package coda (Plummer et al. 2006), we calculate effective sample sizes of the log-posterior (257.1), ${\mu }_{-1}$ (1,578.7), μ₁ (257.6) and ${\sigma }^2$ (2500.0), suggesting adequate convergence. Figure 8 shows both the intensity data and the pixelwise posterior mode of the latent membership vector $\mathit{\theta }$. The QPMCMC requires only 1,977,553,608 target evaluations compared to the 1024 × 20,000,000 $=$ 2.048$\times {10}^{10}$ evaluations required for the analogous parallel MCMC scheme implemented on a conventional computer, representing a 10.36-fold speedup.

4 Discussion

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}\left(\sqrt{P}\right)$. 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 ${\hspace{0.17em}\log \hspace{0.17em}}_2(B)$ qubits with circuit depth $\mathcal{O}(\hspace{0.17em}\log \hspace{0.17em}(B))$. For continuous targets, QPMCMC requires reading $\mathcal{O}(P)$ real vectors ${\mathit{\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=\mathcal{O}(\mathit{PDb})$ classical bits into the quantum machine requires time $\mathcal{O}(\hspace{0.17em}\log \hspace{0.17em}(\mathit{PDb}))$. 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 $\mathcal{O}(\hspace{0.17em}\log \hspace{0.17em}(D))$. Once one has encoded the entire system state, each QPMCMC iteration only requires loading the addresses of the $\mathcal{O}(P)$ proposal states. If one uses b bits to encode each address, then the total time required to load data onto the quantum machine is $\mathcal{O}(\hspace{0.17em}\log \hspace{0.17em}(Pb))$ 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 $\widehat{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-ϵ$. While empirical studies suggest that this error induces little bias, it would be helpful to use this ϵ 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 ϵ by running algorithm $\mathbb{A}$ a mere $\tilde{\mathcal{O}}(\sigma /ϵ)$ times, where $\tilde{\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, De Freitas, and Gordon 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.

Disclosure Statement

There are no competing interests to declare.

Additional information

Funding

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

Appendix A.

Limited Introduction to Quantum Computing

Quantum computers perform operations on unit-length vectors of complex data called quantum bits or qubits. One may write any qubit $\mathit{\psi }$ as a linear combination of the computational basis states $|0〉$ and $|1〉$. In symbols,$$|\psi 〉=\alpha |0〉+\beta |1〉\hspace{1em}\mathrm{ }\text{for}\hspace{1em}\alpha , \beta \in \mathbb{C}\hspace{1em}\mathrm{ }\text{and}\hspace{1em}|\alpha {|}^2+|\beta {|}^2=1 .$$

This formula uses Dirac or bra-ket notation and obscures ideas that are commonplace in the realm of applied statistics. We make the unit-vector specification of $|\psi 〉$ clear by writing$$|0〉=\left[\begin{array}{c}1\\ 0\end{array}\right] ,\hspace{1em}|1〉=\left[\begin{array}{c}0\\ 1\end{array}\right]\hspace{1em}\mathrm{ }\text{or}\hspace{1em}|\psi 〉=\alpha \left[\begin{array}{c}1\\ 0\end{array}\right]+\beta \left[\begin{array}{c}0\\ 1\end{array}\right] .$$

The full machinery of linear algebra is also available within this notation. The conjugate transpose of $|\psi 〉$ is $〈\psi |$. The inner product of $|\psi 〉$ and $|\varphi 〉$ is $〈\varphi |\psi 〉$. The outer product is $|\psi 〉〈\varphi |$. Naturally, we can write $\mathit{\psi }$ as a linear combination of any other such basis elements. Consider instead the vectors$$|+〉=\frac{1}{\sqrt{2}}|0〉+\frac{1}{\sqrt{2}}|1〉\hspace{1em}\mathrm{ }\text{and}\hspace{1em}|-〉=\frac{1}{\sqrt{2}}|0〉-\frac{1}{\sqrt{2}}|1〉 .$$

A little algebra shows that $|+〉$ and $|-〉$ are indeed unit-length and orthogonal to each other. A little more algebra reveals that, with respect to this basis, $\mathit{\psi }$ has coefficients $\alpha \prime$ and $\beta \prime$ given by $(\alpha +\beta )/\sqrt{2}$ and $(\alpha -\beta )/\sqrt{2}$, respectively. A last bit of algebra shows that this representation is consistent with $\mathit{\psi }$’s unit length.

But linear algebra is not the only aspect of quantum computing that should come easily to applied statisticians. In addition to thinking of $\mathit{\psi }$ as a vector, it is also useful to think of $\mathit{\psi }$ as a (discrete) probability distribution over the computational basis states $|0〉$ and $|1〉$. The constraints on coefficients α and β mean that we can think of $|\alpha {|}^2$ and $|\beta {|}^2$ as probabilities that $\mathit{\psi }$ is in state $|0〉$ or $|1〉$, respectively. Accordingly, $|+〉$ and $|-〉$ encode uniform distributions over the computational basis states. In the parlance of quantum mechanics, α, β and $\pm 1/\sqrt{2}$ are all probability amplitudes, and $\mathit{\psi },|+〉$ and $|-〉$ are all superpositions of the computational basis states. Quantum measurement of $\mathit{\psi }$ results in $|0〉$ with probability $|\alpha {|}^2$, but—in the following—we can safely think of this physical procedure as drawing a single discrete sample from $\mathit{\psi }$’s implied probability distribution.

Quantum logic gates perform linear operations on qubits like $\mathit{\psi }$ and take the form of unitary matrices U satisfying ${\mathbf{U}}^{†}\mathbf{U}=\mathbf{I}$ for ${\mathbf{U}}^{†}$ the conjugate transpose of U. One terrifically important single-qubit quantum gate is the Hadamard gate $$\mathbf{H}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}\hfill 1& \hfill 1\\ \hfill 1& \hfill -1\end{array}\right] .$$

One may verify that $\mathbf{H}$ is indeed unitary and that its action maps $|0〉$ to $|+〉$ and $|1〉$ to $|-〉$. In fact, the reverse is also true on account of the symmetry of $\mathbf{H}$. The Hadamard gate thus takes the computational basis states in and out of superposition, facilitating a phenomenon called quantum parallelism. Given a function $f:\{0,1\}\to \{0,1\}$, consider the two-qubit quantum oracle gate ${\mathbf{U}}_f$ which takes $|x〉|y〉$ as input and returns $|x〉|y\oplus f(x)〉$ as output, where ⊕ denotes addition modulo 2. Importantly, the output simplifies to $|x〉|f(x)〉$ for y = 0. We now consider a quantum circuit that acts on two qubits by first applying the Hadamard transform $\mathbf{H}$ to the first qubit and then applying the oracle gate ${\mathbf{U}}_f$ to both. Using the state $|0〉|0〉$ as input, we have(A.1) $$|0〉|0〉\to \frac{1}{\sqrt{2}}|0〉|0〉+\frac{1}{\sqrt{2}}|1〉|0〉\to \frac{1}{\sqrt{2}}|0〉|f(0)〉+\frac{1}{\sqrt{2}}|1〉|f(1)〉 .$$(A.1)

The quantum circuit evaluates $f(·)$ simultaneously over both inputs! Unfortunately, the scientist who implements this circuit cannot directly access ${\mathbf{U}}_f$’s output, and measurement will provide only $|0〉|f(0)〉$ or $|1〉|f(1)〉$ with probability 1/2 each. Unlocking the potential of quantum parallelism requires more ingenuity.

Uncountably many single- and two-qubit quantum gates exist, but the real power of quantum computing stems from the development of complex quantum gates that act on multiple qubits simultaneously. To access this power, we need one more tool that also appears in the statistics literature. The Kronecker or tensor product between an L-by-M matrix A and an N-by-O matrix B is the LN-by-MO matrix(A.2) $$\mathbf{A}\otimes \mathbf{B}=\left[\begin{array}{ccc}{A}_{11}\mathbf{B}& \dots & A_{1M}\mathbf{B}\\ ⋮& \ddots & ⋮\\ A_{L1}\mathbf{B}& \cdots & A_{LM}\mathbf{B}\end{array}\right] .$$(A.2)

In statistics, the Kronecker product features in the definition of a matrix normal distribution and is sometimes helpful when specifying the kernel function of a multivariate Gaussian process (Werner, Jansson, and Stoica 2008). Here, the product is equivalent to the parallel action of individual quantum gates on individual qubits. Simple application of Formula (A.2) shows that(A.3) $$\begin{array}{c}{\mathbf{H}}^{\otimes 2}=\mathbf{H}\otimes \mathbf{H}=\frac{1}{2}\left[\begin{array}{cccc}\hfill 1& \hfill 1& \hfill 1& \hfill 1\\ \hfill 1& \hfill -1& \hfill 1& \hfill -1\\ \hfill 1& \hfill 1& \hfill -1& \hfill -1\\ \hfill 1& \hfill -1& \hfill -1& \hfill 1\end{array}\right]\hspace{1em}\mathrm{ }\text{and}\\ |0〉\otimes |0〉=\left[\begin{array}{c}1\\ 0\end{array}\right]\otimes \left[\begin{array}{c}1\\ 0\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right].\end{array}$$(A.3)

We may also denote the left product ${\mathbf{H}}^{\otimes 2}$ and the right product $|0{〉}^{\otimes 2},|00〉$ or $|0〉|0〉$. One may therefore express (A.1) as a series of transformations applied to the 4-vector on the very right. Letting $|10〉,|01〉$ and $|11〉$ take on analogous meanings to $|00〉$, an immediate result of (A.3) is that(A.4) $${\mathbf{H}}^{\otimes 2}|00〉=\frac{1}{2}(|00〉+|01〉+|10〉+|11〉) .$$(A.4)

The action of ${\mathbf{H}}^{\otimes 2}$ transforms $|00〉$ into a superposition of the states $|00〉,|10〉,|01〉$ and $|11〉$, where the probability of selecting each is a uniform ${(1/2)}^2=1/4$. Writing so many 0s and 1s is tedious, so an alternative notation becomes preferable. Exchanging $|0〉$ for $|00〉,|1〉$ for $|01〉,|2〉$ for $|10〉$, and $|3〉$ for $|11〉$, (A.4) becomes the more succinct$${\mathbf{H}}^{\otimes 2}|00〉=\frac{1}{2}(|0〉+|1〉+|2〉+|3〉) .$$

This formula extends generally to operations over n qubits. Now letting $N=2^n$, we have$${\mathbf{H}}^{\otimes n}|0{〉}^{\otimes n}=\frac{1}{\sqrt{N}}(|0〉+|1〉+\dots +|N-1〉)=\frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}|x〉=:|h〉 ,$$and we call $|h〉$ a uniform superposition over the states $|0〉,\dots ,|N-1〉$. The many-qubit analogue for the quantum parallelism of (A.1) is then$$\left|0{〉}^{\otimes n}\right|0〉\to \left(\frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}|x〉\right)|0〉\to \frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}|x〉|f(x)〉 .$$

Appendix B.

Gumbel-max

We wish to randomly select a single element from the set $\{0,1,\dots ,P\}$ with probability proportional to the unnormalized probabilities ${\mathit{\pi }}^{*}=({\pi }_0^{*},{\pi }_1^{*},\dots ,{\pi }_P^{*})$. That is, there exists a c > 0 such that ${\mathit{\pi }}^{*}=c\mathit{\pi }$, for $\mathit{\pi }$ a valid probability vector, but we only have access to ${\mathit{\pi }}^{*}$. Define $\mathit{\lambda }:=\hspace{0.17em}\log \hspace{0.17em}\mathit{\pi }$ and ${\mathit{\lambda }}^{*}:=\hspace{0.17em}\log \hspace{0.17em}{\mathit{\pi }}^{*}=\hspace{0.17em}\log \hspace{0.17em}\mathit{\pi }+\hspace{0.17em}\log \hspace{0.17em}c$, and assume that $z_0,\dots ,z_P\stackrel{\text{iid}}{\sim }G\mathit{umbel}(0,1)$. Then, the probability density function $g(·)$ and cumulative distribution function $G(·)$ for each individual z_p are$$g(z_p)=\hspace{0.17em}\exp \hspace{0.17em}(-z_p-\hspace{0.17em}\exp \hspace{0.17em}(-z_p))\hspace{1em}\mathrm{ }\text{and}\hspace{1em}G(z_p)=\hspace{0.17em}\exp \hspace{0.17em}(-\hspace{0.17em}\exp \hspace{0.17em}(-z_p)) .$$

Now, defining the random variables ${\alpha }_p^{*}:={\lambda }_p^{*}+z_p,{\alpha }_p:={\lambda }_p+z_p$ and$$\widehat{p}=\text{arg} {\max }_{p=0,\dots ,P}{\alpha }_p^{*} ,$$we have the result$$\mathrm{ }\text{Pr}(\widehat{p}=p)={\pi }_p .$$

In words, the procedure of adding Gumbel noise to unnormalized log-probabilities and taking the index of the maximum produces a random variable that follows the discrete distribution over $\mathit{\pi }$. Moving from left to right:$$\begin{array}{c}\mathrm{ }\text{Pr}(\widehat{p}=p)=\mathrm{ }\text{Pr}({\alpha }_p^{*}>{\alpha }_{p\prime }^{*}, \forall p\prime \ne p)\\ =\mathrm{ }\text{Pr}({\alpha }_p+\hspace{0.17em}\log \hspace{0.17em}c>{\alpha }_{p\prime }+\hspace{0.17em}\log \hspace{0.17em}c, \forall p\prime \ne p)\\ =\mathrm{ }\text{Pr}({\alpha }_p>{\alpha }_{p\prime }, \forall p\prime \ne p)\\ ={\int }_{-\infty }^{\infty }\prod _{p\prime \ne p}\mathrm{ }\text{Pr}({\alpha }_p>{\alpha }_{p\prime }|{\alpha }_p)g({\alpha }_p-{\lambda }_p) \mathrm{d}{\alpha }_p\\ ={\int }_{-\infty }^{\infty }\prod _{p\prime \ne p}G({\alpha }_p-{\lambda }_{p\prime })g({\alpha }_p-{\lambda }_p) \mathrm{d}{\alpha }_p\\ ={\int }_{-\infty }^{\infty }\prod _{p\prime \ne p}\hspace{0.17em}\exp \hspace{0.17em}(-\hspace{0.17em}\exp \hspace{0.17em}({\lambda }_{p\prime }-{\alpha }_p))\\ \hspace{1em}\hspace{0.17em}\exp \hspace{0.17em}(-{\alpha }_p+{\lambda }_p-\hspace{0.17em}\exp \hspace{0.17em}(-{\alpha }_p+{\lambda }_p)) \mathrm{d}{\alpha }_p .\end{array}$$

Recalling that ${\lambda }_{p\prime }=\hspace{0.17em}\log \hspace{0.17em}{\pi }_{p\prime }$, we exponentiate the logarithms, and the integral becomes$$\begin{array}{c}{\pi }_p{\int }_{-\infty }^{\infty }\prod _{p\prime \ne p}\hspace{0.17em}\exp \hspace{0.17em}(-{\pi }_{p\prime }\hspace{0.17em}\exp \hspace{0.17em}(-{\alpha }_p))\hspace{0.17em}\exp \hspace{0.17em}(-{\alpha }_p-{\pi }_p\hspace{0.17em}\exp \hspace{0.17em}(-{\alpha }_p)) \mathrm{d}{\alpha }_p\\ \hspace{1em}={\pi }_p{\int }_{-\infty }^{\infty }\hspace{0.17em}\exp \hspace{0.17em}(-{\alpha }_p)\hspace{0.17em}\exp \hspace{0.17em}(-\sum _{p\prime =0}^P{\pi }_{p\prime }\hspace{0.17em}\exp \hspace{0.17em}(-{\alpha }_p)) \mathrm{d}{\alpha }_p\\ \hspace{1em}={\pi }_p{\int }_{-\infty }^{\infty }\hspace{0.17em}\exp \hspace{0.17em}(-{\alpha }_p)\hspace{0.17em}\exp \hspace{0.17em}(-\hspace{0.17em}\exp \hspace{0.17em}(-{\alpha }_p)) \mathrm{d}{\alpha }_p={\pi }_p ,\end{array}$$where the final equality follows easily from a change of variables.

Appendix C.

Mixing of Parallel MCMC and QPMCMC

To ascertain whether QPMCMC mixes differently compared to conventional parallel MCMC, we run both algorithms for a range of proposal counts to sample from a 10-dimensional standard normal distribution. For each algorithm and proposal setting, we run 100 independent chains for 10,000 iterations and obtain effective sample sizes ESS_d for $d\in \{1,\dots ,10\}$. We then calculate the relative differences between the means and minima of one chain generated using parallel MCMC and one chain generated using QPMCMC; for example:$$\begin{array}{c}\text{Relative}\mathrm{ }\text{difference}\mathrm{ }\text{between}\mathrm{ }\text{means\hspace{0.17em}}{\overline{\mathrm{ }\text{ESS}}}_{(·)}\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}=\frac{\left|{\overline{\mathrm{ }\text{ESS}}}_{\text{pMCMC}}-{\overline{\mathrm{ }\text{ESS}}}_{\text{QPMCMC}}\right|}{{\overline{\mathrm{ }\text{ESS}}}_{\text{pMCMC}}}\end{array}$$

Figure B.1 shows results. In general, the majority relative differences are small. For both statistics, mean relative differences are less than 0.05, regardless of proposal count. Again for both statistics, more than 75% of the independent runs result in relative differences below 0.1. We note that relative differences between means (blue) appear to decrease with the number of proposals, but the same does not hold for relative differences between minima (red).

Fig. B.1 Relative differences between effective sample sizes (ESS) for parallel MCMC (pMCMC) and quantum parallel MCMC (QPMCMC) across a range of proposal counts. We target a 10-dimensional standard normal distribution. For each algorithm and each proposal setting, we generate 100 independent chains of length 10,000 and calculate the mean and minimum ESS across dimensions.

Appendix D.

Note on Simulations

We use R (R Core Team 2021), Python (van Rossum 1995), TensorFlow (Abadi et al. 2016) and TensorFlow Probability MCMC (Lao et al. 2020) in our simulations and the ggplot2 package to generate all figures (Wickham 2016). In R, we use the package coda (Plummer et al. 2006) to calculate effective sample sizes. In Python, we use a built-in function from TensorFlow Probability MCMC. The primary color palette comes from Ram and Wickham (2018).

All simulations rely on the fact that the Grover iterations of (3) manipulate the probability amplitudes in a deterministic manner. For example, the following R code specifies N uniform probability amplitudes (sqrtProbs), performs I Grover iterations and measures a single state according to the resulting probabilities.

sqrtProbs <- rep(1/sqrt(N),N); i <- 1 while (i < = I) { sqrtProbs <- (1-2*marks)*sqrtProbs sqrtProbs <- -sqrtProbs + 2*mean (sqrtProbs) i <- i + 1 } measurement <- sample(size = 1, x = 1:N, prob = sqrtProbsˆ2)

Of course, simulating these iterations on a classical computer requires precomputing the values of marks beforehand. All code is available online at https://github.com/andrewjholbrook/qpMCMC.