Quantum speedup of Monte Carlo methods

Monte Carlo methods use random sampling to estimate numerical quantities which are hard to compute deterministically. One important example is the use in statistical physics of rapidly mixing Markov chains to approximately compute partition functions. In this work, we describe a quantum algorithm which can accelerate Monte Carlo methods in a very general setting. The algorithm estimates the expected output value of an arbitrary randomized or quantum subroutine with bounded variance, achieving a near-quadratic speedup over the best possible classical algorithm. Combining the algorithm with the use of quantum walks gives a quantum speedup of the fastest known classical algorithms with rigorous performance bounds for computing partition functions, which use multiple-stage Markov chain Monte Carlo techniques. The quantum algorithm can also be used to estimate the total variation distance between probability distributions efficiently.


Introduction
Monte Carlo methods are now ubiquitous throughout science, in fields as diverse as statistical physics [1], microelectronics [2] and mathematical finance [3]. These methods use randomness to estimate numerical properties of systems which are too large or complicated to analyse deterministically. In general, the basic core of Monte Carlo methods involves estimating the expected output value μ of a randomized algorithm A. The natural algorithm for doing so is to produce k samples, each corresponding to the output of an independent execution of A, and then to output the averageμ of the samples as an approximation of μ. Assuming that the variance of the random variable corresponding to the output of A is at most σ 2 , the (a) Prior work The topic of quantum estimation of mean output values of algorithms with bounded variance connects to several previously explored directions. First, it generalizes the problem of approximating the mean, with respect to the uniform distribution, of an arbitrary bounded function. This has been addressed by a number of authors. The first asymptotically optimal quantum algorithm for this problem, which uses O(1/ ) queries to achieve additive error , seems to have been given by Heinrich [16]; an elegant alternative optimal algorithm was later presented by Brassard et al. [17]. Using similar techniques to Brassard et al., Wocjan et al. [18] described an efficient algorithm for estimating the expected value of an arbitrary bounded observable. It is not difficult to combine these ideas to approximate the mean of arbitrary bounded functions with respect to non-uniform distributions (see §2).
One of the main technical ingredients in this paper is based on an algorithm of Heinrich for approximating the mean, with respect to the uniform distribution, of functions with bounded L 2 norm [16]. Here, we describe a generalization of this result to non-uniform distributions, using similar techniques. This is roughly analogous to the way that amplitude amplification [19] generalizes Grover's quantum search algorithm [20]. The related problem of quantum estimation of expectation values of observables, an important task in the simulation of quantum systems, has been studied by Knill et al. [21]. These authors give an algorithm for estimating tr(Aρ) for observables A such that one can efficiently implement the operator e −iAt . The algorithm is efficient (i.e. achieves runtimes close to O(1/ )) when the tails of the distribution tr(Aρ) decay quickly. However, in the case where one only knows an upper bound on the variance of this distribution, the algorithm does not achieve a better runtime than classical sampling.
Quantum algorithms have been used previously to approximate classical partition functions and solve related problems. In particular, a number of authors (see [22] and references therein) have considered the complexity of computing Ising and Potts model partition functions. These works in some cases achieve exponential quantum speedups over the best-known classical algorithms. Unfortunately, they in general either produce an approximation accurate up to a specified additive error bound, or only work for specific classes of partition function problems with restrictions on interaction strengths and topologies, or both. Here, we aim to approximate partition functions up to small relative error in a rather general setting.
Using related techniques to the present work, Somma et al. [23] used quantum walks to accelerate classical simulated annealing processes, and quantum estimation of partition functions up to small relative error was addressed by Wocjan et al. [18]. Their algorithm, which is based on the use of quantum walks and amplitude estimation, achieves a quadratic speedup over classical algorithms with respect to both mixing time and accuracy. However, it cannot be directly applied to accelerate the most efficient classical algorithms for approximating partition function problems, which use so-called Chebyshev cooling schedules (discussed in §3). This is essentially because these algorithms are based around estimating the mean of random variables given only a bound on the variance. This was highlighted as an open problem in [18], which we resolve here.
Several recent works have developed quantum algorithms for the quantum generalization of sampling from a Gibbs distribution: producing a Gibbs state ρ ∝ e −βH for some quantum Hamiltonian H [24][25][26][27]. Given such a state, one can measure a suitable observable to compute some quantity of interest about H. Supplied with an upper bound on the variance of such an observable, the procedure detailed here can be used (as for any other quantum algorithm) to reduce the number of repetitions required to estimate the observable to a desired accuracy.

