A Sublinear-Time Quantum Algorithm for Approximating Partition Functions

We present a novel quantum algorithm for estimating Gibbs partition functions in sublinear time with respect to the logarithm of the size of the state space. This is the first speed-up of this type to be obtained over the seminal nearly-linear time algorithm of \v{S}tefankovi\v{c}, Vempala and Vigoda [JACM, 2009]. Our result also preserves the quadratic speed-up in precision and spectral gap achieved in previous work by exploiting the properties of quantum Markov chains. As an application, we obtain new polynomial improvements over the best-known algorithms for computing the partition function of the Ising model, and counting the number of $k$-colorings, matchings or independent sets of a graph. Our approach relies on developing new variants of the quantum phase and amplitude estimation algorithms that return nearly unbiased estimates with low variance and without destroying their initial quantum state. We extend these subroutines into a nearly unbiased quantum mean estimator that reduces the variance quadratically faster than the classical empirical mean. No such estimator was known to exist prior to our work. These properties, which are of general interest, lead to better convergence guarantees within the paradigm of simulated annealing for computing partition functions.


Introduction
The Boltzmann-Gibbs distribution is a paradigmatic tool for modeling systems that obey the principle of maximum entropy. It arises in several fields of research such as statistical mechanics [Geo11; FV17; Sin82], economic modeling [DY00], image processing [GG84;Bes86], statistical learning theory [Cat04], etc. The probability assigned by the Gibbs distribution to each possible configuration of a system is inversely proportional to the exponential of its energy multiplied by the inverse temperature. Mathematically, for a classical Hamiltonian H : Ω → {0, . . . , n} specifying the energy level of each configuration x, the Gibbs distribution at inverse temperature β is given by π β (x) = 1 Z(β) e −βH(x) where the normalization factor is called the partition function. While it is often straightforward to evaluate the partition function at high temperature (when β = 0, it is just the number of possible configurations), the low-temperature regime captures ground state properties that are challenging to compute. For instance, Z(∞) can represent the cardinalities of exponentially large combinatorial structures (such as the number of colorings of a graph [Jer95;ŠVV09], the volume of convex bodies [DFK91;DF91] or the permanent of non-negative matrices [JSV04]) which are generally #P-hard to compute exactly [Val79; DF88; JS93].
The standard approach for evaluating partition functions at low-temperature is to resort to Markov chain Monte Carlo methods [JS96]. A celebrated line of works [VC72; JVV86; DF91; BŠVV08; ŠVV09] has shown how to turn the ability to efficiently sample from the Gibbs distribution into that to efficiently approximate the partition function. At a high level, these works rely on the same two-stage simulated annealing algorithm. First, they compute a short cooling schedule, which is an increasing sequence of inverse temperature 0 = β 1 < · · · < β = ∞ with limited fluctuations in Gibbs distributions between two consecutive values. Next, the partition function is expressed as a telescoping product which is approximated by using a suitable product estimator. As an example, the seminal algorithm of Štefankovič, Vempala and Vigoda [ŠVV09] generates a so-called "Chebyshev cooling schedule" of length 1 = O( log|Ω|) (ignoring logarithmic dependences on the degree n, which is often on the order of n ∼ log|Ω|) where each ratio Z(β i+1 )/Z(β i ) is expressed as the expectation value of a random variable X i with bounded relative second moment E[X 2 i ] ≤ O(E[X i ] 2 ). Such schedules are known to admit a product estimator that requires O( 2 / 2 ) classical Gibbs samples to estimate Z(∞) = Z(0) · E[X 1 ] · · · E[X ] with relative error . Thus, if we let δ denote the spectral gap of a Markov chain generating samples from the considered Gibbs distributions, the overall cost of this classical algorithm is O(log|Ω|/( 2 δ)).
The theory of quantum algorithms provides several directions for accelerating the computation of partition functions. Quantum Markov chains [Sze04;MNRS11] can prepare coherent "qsample" encodings |π β = x π β (x)|x of the Gibbs distribution with a quadratic improvement in spectral gap for the rate of convergence (but an increase dependence on other parameters). Quantum phase estimation [Kit95] and amplitude estimation [BHMT02] lead to quadratically better convergence rates for estimating expectation values [Ter99; AW99; Hei02; BDGT11; Mon15; HM19;Ham21]. Yet, while this may hint at the existence of an O( log|Ω|/( √ δ)) quantum algorithm for estimating partition functions, the best known algorithms [HW20; AHN+21] still require a linear scaling O(log|Ω|/( √ δ)) with the logarithm of the size of the state space. This bottleneck is due to additional challenges posed by current quantum algorithmic techniques. It is for instance significantly harder to prepare the qsample |π β (at low temperature) than to implement the reflection I − 2|π β π β | through it. This obstacle requires using nondestructive procedures [MW05; WA08; WCNA09; TOV+11; ORR13; HW20] to recycle the same qsamples all along the algorithm, and to rely mostly on the reflection operator. Another fundamental limitation faced by current best quantum mean estimators [Mon15; HM19;Ham21] is the presence of biases in the estimates that degrade the convergence guarantee of the product estimators.

