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 . The natural algorithm for doing so is to produce k samples, each corresponding to the output of an independent execution of , 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 is at most σ2, the probability that the value output by this estimator is far from the truth can be bounded using Chebyshev's inequality:

It is therefore sufficient to take k=O(σ2/ϵ2) to estimate μ up to additive error ϵ with, say, 99% success probability. This simple result is a key component in many more complex randomized approximation schemes (e.g. [4,1]).

Although this algorithm is fairly efficient, its quadratic dependence on σ/ϵ seems far from ideal: for example, if σ=1, to estimate μ up to four decimal places, we would need to run over 100 million times. Unfortunately, it can be shown that, without any further information about , the sample complexity of this algorithm is asymptotically optimal [5] with respect to its scaling with σ and ϵ, although it can be improved by a constant factor [6].

We show here that, using a quantum computer, the number of uses of required to approximate μ can be reduced almost quadratically beyond the above classical bound. Assuming that the variance of the output of the algorithm is at most σ2, we present a quantum algorithm which estimates μ up to additive error ϵ, with 99% success probability, using only times.[1] It follows from known lower bounds on the quantum complexity of approximating the mean [7] that the runtime of this algorithm is optimal, up to polylogarithmic factors. This result holds for an arbitrary algorithm used as a black box, given only an upper bound on the variance.

An important aspect of this construction is that the underlying subroutine need not be a classical randomized procedure, but can itself be a quantum algorithm. This enables any quantum speedup obtained by to be used within the overall framework of the algorithm. A particular case in which this is useful is quantum speedup of Markov chain Monte Carlo methods [8]. Classically, such methods use a rapidly mixing Markov chain to approximately sample from a probability distribution corresponding to the stationary distribution of the chain. Quantum walks are the quantum analogue of random walks (e.g. [9] for a review). In some cases, quantum walks can reduce the mixing time quadratically (e.g. [10,11]), although it is not known whether this can be achieved in general [12,13,14]. We demonstrate that this known quadratic reduction can be combined with our algorithm to speed up the fastest known general-purpose classical algorithm with rigorous performance bounds [4] for approximately computing partition functions up to small relative error, a fundamental problem in statistical physics [1]. As another example of how our algorithm can be applied, we substantially improve the runtime of a quantum algorithm for estimating the total variation distance between two probability distributions [15].

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 L2 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−. 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− 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.

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 , we write for the random variable corresponding to the value computed by , with the expected value of denoted . For concreteness, we think of 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}, the final output is then ϕ(x), for some fixed function . If 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 reversible-computation 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 , which is the inverse of the unitary part of . If we do have a description of as a quantum circuit, this can be achieved simply by running the circuit backwards, replacing each gate with its inverse.

Table 1.

Summary of the main quantum algorithms presented in this paper for estimating the mean output value μ of an algorithm . (Algorithm 2, omitted, is a subroutine used in algorithm 3.)

algorithm	precondition	approximation of μ	uses of A and A^−1
1	v(A)∈[0,1]	additive error ϵ	O(1/ϵ)
3	Var(v(A))≤σ^2	additive error ϵ	O~(σ/ϵ)
4	Var(v(A))/(E[v(A)])^2≤B	relative error ϵ	O~(B/ϵ)

We first deal with the special case where the output of is bounded between 0 and 1. Here, a quantum algorithm for approximating 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 . If the final measurement of the algorithm 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 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 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 . 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 up into disjoint intervals depending on size. Write for the ‘truncated’ algorithm which outputs if , and otherwise outputs 0. We estimate μ by applying the above algorithm to estimate for a sequence of intervals which are exponentially increasing in size, and summing the results. As the intervals [p,q) get larger, the accuracy with which we approximate decreases, and values larger than about 1/ϵ are ignored completely. However, the overall upper bound on 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 would be high.