(b) Techniques
We now give an informal description of our algorithms, which are summarized in table 1 (for technical details and proofs, see §2). For any randomized or quantum algorithm A, we write v(A) for the random variable corresponding to the value computed by A, with the expected value of For concreteness, we think of A as a quantum algorithm which operates on n qubits, each initially in the state |0 , and whose quantum part finishes with a measurement of k of the qubits in the computational basis. Given that the measurement returns outcome x ∈ {0, 1} k , the final output is then φ(x), for some fixed function φ : {0, 1} k → R. If A is a classical randomized algorithm, or a quantum circuit using (for example) mixed states and intermediate measurements, a corresponding unitary quantum circuit of this form can be produced using standard reversiblecomputation techniques [28]. As is common in works based on quantum amplitude amplification and estimation [19], we also assume that we have the ability to execute the algorithm A −1 , which is the inverse of the unitary part of A. If we do have a description of A as a quantum circuit, this can be achieved simply by running the circuit backwards, replacing each gate with its inverse.
We first deal with the special case where the output of A is bounded between 0 and 1. Here, a quantum algorithm for approximating μ := E[v(A)] quadratically faster than is possible classically can be found by combining ideas from previously known algorithms [16][17][18]. We append an additional qubit and define a unitary operator W on k + 1 qubits which performs the map |x |0 → |x If the final measurement of the algorithm A is replaced with performing W, then measuring the added qubit, the probability that we receive the answer 1 is precisely μ. Using quantum amplitude estimation [19], the probability that this measurement  returns 1 can be estimated to higher accuracy than is possible classically. Using t iterations of amplitude estimation, we can output an estimateμ such that |μ − μ| = O( √ μ/t + 1/t 2 ) with high probability [19]. In particular, O(1/ ) iterations of amplitude estimation are sufficient to produce an estimateμ such that |μ − μ| ≤ with, say, 99% probability. The next step is to use the above algorithm as a subroutine in a more general procedure that can deal with algorithms A whose output is non-negative, has bounded 2 norm, but is not necessarily bounded between 0 and 1. That is, algorithms for which we can control the expression v(A) 2 The procedure for this case generalizes and is based on the same ideas as a previously known result for the uniform distribution [16].
The idea is to split the output of A up into disjoint intervals depending on size. Write A p,q for the 'truncated' algorithm which outputs v(A) if p ≤ v(A) < q, and otherwise outputs 0. We estimate μ by applying the above algorithm to estimate E[v(A p,q )] for a sequence of O(log 1/ ) intervals which are exponentially increasing in size, and summing the results. As the intervals [p, q) get larger, the accuracy with which we approximate E[v(A p,q )] decreases, and values v(A) larger than about 1/ are ignored completely. However, the overall upper bound on v(A) 2 allows us to infer that these larger values do not affect the overall expectation μ much; indeed, if μ depended significantly on large values in the output, the 2 norm of v(A) would be high.
The final result is that for v(A) 2 = O(1), given appropriate parameter choices, the estimateμ satisfies |μ − μ| = O( ) with high probability, and the algorithm uses AÕ(1/ ) times in total. This scaling is a near-quadratic improvement over the best possible classical algorithm.
We next consider the more general case of algorithms A which have bounded variance, but whose output need not be non-negative, nor bounded in 2 norm. To apply the previous algorithm, we would like to transform the output of A to make its 2 norm low. If v(A) has mean μ and variance upper-bounded by σ 2 , a suitable way to achieve this is to subtract μ from the output of A, then divide by σ . The new algorithm's output would have 2 norm upper-bounded by 1, and estimating its expected value up to additive error /σ would give us an estimate of μ up to . Unfortunately, we of course do not know μ initially, so cannot immediately implement this idea. To approximately implement it, we first run A once and use the outputm as a proxy for μ. Because Var(v(A)) ≤ σ 2 ,m is quite likely to be within distance O(σ ) of μ. Therefore, the algorithm B produced from A by subtractingm and dividing by σ is quite likely to have 2 norm upper-bounded by a constant. We can thus efficiently estimate the positive and negative parts of E[v(B)] separately, then combine and rescale them. The overall algorithm achieves accuracy in timeÕ(σ/ ). For a more precise statement, see theorem 2.5.
A similar idea can be used to approximate the expected output value of algorithms for which we have a bound on the relative variance, namely that Var(v(A)) = O(μ 2 ). In this setting, it turns out thatÕ(1/ ) uses of A suffice to produce an estimateμ accurate up to relative error , i.e. for which |μ − μ| ≤ μ. This is again a near-quadratic improvement over the best possible classical algorithm. See theorem 2.6 for the details.