Contributions
Our main contribution is to develop the first quantum algorithm for approximating Gibbs partition functions with a complexity scaling sublinearly with respect to the logarithm of the size of the state space. Formally, we prove the next theorem in Section 4.
Theorem 4.1 (Informal). There is a quantum algorithm such that, given a Gibbs distribution generated by a Markov chain with spectral gap δ, it computes an estimate Z of the partition function at zero temperature satisfying | Z − Z(∞)| ≤ Z(∞) by using O log 3/4 |Ω| log 3/4 (n)/ + log |Ω| log(n) steps of the quantum walk operator.
Our result achieves state-of-the-art dependence on the other parameters (spectral gap δ, accuracy , degree n of the partition function) up to logarithmic factors. We provide a comparison with prior work in Table 1.
The main technical ingredients of our work are new variants of the quantum phase and amplitude estimation algorithms (Theorems 2.2 and 2.4) that return (nearly) unbiased estimates and restore the qsamples used in the process with high probability. We extend these properties into new quantum mean estimators described in Section 3. In particular, the next result is crucial in the development of our improved product estimator (Theorem 3.3), and may be of independent interest. It shows that the relative error in bias can be decreased exponentially quickly, whereas previous works [Mon15; HM19; Ham21] incurred an inverse linear dependence on .
Proposition 3.1 (Informal). For any parameters t, , σ there is a quantum algorithm such that, given a random variable X with variance at most σ, it computes a mean estimate µ with bias |E[ µ] − µ| ≤ σ and variance Var[ µ] ≤ σ t 2 by using one copy of a qsample |π X (restored at the end of the computation) and O(t log(1/ ) 2 ) applications of the reflection I − 2|π X π X |.
Finally, we provide applications of our work in Section 4.2, where we improved upon the best known algorithms for several approximate counting problems that can be phrased as partition functions. These results are described succinctly in Table 2. Colorings, Ising model, Matchings Independent sets [ŠVV09]