The final result is that for , given appropriate parameter choices, the estimate satisfies with high probability, and the algorithm uses 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 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 to make its ℓ2 norm low. If has mean μ and variance upper-bounded by σ2, a suitable way to achieve this is to subtract μ from the output of , 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 once and use the output as a proxy for μ. Because , is quite likely to be within distance O(σ) of μ. Therefore, the algorithm produced from by subtracting 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 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 . In this setting, it turns out that uses of 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.

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 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/(kT), where T is the temperature and k 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, 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 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 output such that , 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 be a sequence of inverse temperatures. A standard classical approach to design algorithms for approximating partition functions [30,31,32,4,18] is based around expressing Z(βℓ) as the telescoping product If we can compute Z(β0)=|Ω| and can also approximate each of the ratios α:=Z(β)/Z(β) accurately, taking the product will give a good approximation to Z(βℓ). Let π denote the Gibbs (or Boltzmann) probability distribution corresponding to inverse temperature β, where To approximate α, we define the random variable Then one can readily compute that , so sampling from each distribution π allows us to estimate the quantities α. It will be possible to estimate α up to small relative error efficiently if the ratio is low. This motivates the concept of a Chebyshev cooling schedule [4]: a sequence of inverse temperatures β such that for all i. It is known that, for any partition function problem as defined above such that |Ω|=A, there exists a Chebyshev cooling schedule with [4].

It is sufficient to approximate 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 . Given a Chebyshev cooling schedule with , the algorithm thus uses 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 π. Classically, one way of approximately sampling from these distributions is to use a Markov chain which mixes rapidly and has stationary distribution π. 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 π 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 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 |π〉 using quantum walk steps each. Using similar ideas, we can approximately reflect about |π〉 using only quantum walk steps.