(c) Approximating partition functions
In this section, we discuss (with details in §3) how these algorithms can be applied to the problem of approximating partition functions. Consider a (classical) physical system which has state space Ω, together with a Hamiltonian H : Ω → R specifying the energy of each configuration 2 x ∈ Ω. Here, we will assume that H takes integer values in the set {0, . . . , n}. A central problem is to compute the partition function for some inverse temperature β defined by β = 1/(k B T), where T is the temperature and k B is Boltzmann's constant. As well as naturally encapsulating various models in statistical physics, such as the Ising and Potts models, this framework also encompasses well-studied problems in computer science, such as counting the number of valid k-colourings of a graph. In particular, Z(∞) counts the number of configurations x such that H(x) = 0. It is often hard to compute Z(β) for large β but easy to approximate Z(β) ≈ |Ω| for β ≈ 0. In many cases, such as the Ising model, it is known that computing Z(∞) exactly falls into the #P-complete complexity class [29], and hence is unlikely to admit an efficient quantum or classical algorithm. Here, our goal will be to approximate Z(β) up to relative error , for some small . That is, to outputZ such that |Z − Z(β)| ≤ Z(β), with high probability. For simplicity, we will focus on β = ∞ in the following discussion, but it is easy to see how to generalize to arbitrary β.
Let 0 = β 0 < β 1 < · · · < β = ∞ be a sequence of inverse temperatures. A standard classical approach to design algorithms for approximating partition functions [4,18,[30][31][32] is based around expressing Z(β ) as the telescoping product If we can compute Z(β 0 ) = |Ω| and can also approximate each of the ratios α i := Z(β i+1 )/Z(β i ) accurately, taking the product will give a good approximation to Z(β ). Let π i denote the Gibbs (or Boltzmann) probability distribution corresponding to inverse temperature β i , where To approximate α i , we define the random variable Then one can readily compute that E π i [Y i ] = α i , so sampling from each distribution π i allows us to estimate the quantities α i . It will be possible to estimate α i up to small relative error efficiently This motivates the concept of a Chebyshev cooling schedule [4]: It is known that, for any partition function problem as defined above such that |Ω| = A, there exists a Chebyshev cooling schedule with =Õ log A [4].
It is sufficient to approximate E[Y i ] up to relative error O( / ) for each i to get an overall approximation accurate up to relative error . To achieve this, the quantum algorithm presented here needs to use at mostÕ( / ) samples from Y i . Given a Chebyshev cooling schedule with =Õ log A , the algorithm thus usesÕ((log A)/ ) samples in total, a near-quadratic improvement in terms of over the complexity of the fastest known classical algorithm [4].
In general, we cannot exactly sample from the distributions π i . Classically, one way of approximately sampling from these distributions is to use a Markov chain which mixes rapidly and has stationary distribution π i . For a reversible, ergodic Markov chain, the time required to produce such a sample is controlled by the relaxation time τ := 1/(1 − |λ 1 |) of the chain, where λ 1 is the second largest eigenvalue in absolute value [8]. In particular, sampling from a distribution close to π i in total variation distance requires Ω(τ ) steps of the chain.
It has been known for some time that quantum walks can sometimes mix quadratically faster [10]. One case where efficient mixing can be obtained is for sequences of Markov chains whose stationary distributions π are close [11]. Further, for this special case, one can approximately produce coherent 'quantum sample' states |π = x∈Ω π (x)|x efficiently. Here, we can show ( §3) that the Chebyshev cooling schedule condition implies that each distribution in the sequence π 1 , . . . , π −1 is close enough to its predecessor that we can use techniques of Wocjan & Abeyesinghe [11] to approximately produce any state |π i usingÕ( √ τ ) quantum walk steps each. Using similar ideas, we can approximately reflect about |π i using onlyÕ( √ τ ) quantum walk steps.
Approximating E[Y i ] up to relative error O( / ) using our algorithm requires one quantum sample approximating |π i , andÕ( / ) approximate reflections about |π i . Therefore, the total number of quantum walk steps required for each i isÕ( √ τ / ). Summing over i, we get a quantum algorithm for approximating an arbitrary partition function up to relative error usingÕ((log A) √ τ / ) quantum walk steps. The fastest known classical algorithm [4] exhibits quadratically worse dependence on both τ and .
In the above discussion, we have neglected the complexity of computing the Chebyshev cooling schedule itself. An efficient classical algorithm for this task is known [4], which runs in timeÕ((log A)τ ). Adding the complexity of this part, we finish with an overall complexity of We leave the interesting question open of whether there exists a more efficient quantum algorithm for finding a Chebyshev cooling schedule.

(d) Applications
We now sketch several representative settings (for details, see §3) in which our algorithm for approximating partition functions gives a quantum speedup.
-The ferromagnetic Ising model above the critical temperature. This well-studied statistical physics model is defined in terms of a graph The Markov chain known as the Glauber dynamics is known to mix rapidly above a certain critical temperature and to have as its stationary distribution the Gibbs distribution. For example, for any graph with maximum degree O(1), the mixing time of the Glauber dynamics for sufficiently low inverse temperature β is O(n log n) [33]. In this case, as A = 2 n , the quantum algorithm approximates Z(β) to within relative error inÕ(n 3/2 / + n 2 ) steps. The corresponding classical algorithm [4] usesÕ(n 2 / 2 ) steps. -Counting colourings. Here, we are given a graph G with n vertices and maximum degree d.
We seek to approximately count the number of valid k-colourings of G, where a colouring of the vertices is valid if all pairs of neighbouring vertices are assigned different colours.
In the case where k > 2d, the use of a rapidly mixing Markov chain gives a quantum algorithm approximating the number of colourings of G up to relative error in timẽ O(n 3/2 / + n 2 ), as compared with the classicalÕ(n 2 / 2 ) [4]. -Counting matchings. A matching in a graph G is a subset M of the edges of G such that no pair of edges in M shares a vertex. In statistical physics, matchings are studied under the name of monomer-dimer coverings [34]. Our algorithm can approximately count the number of matchings on a graph with n vertices and m edges inÕ(n 3/2 m 1/2 / + n 2 m) steps, as compared with the classicalÕ(n 2 m/ 2 ) [4].
Finally, as another example of how our algorithm can be applied, we improve the accuracy of an existing quantum algorithm for estimating the total variation distance between probability distributions. In this setting, we are given the ability to sample from probability distributions p and q on n elements, and would like to estimate the distance between them up to additive error . A quantum algorithm of Bravyi, Harrow and Hassidim solves this problem using O( √ n/ 8 ) samples [15], while no classical algorithm can achieve sublinear dependence on n [35].
Quantum mean estimation can significantly improve the dependence of this quantum algorithm on . The total variation distance between p and q can be described as the expected value of the random variable

Algorithms
We now give technical details, parameter values and proofs for the various algorithms described informally in §1. Recall that, for any randomized or quantum algorithm A, we let v(A) be the random variable corresponding to the value computed by A. We assume that A takes no input directly, but may have access to input (e.g. via queries to some black box or 'oracle') during its execution. We further assume throughout that A is a quantum algorithm of the following form: apply some unitary operator to the initial state |0 n ; measure k ≤ n qubits of the resulting state in the computational basis, obtaining outcome x ∈ {0, 1} k ; output φ(x) for some easily computable function φ : {0, 1} k → R. We finally assume that we have access to the inverse of the unitary part of the algorithm, which we write as A −1 .
The following simple and well-known result, sometimes known as the powering lemma, will be useful to us in various contexts:

Lemma 2.1 (Powering lemma [36]). Let
A be a (classical or quantum) algorithm which aims to estimate some quantity μ, and whose outputμ satisfies |μ −μ| ≤ except with probability γ , for some fixed γ < 1 2 . Then, for any δ > 0, it suffices to repeat A O(log 1/δ) times and take the median to obtain an estimate which is accurate to within with probability at least 1 − δ.
We will also need the following fundamental result from [19]: Theorem 2.2 (Amplitude estimation [19]). There is a quantum algorithm called amplitude estimation which takes as input one copy of a quantum state |ψ , a unitary transformation U = 2|ψ ψ| − I, a unitary transformation V = I − 2P for some projector P, and an integer t. The algorithm outputsã, an estimate of a = ψ|P|ψ , such that with probability at least 8/π 2 , using U and V t times each.
The success probability of 8/π 2 can be improved to 1 − δ for any δ > 0 using the powering lemma at the cost of an O(log 1/δ) multiplicative factor.

(a) Estimating the mean with bounded output values
We first consider the problem of estimating E[v(A)] in the special case where v(A) is bounded between 0 and 1. The algorithm for this case (described as algorithm 1) is effectively a combination of elegant ideas of Brassard et al. [17] and Wocjan et al. [18]. The former described an algorithm for efficiently approximating the mean of an arbitrary function with respect to the uniform distribution; the latter described an algorithm for approximating the expected value of a particular observable, with respect to an arbitrary quantum state. The first quantum algorithm achieving optimal scaling for approximating the mean of a bounded function under the uniform distribution was due to Heinrich [16].
with probability at least 1 − δ, where C is a universal constant. In particular, for any fixed δ > 0 and any such that 0 ≤ ≤ 1, to produce an estimateμ such that with probability at least

Input: an algorithm
Assume that A is a quantum algorithm which makes no measurements until the end of the algorithm; operates on initial input state |0 n ; and its final measurement is a measurement of the last k ≤ n of these qubits in the computational basis.
(i) Let W be the unitary operator on k + 1 qubits defined by where each computational basis state x ∈ {0, 1} k is associated with a real number φ(x) ∈ [0, 1] such that φ(x) is the value output by A when measurement outcome x is received.
Proof. The complexity claim follows immediately from theorem 2.2. Also observe that W can be implemented efficiently, as it is a controlled rotation of one qubit dependent on the value of φ(x) [18]. It remains to show the accuracy claim. The final state of A, just before its last measurement, can be written as for some normalized states |ψ x . If we then attach an ancilla qubit and apply W, we obtain where P = I ⊗ |1 1|. Therefore, when we apply amplitude estimation, by theorem 2.2, we obtain with probability at least 8/π 2 . The powering lemma (lemma 2.1) implies that the median of O(log 1/δ) repetitions will lie within this accuracy bound with probability at least 1 − δ.
Observe that U = 2|ψ ψ| − I can be implemented with one use each of A and A −1 , and V = I − 2P is easy to implement.
It seems likely that the median-finding algorithm of Nayak & Wu [7] could also be generalized in a similar way, to efficiently compute the median of the output values of any quantum algorithm. As we will not need this result here, we do not pursue this further.

(b) Estimating the mean with bounded 2 norm
We now use algorithm 1 to give an efficient quantum algorithm for approximating the mean output value of a quantum algorithm whose output has bounded 2 norm. In what follows, for  All of the uses of algorithm 1 succeed, except with probability at most 1 5 in total. To estimate the total error in the case where they all succeed, we write and use the triangle inequality term by term to obtain Let p(x) denote the probability that A outputs x. We have and similarly So the total error is at most We apply Cauchy-Schwarz to the first part of each term in the sum where the second inequality follows from Inserting this bound and using E[v(A 0,1 )] ≤ 1, we obtain Inserting the definitions of t 0 and k, we get an overall error bound using 0 < < 1 2 in the second inequality. For a sufficiently large constant D, this is upper-bounded by ( v(A) 2 + 1) 2 as claimed. (1), to achieve additive error the number of uses of A that we need is O((1/ ) log 3/2 (1/ ) log log(1/ )). By the powering lemma, we can repeat algorithm 2 O(log 1/δ) times and take the median to improve the probability of success to 1 − δ for any δ > 0.

(c) Estimating the mean with bounded variance
We are now ready to formally state our algorithm for estimating the mean output value of an arbitrary algorithm with bounded variance, as algorithm 3. For clarity, some of the steps are reordered as compared with the informal description in §1. Recall that, in the classical setting, if we wish to estimate E[v(A)] up to additive error for an arbitrary algorithm A such that Input: an algorithm A such that Var(v(A)) ≤ σ 2 for some known σ , and an accuracy such that < 4σ . . We therefore assume that |m − μ | ≤ 3. In this case, we have by linearity of expectation, using a union bound we have that σμ approximates E[v(A)] up to additive error with probability at least 2 3 .

(d) Estimating the mean with bounded relative error
It is often useful to obtain an estimate of the mean output value of an algorithm which is accurate up to small relative error, rather than the absolute error achieved by algorithm 3. Assume that we have the bound on the relative variance that Var Proof. The complexity bounds follow from lemma 2.4; we now analyse the claim about accuracy. m is a random variable whose expectation is E[v(A)] and whose variance is Var(v(A))/ 32B .  By Chebyshev's inequality, we have We Multiplying bym and taking a union bound, we get an estimate of E[v(A)] which is accurate up to except with probability at most 1 4 .
Once again, using the powering lemma, we can repeat algorithms 3 and 4 O(log 1/δ) times and take the median to improve their probabilities of success to 1 − δ for any δ > 0. Algorithm 4 can be extended to work for subroutines A which output both positive and negative values in a similar way to algorithm 3, by modifying step (ii) of the algorithm to estimate and recombine the positive and negative parts of the output of A/|m|. We omit the details as this variant is not required for the applications below.
To see that algorithms 3 and 4 are close to optimal, we can appeal to a result of Nayak & Wu [7]. Let A be an algorithm which picks an integer x between 1 and N uniformly at random, for some large N, and outputs f (x) for some function f : {1, . . . , N} → {0, 1}. Then E[v(A)] = |{x : f (x) = 1}|/N. It was shown by Nayak & Wu [7] that any quantum algorithm which computes this quantity for an arbitrary function f up to (absolute or relative) error must make at most Ω(1/ ) queries to f in the case that |{x : f (x) = 1}| = N/2. As the output of A for any such function has variance 1 4 , this implies that algorithms 2 and 4 are optimal in the black-box setting in terms of their scaling with , up to polylogarithmic factors. By rescaling, we get a similar near-optimality claim for algorithm 3 in terms of its scaling with σ .

Partition function problems
In this section, we formally state and prove our results about partition function problems. We first recall the definitions from §1. A partition function Z is defined by Z(β) = x∈Ω e −β H(x) , where β is an inverse temperature and H is a Hamiltonian function taking integer values in the set {0, . . . , n}. Let 0 = β 0 < β 1 < · · · < β = ∞ be a sequence of inverse temperatures and assume that we can easily compute Z(β 0 ) = |Ω|. We want to approximate Z(∞) by approximating the ratios α i := Z(β i+1 )/Z(β i ) and using the telescoping product Finally, a sequence of Gibbs distributions π i is defined by

(a) Chebyshev cooling schedules
We start by motivating, and formally defining, the concept of a Chebyshev cooling schedule [4].
To approximate α i , we define the random variable The following result was shown by Dyer & Frieze [31] (see [4] for the statement here).
Thus, a classical algorithm can approximate Z(∞) up to relative error using O(B 2 / 2 ) samples in total, assuming that Z(0) can be computed without using any samples and that we have To characterize the latter constraint, observe that we have This motivates the following definition:
Assume that we have a sequence of estimatesα i such that, for all i, |α i − α i | ≤ ( /2 )α i with probability at least 1 − 1/(4 ). We output as a final estimateZ = Z(0)α 0α1 · · ·α −1 . By a union bound, all of the estimatesα i are accurate to within ( /2 )α i , except with probability at most 1 4 . Assuming that all the estimates are indeed accurate, we have with probability at least 3 4 . Using these ideas, we can formalize the discussion in §1. (b) Approximate sampling It is unfortunately not always possible to exactly sample from the distributions π i . However, one classical way of approximately sampling from each of these distributions is to use a (reversible, ergodic) Markov chain which has unique stationary distribution π i . Assume the Markov chain has relaxation time τ , where τ := 1/(1 − |λ 1 |), and λ 1 is the second largest eigenvalue in absolute value. Then one can sample from a distributionπ i such that π i − π i ≤ using O(τ log(1/( π min,i ))) steps of the chain, where π min,i = min x |π i (x)| [8]. We would like to replace the classical Markov chain with a quantum walk, to obtain a faster mixing time. A construction due to Szegedy [37] defines a quantum walk corresponding to any ergodic Markov chain, such that the dependence on τ in the mixing time can be improved to O( √ τ ) [12]. Unfortunately, it is not known whether in general the dependence on π min,i can be kept logarithmic [12,14]. Indeed, proving such a result is likely to be hard, as it would imply a polynomial-time quantum algorithm for graph isomorphism [13].
Nevertheless, it was shown by Wocjan & Abeyesinghe [11] (improving previous work on using quantum walks for classical annealing [23]) that one can achieve relatively efficient quantum sampling if one has access to a sequence of slowly varying Markov chains. [11]). Let M 0 , . . . , M r be classical reversible Markov chains with stationary distributions π 0 , . . . , π r such that each chain has relaxation time at most τ . Assume that | π i |π i+1 | 2 ≥ p for some p > 0 and all i ∈ {0, . . . , r − 1}, and that we can prepare the state |π 0 . Then, for any > 0, there is a quantum algorithm which produces a quantum state |π r such that |π r − |π r |0 a ≤ , for some integer a. The algorithm uses O(r √ τ log 2 (r/ )(1/p) log(1/p)) steps in total of the quantum walk operators W i corresponding to the chains M i .