Proof overview
We give a high-level description of the algorithms developed in this paper. We first explain the main differences with previous works [WCNA09; Mon15; HW20; AHN+21] that allow us to obtain a sublinear algorithm for estimating partition functions. Next, we present new variants of phase estimation, amplitude estimation and quantum mean estimation that are needed for implementing our approach.
Estimating partition functions (Section 4). The current bottleneck in quantizations of the classical simulated annealing algorithm of Štefankovič, Vempala and Vigoda [ŠVV09] lies in estimating the expectation of a Chebyshev cooling schedule for the telescoping product displayed in Equation (2). Abstractly, this problem amounts to estimating with relative error the expectation of the product X = X 1 · · · X of independent random variables with bounded relative variance . By Bernoulli's inequality, one can show that estimating the expectation of each inner random variable X i with relative error O( / ) and taking the product of the estimates is sufficient for that purpose. Montanaro [Mon15] developed a quantum algorithm for computing each such inner estimate in time O( / ), quadratically faster than it is possible classically, leading to an overall complexity of O( 2 / ). Nevertheless, a different classical estimator from Dyer and Frieze [DF91] provides the same quadratic complexity with respect to the schedule length: replace the / -error inner estimates with unbiased estimates of variance . The latter estimates can simply be obtained by averaging O( / 2 ) samples due to the bounded variance property. The overall variance of this new product estimator is shown to be O(( E[X]) 2 ) which, by Chebyshev's inequality and the fact that it is unbiased, implies that the overall relative error is . We follow a similar approach in our work to speed up the estimation of E[X]. Our main contribution is a new quantum estimator that reduces the variance quadratically faster than classically, while keeping the estimates (nearly) unbiased. This allows us to compute each inner estimate in time O( √ / ), leading to an overall complexity of O( 3/2 / ). Our estimator requires new unbiased variants of the phase and amplitude estimation algorithms, which are sketched in the next paragraphs. Additionally, these algorithms are made nondestructive (i.e. they restore a copy of their starting state) to ensure the reusability of the qsamples which are an expensive resource in our applications.
Unbiased phase estimation (Section 2.1). The phase estimation problem asks the question of estimating the eigenphase θ ∈ [0, 1) of a quantum state |ψ satisfying U |ψ = e 2πiθ |ψ for a given unitary operator U . Kitaev's algorithm [Kit95] returns an estimate θ within distance 1/t from θ with constant success probability by using O(t) controlled-U operations. By standard success amplification techniques, the variance of this estimate can be guaranteed to be of the order of 1/t 2 + at the cost of an O(log 1/ ) overhead in complexity. One simple idea 2 to make this estimator unbiased is to consider the phase-shifted unitary U ϕ |ψ = e 2πi(θ+ϕ) |ψ for a random ϕ, run phase estimation on it and subtract the shift ϕ from the estimate. Although it provides an unbiased estimate of e 2πiθ (and, therefore, of sin 2 (θ)), we do not know how to reduce its variance without introducing a new bias. We overcome this difficulty by exploiting more properties of the output distribution of phase estimation and looking for an unbiased estimate of θ directly. First, we show that it is possible to compute with high probability the exact value of the first log t bits of θ + ϕ = 0.b 1 b 2 . . . for a well-chosen phase shift ϕ (steps 1-5 in Algorithm 1). This relies on the ability to pinpoint exactly the interval k t , k+1 t for k ∈ {0, . . . , t − 1} containing the phase θ + ϕ, as long as the latter is not too close to an interval endpoint (which is avoided by the choice of ϕ). Secondly, we look for a low-variance unbiased estimate of the remaining bits 0.0 log t b log t+1 b log t+2 . . . (steps 6-7 in Algorithm 1). For that, we consider the exponentiated unitary U t ϕ which amplifies the eigenphase to λ = 0.b log t+1 b log t+2 . . . . Using phase-to-amplitude conversion techniques, we transform it into a new unitary V that encodes λ into the amplitude of a new auxiliary quantum register. This register can simply be measured to obtain an unbiased binary estimate λ ∈ {−1, 1} of λ. The variance of the correction term λ/t added to the previously computed bits 0.b 1 b 2 . . . b log t is thus of order O(1/t 2 ). We also have to be careful that λ is bounded away from 0 and 1 in order for the conversion techniques to succeed. This is again guaranteed by carefully choosing the phase shift. The final algorithm is summarized in Theorem 2.2.
Nondestructive amplitude estimation (Appendix B). The amplitude estimation algorithm [BHMT02] approximates the measurement probability p = Π|ψ 2 for a state |ψ and a projector Π, given access to the reflection operators R ψ = I − 2|ψ ψ| and R Π = I − 2Π. This algorithm relies on the observation that the Grover operator G = −R ψ R Π has two conjugate eigenvalues e ±2i arcsin( √ p) with |ψ being a superposition over the two corresponding eigenvectors. In its original form [BHMT02], the algorithm consists of running phase estimation on the latter state and applying the function x → sin 2 (πx) on the measured register. Although the post-measurement state is still a superposition over the two above eigenvectors, this process can introduce a relative phase that makes it far from the initial state |ψ . Harrow and Wei [HW20] argued that amplitude estimation can be made nondestructive (i.e. it restores the state |ψ ) by using an uncomputation trick inspired by the Marriott-Watrous rewinding technique [MW05]. We discuss this result in Appendix B, where we provide a new proof of it that avoids some complications in the analysis.
Nondestructive unbiased amplitude estimation (Section 2.2). We aim at combining the above two properties in order to obtain a variant of amplitude estimation that is both nondestructive and unbiased. Unfortunately, the properties of the Grover operator do not allow us to simply use the unbiased phase estimation on it. Indeed, it would require the ability to shift its two eigenphases by opposite values simultaneously. Moreover, it would provide an unbiased estimate of arcsin( √ p) instead of p. We overcome these two obstacles by first applying the function x → sin 2 (πx) on the eigenphases of the Grover operator (via amplitudeto-phase conversion techniques [GAW17]), which produces a new unitary G with "merged" eigenvalue e 2πip and eigenvector |ψ . We next apply the unbiased phase estimation on G and |ψ , which is inherently nondestructive since the latter state is now an eigenvector. Yet, this unbiased estimator poses a new problem: its variance is worse than that of the original amplitude estimation algorithm (roughly, 1/t 2 vs. p/t 2 + 1/t 4 ), which is due to encoding p instead of arcsin( √ p) into the phase. Our solution to this problem is to first compute a rough biased amplitude estimate q = Θ(max(p, 1/t 2 )) (via standard nondestructive amplitude estimation) and to apply the new function x → sin 2 (πx)/q on the eigenphases of the Grover operator (via linear amplitude amplification [Low17]). Since the obtained operator will require time O(1/ √ q) to be simulated, we are left with O( √ qt) steps available to perform unbiased phase estimation on it.
The variance is then shown to be of the order of p/t 2 when p = O(1/t 2 ) and q/t 2 otherwise. A side tool of independent interest (needed when p is smaller than O(1/t 2 )) is a nondestructive coin flipping algorithm (Proposition 2.3) for sampling from the Bernoulli variable of parameter Π|ψ 2 .
Nondestructive unbiased mean estimator (Section 3). Our final step is to generalize the amplitude estimation algorithm to estimate the expectation of an arbitrary random variable X whose distribution is encoded into a qsample |π X = x Pr[X = x]|x (amplitude estimation handles the case of Bernoulli distributions). The estimator shall meet the three requirements of being unbiased, nondestructive and low-variance in order to be applied efficiently to the partition function problems. We build upon the quantum estimators from [Hei02; Mon15; HM19; Ham21] that already satisfy the low-variance property. In the latter works, the estimation of E[X] is reduced to estimating the means of a series of Bernoulli distributions whose parameters are determined by certain quantiles of the input distribution. We follow the same approach, by instead using the nondestructive unbiased amplitude estimation algorithm in the reduction. This is not enough yet to make the estimator nondestructive and unbiased. First, we modify a recentering step in [Mon15; Ham21] that computed a rough mean estimate by averaging classical samples (obtained from destructive measurements of copies of |π X ). We develop a new median estimator for this purpose that needs not destroy any copy of |π X (Proposition 3.2). Secondly, the existing estimators require truncating the random variable X by replacing its outcomes larger than a certain threshold value with zero. We suppress most of the truncation-induced bias by adding extra Bernoulli distributions to the reduction. The sum of these added estimates has low variance and is equal (in expectation) to most of the previously truncated tail of X. The final estimator is summarized in Proposition 3.1.