Approximating up to relative error O(ϵ/ℓ) using our algorithm requires one quantum sample approximating |π〉, and approximate reflections about |π〉. 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 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 . 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.

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 G=(V,E) by the Hamiltonian , where |V |=n and z∈{±1}. 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 [33]. In this case, as A=2, the quantum algorithm approximates Z(β) to within relative error ϵ in steps. The corresponding classical algorithm [4] uses 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 time , as compared with the classical [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 steps, as compared with the classical [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 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 R(x)=(|p(x)−q(x)|)/(p(x)+q(x)), where x is drawn from the distribution r=(p+q)/2 [15]. For each x, R(x) can be computed up to accuracy ϵ using iterations of amplitude estimation. Wrapping this within O(1/ϵ) iterations of the mean-estimation algorithm, we obtain an overall algorithm running in time . See §4 for details.

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 , we let be the random variable corresponding to the value computed by . We assume that 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 is a quantum algorithm of the following form: apply some unitary operator to the initial state |0〉; measure k≤n qubits of the resulting state in the computational basis, obtaining outcome x∈{0,1}; output ϕ(x) for some easily computable function . We finally assume that we have access to the inverse of the unitary part of the algorithm, which we write as .

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 be a (classical or quantum) algorithm which aims to estimate some quantity μ, and whose output satisfies except with probability γ, for some fixed . Then, for any δ>0, it suffices to repeat 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 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 multiplicative factor.

Estimating the mean with bounded output values

We first consider the problem of estimating in the special case where 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].

Theorem 2.3

Let |ψ〉 be defined as in algorithm 1 and set U=2|ψ〉〈ψ|−I. Algorithm 1 uses copies of the state uses U times and outputs an estimate such that 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 1−δ, it suffices to take . To achieve with probability at least 1−δ, it suffices to take t=O(1/ϵ).

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 , just before its last measurement, can be written as for some normalized states |ψ〉. If we then attach an ancilla qubit and apply W, we obtain We have where P=I⊗|1〉〈1|. Therefore, when we apply amplitude estimation, by theorem 2.2, we obtain an estimate of such that with probability at least 8/π2. The powering lemma (lemma 2.1) implies that the median of 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 and , 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.

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 any algorithm , let , , , be the algorithms defined by executing to produce a value and:

— : If , output , otherwise output 0;

— : If , output , otherwise output 0;

— : If , output , otherwise output 0.

In addition, for any algorithm and any function , let be the algorithm produced by evaluating and computing . Note that algorithm 1 can easily be modified to compute rather than , for any function , by modifying the operation W.

Our algorithm (algorithm 2) and correctness proof are a generalization of a result of Heinrich [16] for computing the mean with respect to the uniform distribution of functions with bounded L2 norm, and are based on the same ideas. Write .

Lemma 2.4

Let U=2|ψ〉〈ψ|−I. Algorithm 2 uses copies of |ψ〉, uses U times and estimates up to additive error with probability at least .

We first show the resource bounds. Algorithm 1 is run times, each time with parameter . By theorem 2.3, each use of algorithm 1 consumes copies of |ψ〉 and uses U times. The total number of copies of |ψ〉 used is , and the total number of uses of U is .

All of the uses of algorithm 1 succeed, except with probability at most 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 outputs x. We have By theorem 2.3, 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 , we obtain Inserting the definitions of t0 and k, we get an overall error bound using in the second inequality. For a sufficiently large constant D, this is upper-bounded by as claimed. ▪

Observe that, if , to achieve additive error ϵ the number of uses of that we need is . By the powering lemma, we can repeat algorithm 2 times and take the median to improve the probability of success to 1−δ for any δ>0.

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 up to additive error ϵ for an arbitrary algorithm such that , we need to use Ω(σ2/ϵ2) times [5].

Theorem 2.5

Let U=2|ψ〉〈ψ|−I. Algorithm 3 uses copies of |ψ〉, uses U times and estimates up to additive error ϵ with success probability at least .

First, observe that is quite close to with quite high probability. As , by Chebyshev's inequality, we have . We therefore assume that . In this case, we have where the first inequality is the triangle inequality. Thus , which implies that and .

The next step is to use algorithm 2 to estimate and with accuracy ϵ/(32σ) and failure probability . By lemma 2.4, if the algorithm succeeds in both cases, the estimates are accurate up to ϵ/(8σ). We therefore obtain an approximation of each of and up to additive error ϵ/(2σ). As we have by linearity of expectation, using a union bound we have that approximates up to additive error ϵ with probability at least . ▪

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 , where we normally think of B as small, e.g. B=O(1). Classically, it follows from Chebyshev's inequality that the simple classical algorithm described in the Introduction approximates up to additive error with O(B/ϵ2) uses of . In the quantum setting, we can improve the dependence on ϵ near-quadratically; we describe this as algorithm 4 below.

Theorem 2.6

Let U=2|ψ〉〈ψ|−I. Algorithm 4 uses copies of |ψ〉, uses U times and outputs an estimate such that .

The complexity bounds follow from lemma 2.4; we now analyse the claim about accuracy. is a random variable whose expectation is and whose variance is . By Chebyshev's inequality, we have We can thus assume that . In this case, when we apply algorithm 2 to , we receive an estimate of which is accurate up to additive error except with probability , where we use . Multiplying by and taking a union bound, we get an estimate of which is accurate up to ϵ except with probability at most . ▪

Once again, using the powering lemma, we can repeat algorithms 3 and 4 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 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 . 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 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 . Then . 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 for any such function has variance , 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 , where β is an inverse temperature and H is a Hamiltonian function taking integer values in the set {0,…,n}. Let be a sequence of inverse temperatures and assume that we can easily compute Z(β0)=|Ω|. We want to approximate by approximating the ratios α:=Z(β)/Z(β) and using the telescoping product Finally, a sequence of Gibbs distributions π is defined by π(x)=(1/Z(β)) e−.

Chebyshev cooling schedules

We start by motivating, and formally defining, the concept of a Chebyshev cooling schedule [4]. To approximate α, we define the random variable Y (x)=e−(. Then The following result was shown by Dyer & Frieze [31] (see [4] for the statement here).

Theorem 3.1

Let Y 0,…,Y ℓ−1 be independent random variables such that for all i, and write . Let be the average of 16Bℓ/ϵ2 independent samples from Y , and set . Then .

Thus, a classical algorithm can approximate 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 so This motivates the following definition:

Definition 3.2 (Chebyshev cooling schedules [4])

Let Z be a partition function. Let β0,…,βℓ be a sequence of inverse temperatures such that . The sequence is called a B-Chebyshev cooling schedule for Z if for all i, for some fixed B.

Assume that we have a sequence of estimates such that, for all i, with probability at least 1−1/(4ℓ). We output as a final estimate . By a union bound, all of the estimates are accurate to within (ϵ/2ℓ)α, except with probability at most . Assuming that all the estimates are indeed accurate, we have for ϵ<1. Thus, with probability at least .

Using these ideas, we can formalize the discussion in §1.

Theorem 3.3

Let Z be a partition function with |Ω|=A. Assume that we are given a B-Chebyshev cooling schedule for Z. Further assume that we have the ability to exactly sample from the distributions π, i=1,…,ℓ−1. Then there is a quantum algorithm which outputs an estimate such that using samples in total.

For each i=1,…,ℓ−1, we use algorithm 4 to estimate up to additive error with failure probability 1/(4ℓ). As the β form a B-Chebyshev cooling schedule, , so . By theorem 2.6, each use of algorithm 4 requires samples from π to achieve the desired accuracy and failure probability. The total number of samples is thus as claimed. ▪

Approximate sampling

It is unfortunately not always possible to exactly sample from the distributions π. 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 π. 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 such that using steps of the chain, where [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 [12]. Unfortunately, it is not known whether in general the dependence on 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.

Theorem 3.4 (Wocjan & Abeyesinghe [11])

Let M0,…,M be classical reversible Markov chains with stationary distributions π0,…,π such that each chain has relaxation time at most τ. Assume that |〈π|π〉|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 such that for some integer a. The algorithm uses steps in total of the quantum walk operators W corresponding to the chains M.

In addition, one can approximately reflect about the states |π〉 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.

Theorem 3.5 (Wocjan & Abeyesinghe [11], see [18] for version here)

Let M0,…,M be classical reversible Markov chains with stationary distributions π0,…,π such that each chain has relaxation time at most τ. For each i, there is an approximate reflection operator such that where |ϕ〉 is arbitrary, and |ξ〉 is a vector with ∥|ξ〉∥≤ϵ. The algorithm uses steps of the quantum walk operator W corresponding to the chain M.

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 |〈π|π〉|2 are large for all i. We go via the χ2 divergence As noted in [4], one can calculate that Therefore, if the β values form a Chebyshev cooling schedule, χ2(π,π)≤B−1 for all i. For any distributions ν, π, we also have by applying Jensen's inequality to the function . So, for all i, |〈π|π〉|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 π, given access to samples from π, if χ2(π,π)=O(1); this is the reverse of the condition (3.1). Here, we do not need to reverse the schedule as the precondition |〈π|π〉|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 for B=O(1). Assume that for every inverse temperature β we have a reversible ergodic Markov chain M with stationary distribution π and relaxation time upper-bounded by τ. Further assume that we can sample directly from M0. Then, for any δ>0 and there is a quantum algorithm which uses steps of the quantum walks corresponding to the M chains and outputs such that .

For each i, we use algorithm 4 to approximate α up to relative error ϵ/(2ℓ), with failure probability γ, for some small constant γ. This would require R reflections about the state |π〉, for some R such that , and copies of |π〉.

Instead of performing exact reflections and using exact copies of the states |π〉, we use approximate reflections and approximate copies of |π〉. By theorem 3.5, walk operations are sufficient to reflect about |π〉 up to an additive error term of order ϵ. By theorem 3.4, as we have a Chebyshev cooling schedule, a quantum state such that can be produced using steps of the quantum walks corresponding to the Markov chains M0,…,M.

We choose ϵ=γ/R, ϵ=γ. Then the final state of algorithm 4 using approximate reflections and starting with the states rather than |π〉 can differ from the final state of an exact algorithm by at most Rϵ+ϵ=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 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 required is thus and the total number of steps needed to perform the reflections is . Adding the two, substituting the value of R, and using , we get an overall bound of as claimed. ▪

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].

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:

Theorem 3.7 (Štefankovič et al. [4])

Let Z be a partition function. Assume that for every inverse temperature β we have a Markov chain M with stationary distribution π and relaxation time upper-bounded by τ. Further assume that we can sample directly from M0. Then, for any δ>0 and any B=O(1), we can produce a B-Chebyshev cooling schedule of length with probability at least 1−δ, using at most steps of the Markov chains.

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.

Corollary 3.8

Let Z be a partition function and let ϵ>0 be a desired precision such that . Assume that for every inverse temperature β, we have a Markov chain M with stationary distribution π and relaxation time upper-bounded by τ. Further assume that we can sample directly from M0. Then, for any δ>0, there is a quantum algorithm which uses steps of the M chains and their corresponding quantum walk operations, and outputs such that .

The best comparable classical result known is [4]. We therefore see that we have achieved a near-quadratic reduction in the complexity with respect to both τ and ϵ, assuming that . Otherwise, we still achieve a near-quadratic reduction with respect to ϵ.

Some partition function problems

In this section, we describe some representative applications of our results to problems in statistical physics and computer science.

The ferromagnetic Ising model

This well-studied statistical physics model is defined in terms of a graph G=(V,E) by the Hamiltonian , where |V |=n and z∈{±1}. 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}, each of whose transitions involves only updating individual sites, and whose stationary distribution is the Gibbs distribution π(z)=(1/Z(β)) e−. 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.

Theorem 3.9 (Mossel & Sly [33])

For any integer d>2, and inverse temperature β>0 such that the mixing time of the Glauber dynamics on any graph of maximum degree d is .

(More precise results than theorem 3.9 are known for certain specific graphs such as lattices [42].) As we have A=2, in the regime where , the quantum algorithm approximates Z(β) to within ϵ relative error in steps. The fastest known classical algorithm with rigorously proved performance bounds [4] uses time . We remark that an alternative approach of Jerrum & Sinclair [29], which is based on analysing a different Markov chain, gives a polynomial-time classical algorithm which works for any temperature, but is substantially slower.

Counting colourings

Here, we are given as input 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, and k=O(1). In physics, this problem corresponds to the partition function of the Potts model evaluated at zero temperature. It is known that the Glauber dynamics for the Potts model mixes rapidly in some cases [43]. One particularly clean result of this form is work of Jerrum [44] showing that this Markov chain mixes in time if k>2d. As here A=k, we obtain a quantum algorithm approximating the number of colourings of G up to relative error ϵ in steps, as compared with the classical [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 often known as monomer–dimer coverings [34]. To count the number of matchings, we consider the partition function , where is the set of matchings of G. We have , while , as in this case the sum is zero everywhere except the empty matching (00=1). Therefore, in this case, we seek to approximate Z(0) using a telescoping product which starts with . In terms of the cooling schedule , we have As we have reversed our usage of the cooling schedule, rather than looking for it to be a B-Chebyshev cooling schedule, we instead seek the bound Z(2β−β)Z(β)/Z(β)2≤B to hold for all i=0,…,ℓ−1. That is, the roles of β and β have been reversed as compared with definition 3.2. However, the classical algorithm for printing a cooling schedule can be modified to output a ‘reversible’ schedule where this constraint is satisfied too, with only a logarithmic increase in complexity [4]. In addition, it was shown by Jerrum & Sinclair [45,46] that, for any β, there is a simple Markov chain which has stationary distribution π, where and which has relaxation time τ=O(nm) on a graph with n vertices and m edges. Finally, in the setting of matchings, A=O(n!2). Putting these parameters together, we obtain a quantum complexity , as compared with the lowest known classical bound [4].

Estimating the total variation distance

Here, we give the technical details of our improvement of the accuracy of a quantum algorithm of Bravyi et al. [15] 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 up to additive error ϵ. Classically, estimating ∥p−q∥ up to error, say, 0.01 cannot be achieved using O(n) samples for any α<1 [35], but in the quantum setting the dependence on n can be improved quadratically:

Theorem 4.1 (Bravyi et al. [15])

Given the ability to sample from p and q, there is a quantum algorithm which estimates ∥p−q∥ up to additive error ϵ, with probability of success 1−δ, using samples.

Here, we will use theorem 2.3 to improve the dependence on ϵ and δ of this algorithm. We will approximate the mean output value of a subroutine previously used in [15] (algorithm 5).

If the estimates , in this subroutine were precisely accurate, the expected output of the subroutine would be We now bound how far the expected output of the algorithm is from this exact value. By linearity of expectation, where d(x)=|p(x)−q(x)|/(p(x)+q(x)) and . Note that is a random variable. Split [n] into ‘small’ and ‘large’ parts according to whether r(x)≤ϵ/n. Then using that . From theorem 2.2, for any δ>0, we have except with probability at most δ, using samples from p. If for some 0≤η≤1, this implies that except with probability at most δ. A similar claim also holds for . We now use the following technical result from [15]:

Proposition 4.2

Consider a real-valued function f(p,q)=(p−q)/(p+q), where 0≤p,q≤1. Assume that for some . Then .

By proposition 4.2, for all x such that , we have , except with probability at most 2δ. We now fix . Then, for all x such that p(x)+q(x)≥2ϵ/n, except with probability at most 2δ. Thus, for all x such that r(x)≥ϵ/n, . Taking δ=ϵ, we have for any ϵ, using samples. It therefore suffices to use samples to achieve . As the output of this subroutine is bounded between 0 and 1, to approximate up to additive error ϵ/2 with failure probability δ, it suffices to use the subroutine times by theorem 2.3. So the overall complexity is . For small ϵ and δ, this is a substantial improvement on the complexity stated in [15].