Theorem 3.4 (Wocjan & Abeyesinghe
In addition, one can approximately reflect about the states |π i more efficiently still, with a runtime that does not depend on r. This will be helpful because algorithm 4 uses significantly more reflections than it does copies of the starting state. In our setting, we can easily create the quantum state |π 0 , which is the uniform superposition over all configurations x. We now show that the overlaps | π i |π i+1 | 2 are large for all i. We go via the χ 2 divergence As noted in [4], one can calculate that Therefore, if the β i values form a Chebyshev cooling schedule, χ 2 (π i+1 , π i ) ≤ B − 1 for all i. For any distributions ν, π , we also have 1 by applying Jensen's inequality to the function x → 1/ √ x. So, for all i, | π i |π i+1 | 2 ≥ 1/B. Note that in [4], it was necessary to introduce the concept of a reversible Chebyshev cooling schedule to facilitate 'warm starts' of the Markov chains used in the algorithm. That work uses the fact that one can efficiently sample from π i+1 , given access to samples from π i , if χ 2 (π i , π i+1 ) = O(1); this is the reverse of the condition (3.1). Here, we do not need to reverse the schedule as the precondition | π i |π i+1 | 2 ≥ Ω(1) required for theorem 3.4 is already symmetric.
We are now ready to formally state our result about approximating partition functions. We assume that is relatively small to simplify the bounds; this is not an essential restriction. Theorem 3.6. Let Z be a partition function. Assume we have a B-Chebyshev cooling schedule β 0 = 0 < β 1 < β 2 < · · · < β = ∞ for B = O (1) Proof. For each i, we use algorithm 4 to approximate α i up to relative error /(2 ), with failure probability γ , for some small constant γ . This would require R reflections about the state |π β i , for some R such that R = O(( / ) log 3/2 ( / ) log log( / )), and O(log( / ) log log( / )) copies of |π β i .
Instead of performing exact reflections and using exact copies of the states |π i , we use approximate reflections and approximate copies of |π i . By theorem 3.5, O( √ τ log(1/ r )) walk operations are sufficient to reflect about |π i up to an additive error term of order r . By theorem 3.4, as we have a Chebyshev cooling schedule, a quantum state |π i such that |π i − |π i |0 b ≤ s can be produced using O( √ τ log 2 ( / s )) steps of the quantum walks corresponding to the Markov chains M 0 , . . . , M i . We choose r = γ /R, s = γ . Then the final state of algorithm 4 using approximate reflections and starting with the states |π i rather than |π i can differ from the final state of an exact algorithm by at most R r + s = 2γ in 2 norm. This implies that the total variation distance between the output probability distributions of the exact and inexact algorithms is at most 2γ , and hence by a union bound that the approximation is accurate up to relative error /(2 ) except with probability 3γ . For each i, we then take the median of O(log( /δ)) estimates to achieve an estimate which is accurate up to relative error /(2 ) except with probability at most δ/ . By a union bound, all the estimates are accurate up to relative error /(2 ) except with probability at most δ, so their product is accurate to relative error except with probability at most δ.
The total number of steps needed to produce all the copies of the states |π i required is thus O · √ τ (log 2 ) · log log log · log δ and the total number of steps needed to perform the reflections is O( · √ τ (log R) · R · log( /δ)).
Adding the two, substituting the value of R, and using We remark that, in the above complexities, we have chosen to take the number of quantum walk steps used as our measure of complexity. This is to enable a straightforward comparison with the classical literature, which typically uses a random walk step as its elementary operation for the purposes of measuring complexity [4]. To implement each quantum walk step efficiently and accurately, two possible approaches are to use efficient state preparation [38] or recently developed approaches to efficient simulation of sparse Hamiltonians [39].

(c) Computing a Chebyshev cooling schedule
We still need to show that, given a particular partition function, we can actually find a Chebyshev cooling schedule. For this, we simply use a known classical result: We remark that a subsequent algorithm [40] improves the polylogarithmic terms and the hidden constant factors in the complexity. However, this algorithm assumes that we can efficiently generate independent samples from distributions approximating π β for arbitrary β. The most efficient general algorithm known [4] for approximately sampling from arbitrary distributions π β uses 'warm starts' and hence does not produce independent samples.
Combining all the ingredients, we have the following result. The best comparable classical result known isÕ((log A)τ/ 2 ) [4]. We therefore see that we have achieved a near-quadratic reduction in the complexity with respect to both τ and , assuming that ≤ 1/ √ τ . Otherwise, we still achieve a near-quadratic reduction with respect to .

(d) Some partition function problems
In this section, we describe some representative applications of our results to problems in statistical physics and computer science.

(i) The ferromagnetic Ising model
This well-studied statistical physics model is defined in terms of a graph G = (V, E) by the Hamiltonian H(z) = − (u,v)∈E z u z v , where |V| = n and z ∈ {±1} n . A standard method to approximate the partition function of the Ising model uses the Glauber dynamics. This is a simple Markov chain with state space {±1} n , each of whose transitions involves only updating individual sites, and whose stationary distribution is the Gibbs distribution π β (z) = (1/Z(β)) e −βH(z) . This Markov chain, which has been intensively studied for decades, is known to mix rapidly in certain regimes [41]. Here, we mention just one representative recent result.