Phase estimation
In this section, we describe a variant of phase estimation that returns a (nearly) unbiased estimate. We use the next notation for separating the most and least significant bits of a phase θ ∈ [0, 1].
Suppose that we can call a controlled unitary U , and that we have a copy of an eigenstate |ψ , satisfying U |ψ = e 2πiθ |ψ . Then, we can implement the phase-shift operation e 2πiϕ U |ψ = e 2πi(θ+ϕ) |ψ for any ϕ with one call to U , and the bit-shift operation U 2 τ |ψ = e 2πiθ>τ |ψ with 2 τ calls to U . We make use of these two operations to design the following unbiased, nondestructive, low-variance phase estimation algorithm.

4:
if there is an index i such that p(i), p(i + 1) ≥ 0.17 and (i/t ) >τ ∈ (4/7, 5/7) then 5: Stop the for loop and set θ ϕ = (i/t ) ≤τ . 6: Let U ϕ = e 2πiϕ U . Apply phase-to-amplitude conversion on the exponentiated unitary U t ϕ with precision /4 to obtain a unitary V such that 7: Compute the state V |ψ |0 ⊗a and measure its last a qubits in the standard basis. If the result is the all-zero string then set b = 0 else set b = 1/(2π).
Proof. We decompose the analysis of Algorithm 1 into three parts. First, we show that the two phases ϕ ∈ (0, 1/t) and θ ϕ ∈ (0, 1) obtained at the end of steps 1-5 satisfy the equations with probability at least 1 − /8. We refer to this part as exact phase estimation since θ ϕ is exactly equal to the first τ bits of the shifted phase θ + ϕ. Next, conditioned on any pair (ϕ, θ ϕ ) satisfying the above equation, we prove that the estimate b obtained at the end of steps 6-7 satisfies We refer to this part as unbiased phase estimation. Finally, we analyze the expectation and variance of the final estimate constructed at step 8. Exact phase estimation (steps 1-5). For each phase shift ϕ and integer i ∈ {0, . . . , t − 1}, let ϕ i denote the distance between the angles θ + ϕ and i/t along the unit circle. The two smallest distances are achieved by ϕ i , ϕ i +1 ≤ 1/t where i = t (θ + ϕ) ≤3τ , whereas the other indices j / ∈ {i , i + 1} are at distance at least ϕ j > 1/t . By Hoeffding's inequality and a union bound over the t possible indices, at each round of step 3 the empirical frequency vector p is at distance p − p ∞ ≤ 0.05 from the probability vector p with probability at least 1 − /72. Since step 3 is executed at most 9 times, all of the frequency vectors computed by the algorithm satisfy this bound with probability at least 1 − /8. On the other hand, there is at least one value of ϕ (among the 9 possible ones) such that ϕ i , ϕ i +1 ≤ 5/(8t ) and (θ + ϕ) >τ ∈ (1/2, 3/4). The algorithm will correctly detect such a value because of the gaps proved in Corollary A.2 between the probabilities p(i ), p(i + 1) and the probabilities p(j) for j / ∈ {i , i + 1}. This proves Equation (3).

Amplitude estimation
In this section, we describe a variant of amplitude estimation that returns a (nearly) unbiased estimate and restores the initial state with high probability. We start with a simple procedure, adapted from the Marriott-Watrous scheme [MW05], to sample nondestructively from a Bernoulli variable of parameter Π|ψ 2 given a projector Π and a single copy of a state |ψ .   Proof. Let p = Π|ψ 2 denote the probability to be estimated and define r = 2p(1 − p). By the standard Marriott-Watrous analysis [MW05], the residual state at the end of step 2 is |ψ with probability 1 − r (in which case steps 3-5 need not be executed) and the probability that one iteration of steps 4-5 succeeds in restoring |ψ is r. Thus, the expected number of iterations needed to restore |ψ is r 2 ∞ i=0 (i + 1)(1 − r) i = 1.
We now extend the above algorithm to a low-variance unbiased estimator.
Proof. By a simple application of Proposition B.2, with probability at least 1 − /8, the amplitude estimate q computed at step 1 satisfies We analyze the rest of the algorithm conditioned on the above inequalities being true. Similarly to the proof of Theorem 2.2, by the laws of total expectation and total variance, when removing the condition the expectation and variance of the estimate p will differ respectively by at most /4 and 7 /8 from the values computed below. The amplification parameter τ at step 2 satisfies τ Π|ψ = τ √ p ≤ 1/2 as required by Lemma A.3. Thus, the norm of the post-amplification state |ψ projected along Π is We split the analysis of the subsequent steps of Algorithm 3 into two cases depending on whether the biased estimate q is smaller or larger than 1/t 2 .

Mean estimation
We construct first a quantum estimator that decreases the variance quadratically faster than classically while keeping the estimate (nearly) unbiased. Moreover, it requires a single qsample, which is not destroyed with high probability. The main ingredient is the nondestructive amplitude estimation algorithm constructed in the previous section. The estimator is based on expressing the mean as a linear combination of carefully chosen Bernoulli random variables (with parameters µ + j and µ − j defined at steps 4 and 7 of Algorithm 4) which are related to different truncations of the input random variable. We use a more sophisticated truncation pattern than in previous work [Hei02; Mon15; Ham21; CHJ22] to ensure that the estimate is nearly unbiased. Our estimator requires knowing a rough estimate m of the mean and an upper-bound σ on the variance of the random variable. We explain how to find m nondestructively in Proposition 3.2, whereas σ will be provided directly by our applications. Use a controlled rotation to transform |π X into the state and denote µ + j = E X a j 1X ∈(a j−1 ,a j ] = (I ⊗ |1 1|)|π + j 2 .

7:
Repeat steps 4-6 with the state Proposition 3.1 (Unbiased mean estimator). Let X be a finite random variable with mean µ and variance σ 2 encoded as a qsample |π X . Given two reals ∈ (0, 1) and t ≥ 1 and two estimates m, σ such that |µ − m| ≤ 17σ and σ ≥ σ, the unbiased mean estimator Unbiased(X, t, m, σ, ) (Algorithm 4) outputs a mean estimate µ such that The algorithm needs one copy of |π X , which is restored at the end of the computation with probability at least 1 − , and O(t log(1/ ) 2 ) applications of the reflection I − 2|π X π X | in expectation.
Proof. We start by giving several properties of the shifted random variableX = X − m. First, its second moment is at most since, by assumption, the approximate median m satisfies |µ − m| ≤ 17σ. Next, the expected value of |X| over any interval (a j−1 , a j ] is at most E |X|1 |X|∈(a j−1 ,a j ] ≤ E X 2 1 |X|∈(a j−1 ,a j ] /a j−1 .
Hence, the normalized means µ + j and µ − j satisfy when j > 0.
Lastly, the expectation of |X| over the interval [a k , +∞) is at most where we used Equation (6) and the assumption that σ ≥ σ.
We now analyze the expectation and variance of the final estimate µ. Each of the µ + j and µ − j estimates introduces a bias of at most by Theorem 2.4 and there is another bias due to the truncated part E[X1 |X| ≥ a k ] analyzed in Equation (8). Overall, we have |E[ µ] − µ| ≤ k j=0 a j + 290σ 2 k ≤ (2 k+1 + 290 2 k ) σ ≤ σ. Since the estimates µ j are independent random variables, the variance is Var[ µ] = k j=0 a 2 j Var[ µ j ]. Thus, by Theorem 2.4 and Equation (7), By Theorem 2.4, the algorithm uses O kT log(1/ ) = O(t log(1/ ) 2 ) reflections through |π X and it preserves the copy of |π X with probability at least 1 − 2(k + 1) ≥ 1 − .
We now explain how to find an estimate m that satisfies the first requirement |µ − m| ≤ 17σ needed to apply the above estimator. If we had access to classical samples then it would suffice to compute the average of O(log 1/η) samples to get the above inequality with probability 1 − η. Nevertheless, we do not know how to extract a classical sample from a qsample |π X without destroying it. Since the mean is always within one standard deviation of the median, we can instead attempt to estimate the median. There exist quantum algorithms [NW99;Nay99;Ham21] for computing the quantiles of a distribution, however they also require measuring qsamples. We follow a new approach that incurs a logarithmic overhead in the size of the support but is nondestructive. Our algorithm consists of performing a binary search over the support of the random variable, where we test if a value x is close to the median by estimating the probability Pr[X ≥ x] with the nondestructive amplitude estimation algorithm.
1: Let x 1 < · · · < x n denote the values in the support of X. 2: Set a = 1 and b = n + 1. Run the following binary search until a = b:
Proposition 3.2 (Nondestructive median estimator). Let X be a finite random variable with support size n, mean µ and variance σ 2 . Given a real η ∈ (0, 1), the nondestructive median estimator Median(X, η) (Algorithm 5) outputs an approximate median m such that with probability at least 1 − η. The algorithm needs one copy of |π X , which is restored at the end of the computation with probability at least 1 − η, and O(log(n) log(log(n)/η)) applications of the reflection I − 2|π X π X |.
Proof. By Theorem B.2 and a union-bound argument, all the estimates p k computed by Algorithm 5 satisfy | p k − p k | ≤ √ p k /(3 √ 2) + 1/18 with probability at least 1 − η. In particular, if p k ≥ 1/2 then p k ≥ p k − √ p k /(3 √ 2) − 1/18 ≥ 5/18 and if p k ≤ 1/18 then p k ≤ p k + √ p k /(3 √ 2) + 1/18 ≤ 1/6. Thus, the binary search stops with a = b satisfying p a < 1/2 ≤ 9p a−1 with probability at least 1 − η (note that a ≥ 2 since p 1 = 1). Consequently, if we define the quantile function Q X (p) = sup{x ∈ R : Pr[X ≥ x] ≥ p} for all p ∈ [0, 1], we have Q X (1/2) ≤ m ≤ Q X (1/18). Finally, the result |µ − m| ≤ 17σ follows from the quantile We describe our final estimator for estimating the mean of a product of independent random variables X 1 , . . . , X with relative error . We make the standard assumptions that the relative second moment E[X i ] 2 and the inverse fidelity 1/| π X i |π X i+1 | 2 are upper-bounded by some known value B, which is a constant in our applications (Section 4). Our estimator uses a similar approach to that of Dyer and Frieze [DF91], by taking the product of (nearly) unbiased estimates µ 1 , . . . , µ of E[X 1 ], . . . , E[X ] respectively, each with variance O( 2 / ). We obtain a quantum speed-up by using the quantum unbiased estimator of Proposition 3.1 to reduce the variances quadratically faster than classically.

3:
Compute an estimate m i of the median of X i by using the nondestructive median estimator Median( X i , 1/(11 )). Set σ i = m i B.

4:
Compute an estimate µ i of the expectation E[ X i ] by using the unbiased mean estimator
Proof. By a union-bound argument, with probability at least 10/11, all the median estimates m i computed at step 3 are within distance 17 Var[ X i ] of the means E[ X i ]. We condition on this event happening for the rest of the proof.
We start by analysing steps 2-5 for a given random variable X i . First, notice that we can construct a qsample for X i by taking |π X i = U average (|π X i ⊗K ⊗ |0 ) where U average is a unitary computing the average of the first K registers in the last register. The expectation is unchanged E[X i ] = E[ X i ], whereas the relative second moment and the fidelity are since | π X i |π X i+1 | 2 = | π X i |π X i+1 | 2K by properties of the fidelity. The next equation shows that the distance between the rough estimate m i computed at step 3 and the expectation is bounded by a multiple of the standard deviation (by Proposition 3.2), which in turn is a small multiple of the expectation (since the relative variance is small by Equation (9)).
Finally, the expectation of each estimate µ i computed at step 4 is bounded as where the first inequality is by Theorem 2.4 and the second is by σ i = B m i and Equation (10). Consider now the product random variable X = X 1 · · · X . The expectation of the final estimate µ is at most where the second step is by Proposition 3.1 and the fourth step uses that m i /E[ µ i ] ≤ 3 by Equations (10) and (11). We conclude by using the triangle and Chebyshev inequalities that ( E[ µ]) 2 ≥ 3/4. The algorithm needs K copies of the initial qsample |π X 1 to anneal into K copies of |π X i at the i-th stage of steps 2-5. Each transition from |π X i ⊗K to |π X i+1 ⊗K requires O(B K ) measurements in expectation by Lemma A.6. Furthermore, at each stage, step 3 uses O(log(n) log( )) reflections through |π X i by Proposition 3.2 and step 4 uses O(B √ / ) reflections by Proposition 3.1. Overall, the number of reflections through |π X 1 , . . . , |π X is

Partition function estimation
In this section, we showcase an application of our newly-constructed quantum mean estimation algorithm. Specifically, we show how it can be used to speed up existing quantum algorithms for estimating partition functions. In Section 4.1, we elaborate on the generic algorithm for partition function estimation and how our results provide a speed up and, in Section 4.2, we discuss how this gives rise to more efficient quantum algorithms for several applications.

General algorithm for partition function estimation
In this subsection, we first describe our partition function estimation algorithm in high level, and define the required notation along the way. We follow the exposition in [HW20], but use slightly different notation to match the rest of this document more closely.
Let Ω be a state space, and let H : Ω → {0, . . . , n} be a Hamiltonian, i.e., a function associating a numerical value to each of the states. Our goal will be to compute the number of states whose value is 0, i.e., |H −1 (0)|. To that end, we define the partition function Z : and we observe that Z(0) = |Ω|, and Z(∞) = |H (−1) (0)|. The partition function arises frequently in statistical physics, where β is referred to as the inverse temperature. Next, we define a sequence of inverse temperatures 0 = β 0 < β 1 < · · · < β = ∞, and express Z(∞) as the product Since throughout the sequence of β i 's, we are increasing the inverse temperature, and hence decreasing the temperature, this sequence is referred to as the cooling schedule. Its length is called the schedule length. The core idea for estimating Z(∞) is to evaluate each of the factors in the product on the right-hand side of Equation (12) individually. Thus, let i ∈ [ ]. We endow Ω with a probability distribution π i , called the Gibbs distribution in statistical physics, and define a random variable X i : Ω → R on it, with It follows immediately that The idea is now to devise a procedure that samples from the Gibbs distribution and obtains an estimator In order to be able to sample from the Gibbs distribution, we encode it in a quantum state. To that end, we define the Gibbs state |π i ∈ C Ω at inverse temperature β i , defined as In typical applications (see the next subsection), it can be much easier to reflect through the Gibbs state than to prepare it. Therefore, we estimate E[X i ] using a nondestructive mean estimation algorithm, which we constructed in the previous sections. This has the benefit of having access to the Gibbs state even after the computation of µ i , which makes it possible to reuse it for computing the next factor in the product on the right-hand side of Equation (12). It remains to choose the cooling schedule, i.e., the inverse temperatures β 1 , . . . , β −1 . For reasons that are sketched on high level, we want our cooling schedule to have the following two properties.
1. B-Chebyshev. For any B ≥ 1, a cooling schedule is called B-Chebyshev, if for any two subsequent β i−1 and β i , we have This requirement can be viewed as a tail bound -it tells us that X i cannot have too much weight above E[X i ]. Indeed, by Markov's inequality, It ensures that we can not only give a good additive estimate of E[X i ], but also a good multiplicative one.
2. B-slowly varying. For any B ≥ 1, a cooling schedule is called B-slowly varying if This requirement can be seen as a proximity requirement -it tells us that the Gibbs states |π i and |π i+1 have some non-negligible overlap. It ensures that we can anneal, i.e., transform, |π i into |π i+1 without having to do too much work, and hence circumvents having to prepare |π i+1 from scratch for each of the factors on the right-hand side of Equation (12).
When we set B = e 2 , it is shown in [ŠVV09] that there exists a cooling schedule of length = O( log |Ω| log(n)), which is B-Chebyshev. Subsequently, it was shown [HW20] that one can modify the construction so that the resulting cooling schedule is not only B-Chebyshev, but also slowly varying. Moreover, computing this cooling schedule can be done on the fly, with associated cost scaling approximately linearly in the schedule length.
Combining previous work with our product estimator, Theorem 3.3, gives rise to the following result.
Proof. We run the algorithm from [HW20] to compute the cooling schedule on the fly. The total number of Gibbs state reflections used in this step is O( log |Ω| log(n)) [HW20, Theorem 13], and the resulting cooling schedule is both e 2 -slowly varying and e 2 -Chebyshev, and of length O( log |Ω| log(n)). For computing the product in Equation (12), we use the product estimator, Theorem 3.3. The total number of Gibbs state reflections used in this algorithm is O(log 3/4 |Ω| log 3/4 (n)/ + log |Ω| log(n)). With standard failure probability reduction techniques, we obtain that the number of Gibbs state reflections performed by the resulting algorithm scales as the sum of the two complexities, which yields the desired result.
The core of our improvement lies in our product estimator -the one used in [HW20] is quadratic in the schedule length, whereas ours is subquadratic. Thus, more generally, if we are in a setting where a cooling schedule that is shorter than O( log |Ω| log(n)) suffices, then our resulting algorithm scales with the schedule length to the power of 3/2. A typical setting where this is the case is when we a priori know a lower bound on |H (−1) (0)| -then we can terminate the annealing at a lower inverse temperature, and hence obtain a shorter cooling schedule. We leave working out the details in this setting for future work.

Applications
The quantum partition function estimator can be applied to several counting problems. See [HW20] for more possible applications. In all applications we list in this section, n scales at most polylogarithmically in |Ω|, so we can hide any dependence on log(n) in the tilde of the big-O notation.
Counting k-colorings Let G = (V, E) be an undirected simple graph with degree at most ∆, and suppose we have k colors. We would like to count the number of ways we can color the vertices such that no two adjacent ones have the same color. We let Ω be all the possible colorings, which means Z(0) = |Ω| = k |V | . We define H(x) to be the number of monochromatic edges in G. Then, Z(∞) = |H (−1) (0)| is exactly the number of k-colorings. Jerrum [Jer95] constructed a Markov chain whose stationary distribution is the Gibbs state at infinite temperature, and analyzed that its mixing time is O(|V | log |V |). With a metropolis sampling correction this can be turned into Gibbs states at arbitrary temperatures, and thus performing approximate Gibbs state reflections takes time O(1/ |V |). Multiplying by the complexity in Theorem 4.1, we can obtain an -multiplicative approximation to Z(∞) in time O(log |Ω| 3/4 |V |/ ) = O(|V | 5/4 / ).
Ising model Let G = (V, E) be an undirected simple graph of maximum degree ∆, and assign a sign ±1 to each of the vertices. We let Ω be the set of all assignments of signs to the vertices, which readily implies that |Ω| = 2 |V | . We let H(x) be the number of edges whose endpoints have the same sign. Then, Mossel and Sly [MS13] construct a Markov chain whose stationary distribution is the Gibbs state at arbitrary temperature. They also show that it mixes in time O(n log(n)) as long as (∆ − 1) tanh(β) < 1. Thus, for any of these inverse temperatures β, we can evaluate Z(β) up to relative error in time O(|V | 5/4 / ).
Counting matchings Let G = (V, E) be an undirected simple graph. A matching M is a subset of edges such that every vertex v ∈ V is incident to at most edge in M . We let Ω be the set of all matchings, and H(M ) = |M |, i.e., the number of edges in the matching. Then Z(0) = |Ω| is the number of matchings, and Z(∞) = 1. Now, the setting is slightly different than before, since the partition function is easy to calculate at zero temperature, rather than at infinite temperature. The core idea is to run the same algorithm as outlined before, but anneal in the opposite direction. Previous works [HW20; Mon15] justify this idea in more detail. Jerrum and Sinclair [JS89] constructed a Markov chain whose stationary distribution is the Gibbs state at zero temperature, and analyzed that the mixing time of this Markov chain is O(|V ||E|). Similarly as before, the Markov chain can be easily adapted by a metropolis sampling correction so that its stationary distribution is a Gibbs state at arbitrary temperature, and so we can perform Gibbs state reflections in time O( |V ||E|). The state space of this random walk is of size |Ω| = O(|V |!2 |V | ), which implies that log |Ω| = O(|V | log |V |). Thus, combining with the complexity in Theorem 4.1, we can find a -multiplicative estimate of the number of matchings in time O(|V | 5/4 |E| 1/2 / ).
Counting independent sets Let G = (V, E) be an undirected simple graph, with degree at most ∆. An independent set is a subset of vertices I ⊆ V such that none of the edges has both endpoints in I. Similarly as before, we let Ω be the set of all independent sets, and we let H(I) = |I|. Then, as in the setting where we were counting matchings, we know Z(∞) = 1, and we want to compute Z(0) = |Ω|. Vigoda [Vig01] constructed a Markov chain whose stationary distribution is the Gibbs state at arbitrary temperature, and analyzes the mixing time to be O(|V | log |V |). Furthermore, |Ω| ≤ 2 |V | , and hence our algorithm runs in time O(|V | 5/4 / ).
Computing the volume of a convex body We can also use our techniques to speed up one of the costly subroutines in the volume estimation algorithm of convex bodies from [CCH+19]. To that end, let r, R > 0, and let K ⊆ R d be a convex set with B(0, r) ⊆ K ⊆ B(0, R). Let D = R/r, and without loss of generality assume that B(0, 1) ⊆ K ⊆ B(0, D).
The algorithm from [CCH+19] first estimates the volume of K with O(d 3 + d 2.5 / ) queries to O K , and then estimates the ratio between the volume of K and K with another O(d 2.5 / ) queries. The first part can be sped up with our methods to O(d 3 + d 2 / ), but it's unclear if our techniques help to improve the second part. We leave this for further research. and ∆ i = min{|θ − i/t|, |1 + θ − i/t|, |1 − θ + i/t|} is the distance between θ and i/t along the unit circle. The algorithm needs one copy of |ψ , which is restored at the end of the computation, and O(t) controlled-U operations.
The next two results provide efficient quantum algorithms for converting between phase and amplitude encodings of a real parameter.
Finally, we explain how to anneal from one state to the next, given that there is a nonnegligible overlap between the two. The idea is due to Marriott and Watrous [MW05] and is similar to the coin flipping algorithm (Proposition 2.3).
Algorithm 7. Annealing from |ψ to |φ , Anneal(|ψ , |φ ) 1: Start with the state |ψ . the amplification parameter η. The probability of measuring |ψ |0 in the first two registers of the resulting state is at least 1 − δ, in which case the last register collapses to 1 √ i π(i) 2 i π(i)|i .
We obtain below the same statement as [HW20] for performing nondestructive amplitude estimation. Furthermore, we show that the amplitude estimate is zero with high probability when the amplitude is sufficiently small (as is the case for the regular amplitude estimation [BHMT02]).