Quantum Monte Carlo Integration: The Full Advantage in Minimal Circuit Depth

This paper proposes a method of quantum Monte Carlo integration that retains the full quadratic quantum advantage, without requiring any arithmetic or quantum phase estimation to be performed on the quantum computer. No previous proposal for quantum Monte Carlo integration has achieved all of these at once. The heart of the proposed method is a Fourier series decomposition of the sum that approximates the expectation in Monte Carlo integration, with each component then estimated individually using quantum amplitude estimation. The main result is presented as theoretical statement of asymptotic advantage, and numerical results are also included to illustrate the practical benefits of the proposed method. The method presented in this paper is the subject of a patent application [Quantum Computing System and Method: Patent application GB2102902.0 and SE2130060-3].


Introduction
Monte Carlo integration (MCI) solves the problem of estimating the mean of some probability distribution by averaging samples therefrom. This could either be an average of the samples themselves, or of the samples with some function applied first. It is known that a quadratic quantum speed-up in MCI is available [2] by calling quantum amplitude estimation (QAE) [3] as † Contact: Steven.Herbert@Quantinuum.com After posting a pre-print of this paper, we were made aware of an existing work, Ref. [1], which proposes a somewhat similar decomposition of the Monte Carlo integralalthough we note that Ref. [1] lacks the detailed algorithmic specification and rigorous proof of performance that are central to the present work. a sub-routine. Such is the ubiquity of Monte Carlo methods throughout the physical, biological, data and information sciences, that this has, in turn, spawned a great deal of interest in quantum Monte Carlo integration (QMCI) -most notably in applications related to finance such as option pricing [4][5][6][7][8][9][10][11][12][13][14].
However, the theoretical possibility of using quantum computing to perform MCI more efficiently has not masked a number of clear obstacles that are in the way of actually doing so in practice. Most notably the quantum advantage in MCI is an advantage in query complexity: the number of uses of a quantum state that encodes the probability distribution in question will be fewer than the number of classical samples from that distribution -but this immediately raises the question of how to prepare such a quantum state. Recently it was shown that there is no advantage when the oft-cited Grover-Rudolph method [15] is used to prepare states encoding probability distribution when a full audit of the computational complexity is undertaken [16].
Even should this problem be resolved, there are a number of reasons to be skeptical about the possibility of QMCI yielding a useful quantum advantage in the NISQ era. Most notably, QAE, as originally described, requires use of the quantum phase estimation (QPE), which is unlikely to be practical on near-term quantum hardware. On this front, though, there are reasons for optimism -it has been shown that QPE can be omitted in favour of classical post-processing to extract an estimate of the statistical quantity of interest from measurements of the quantum state [17][18][19][20]. Moreover, recent efforts have been made to further cut the required circuit depth by interpolating between classical and quantum MCI [21,22].
To a certain extent, however, these advances have served only to reveal further issues with NISQ compatible QMCI: specifically that the "natural" quantity to estimate on a quantum computer is: where X are samples from p(x). Any other quantity, including the mean of p(x) itself, will incur a significant amount of arithmetic to be performed on the quantum computer. Indeed, this cost is not unlike that of using QPE -not a problem asymptotically, but almost certainly prohibitive in the NISQ era. The foremost proposal to work around the need for quantum arithmetic is showcased in Option Pricing using Quantum Computers [23] in which a technique proposed by Woerner and Egger [10] is used. Specifically, the technique shifts and rescales the support of the probability distribution such that for all x : p(x) = 0, sin 2 (x+π/4) ≈ x+1/2, thus enabling an estimate of the mean to be extracted. However, this method has an adverse effect on the rate of convergence, which can be seen as a false economy: the requirement of quantum arithmetic has been removed to reduce the circuit depth, but this comes at a cost of reduced convergence rate, and so correspondingly deeper circuits will be required to achieve the specified estimation accuracy. In this article we propose an alternative solution that still eliminates the need for quantum arithmetic, but also upholds the full quadratic quantum advantage. Instead of squeezing the support of the probability distribution to correspond to an approximately linear region of a trigonometric function, instead we propose that the function applied to the samples is extended as a periodic, piecewise function. This periodic function can be decomposed as a Fourier series, whose various components have means that can be estimated individually and then recombined. With judicious choices about how to do this Fourier series decomposition, we show that the quadratic quantum advantage is indeed retained.
Furthermore, by performing the QMCI in this way, there is no "favoured" quantity to compute (in the sense of (1)), and computation of the mean of any function applied to the samples is equally easy in principle. This immediately suggests a way in which an actual computational advantage, as opposed a theoretical advantage in query complexity, may be realised using QMCI. Quantum advantage can be realised using the proposed method if: (a) there exists some probability distribution that can be exactly prepared as a quantum state in shallow-depth; (b) a function is then applied to the random samples such that the computed mean is a quantity that cannot be analytically evaluated. This is a clear and useful distinction over, for example, the method proposed by Stamatopoulos et al where the distribution of interest must be directly encoded into the quantum state [23].
Finally, it is worth highlighting that our proposed method of QMCI can also be used to compute the expectation of the product of functions of two correlated random variables, again with the full quadratic quantum advantage and without requiring quantum arithmetic. This may find application in, for example, the computation of correlations and covariances of random variables from multivariate probability distributions. Table 1 summarises the benefits of the method of QMCI proposed here compared to existing alternatives. Convergence is measured in meansquared error (MSE) as a function of the number of samples from the probability distribution, q. The table includes classical MCI for comparison, but it should be noted that there are classical alternatives, for example quadrature and quasi-Monte Carlo methods. Ostensibly, these appear to have favourable performance (achieving MSE convergence proportional to q −2 ), however they suffer the curse of dimensionality: the complexity grows exponentially with the number of dimensions. Hence it is noteworthy that our proposed method applies to random variables marginalised out from multivariate probability distributions as, for a sufficiently large number of dimensions, classical MCI will be the most appropriate classical algorithm to benchmark against.
The remainder of this paper is organised as follows: in Section 2 the problem at hand is specified, the relevant definitions given, and precise details of the motivation presented; in Section 3 the main theoretical result is given; in Section 4 the result pertaining to the expectation of products of random variables is given as a Corollary, whose proof is then given in Appendix A; in Section 5 a simulated example is given as a numerical proof of principle; and in Section 6 the results of the paper are discussed. The method presented in this paper is the subject of a patent application [24].

Preliminaries and Motivation
Suppose we have the d-dimensional discrete probability distribution p (that is a probability distribution over some finite subset of R d ). The points of non-zero probability mass are spaced at constant intervals and for simplicity, we let the number of points along each axis be a power of 2, say 2 N j for the j th dimension (we can "zero-pad" so there is no loss of generality in this assumption). The probability distribution is encoded in a quantum state, |p , of the form: (1) , . . . , x (d) ) |x (1) . . . x (d) (2) and we assume that the quantum algorithm has access to a circuit, P , that prepares |p , ie, |p = P |0 (where |0 denotes the tensor product of an appropriate number of qubits in the |0 state).
QMCI estimates the expectation of some bounded function, f (x), of samples from the marginal distribution of a single dimension of the multivariate probability distribution. This dimension must therefore be specified -let this be the i th dimension -so that QMCI returns an estimate of the mean of the probability distribution implicitly defined by f (x) applied to samples, X, from p(x (i) ), that is an estimate of: To achieve this, the support of p(x (i) ) is also taken as an input to the QMCI algorithm: this is given as the value that the first point of the distribution corresponds to, x l and the spacing interval, ∆. It is also useful to express the maximum value of the support of the distribution, In order to obtain a quantum advantage using the QMCI method in this paper, we further require that f (x) is continuous, has continuous first derivative and bounded piecewisecontinuous second and third derivatives over the support of p(x (i) ) ('the support of p( Thus A |0 = cos θ a |Ψ 0 |0 + sin θ a |Ψ 1 |1 where |Ψ 0 and |Ψ 1 are just some quantum states and sin 2 θ a = x (i) f (x (i) )p(x (i) ), ie, the expectation we are seeking (see Ref. [17] for a more detailed exposition of this), which can be estimated by calling some QAE algorithm as a subroutine 1 .

Definition 1. A QAE algorithm takes as an input a circuit
A where A |0 = cos θ a |Ψ 0 |0 + sin θ a |Ψ 1 |1 , and returns an estimate of s = sin 2 θ a . The allowed number of uses, denoted "q", of the circuit A is given as an input to the QAE algorithm.
From this definition, we can write any QAE algorithm as a function in pseudocode form: whereŝ expresses the fact that this is an estimate of s. It is also useful to formally define the MSE, The convergence of any QAE algorithm can be expressed as:ˆ for some λ, and the method proposed in this paper assumes that the convergence rate of the QAE algorithm as a function of the number of uses of A is known. That is, there exists a constant k 1 and we know the constant λ such that: A classical MCI algorithm can be used to perform QAE, in the sense of the above definition, and in this case λ = 1, whereas typical QAE algorithms have λ = 2. However, the method we propose here works for any QAE algorithm, including a proposed method that interpolates between the classical and quantum (so 1 ≤ λ ≤ 2) [21,22,25]. QAE algorithms use the circuit A (and hence P ) to build a "Grover iterate" circuit, where S 0 and S χ do not depend on P or A. Further details do not directly concern us hereapart from that circuit depth of QMCI is usually counted in number of sequential uses of Q [3,17]. 1 Strictly speaking, it is required that 0 ≤ f (x) ≤ 1 for all x, and so if this is not the case then some transformation must be performed. For instance, in the case where f (x) = x (that is the mean is being estimated) then a simple affine transformation suffices, but more generally this may be trickier.

Problems with existing methods
We can see that performing QAE when A = R(P ⊗ I) returns an estimate of s = f (x)p(x), and thus does yield a quantum speed-up. However, this comes at the cost of using the circuit R, which encodes the function applied to the random samples, and in general this will be prohibitively complex -even when QPE-free forms of QAE are used. Even the (classically) trivial case where f (x) = x (ie, we are just finding the mean of p(x)) would entail significant amounts of arithmetic to be performed quantumly. An important exception is the special case where f (x) = sin 2 (mx + c) for some constants m and c, which can be achieved by a bank of R y rotation gates, as shown in Fig. 1.

QMCI of a Single Random Variable
The main result of the paper requires QAE to be performed on the quantum circuit A(P, i, β, n, ω) given in Fig. 1. P is a circuit that prepares from |0 a quantum state encoding a multivariate probability distribution, as in (2), and i denotes the dimension of interest (we assume that the support of the i th dimension is known, in terms of the value of the first point of non-zero probability mass, x l and spacing between points of probability mass, ∆); n is an integer and ω and β real numbers.

Proposition 2.
Let QAE(A, q) be a QAE algorithm in the sense of Definition 1, which has MSE convergence parameterised by λ as in (6).
i. 1−2QAE(A(P, i, 0, n, ω), q) is an estimate of: ii. 1 − 2QAE(A(P, i, π/2, n, ω), q) is an estimate of: Proof. The function of circuit A is to apply an R y rotation to the last qubit that is proportional to some constant plus the value of the i th register of P . Recall that the circuit P , when applied to |0 , encodes a multivariate probability distribution in the computational basis states of a quantum state, as defined in (2), thus from Fig. 1 we can see that, when β = 0, A |0 prepares the state: and thus the amplitude of the |1 component of the final qubit is: where p(x (i) ) is the marginal distribution -note that the marginalisation has been achieved automatically. Thus QAE(A(P, i, 0, n, ω), q) estimates the value of 0.5(1 − x (i) p(x (i) ) cos(nωx (i) )) with MSE ∈ Θ(q −λ ), by (6).
So it follows that 1 − 2QAE(A(P, i, 0, n, ω), q) therefore estimates x (i) p(x (i) ) cos(nωx (i) ), and as this differs from the original estimate only by being shifted and scaled by a constant factor, the estimation accuracy remains such that MSE ∈ Θ(q −λ ).
Turning now to the second part of the Proposition, by the same reasoning as before, when β = π/2, A |0 prepares the state: and thus the amplitude of the |1 component of the final qubit is: Proposition 2 now allows us to proceed to the main result of the paper.
) and f is a function that is continuous in value and first derivative and whose second and third derivatives are piecewise-continuous and bounded, can be estimated with MSE ∈ Θ(q −λ ), where q is the number of uses of a circuit preparing the quantum state |p (ie, as defined in (2)) and λ is the convergence rate of some QAE subroutine (ie, as defined in (7)) which operates on circuits of the form defined in Fig. 1.
Remark. This result substantiates the "best of all worlds" claim made in Section 1: the convergence of the QMCI estimate is equal to that of the underlying QAE algorithm; whilst the circuit in Fig. 1 requires only a single bank of rotation gates to construct A from P , which is minimal. That is, if the probability distribution of interest is over 2 N i points then at least N i two-qubit gates are required in order that the amplitude of a further qubit is equal to the mean of some function the probability distribution (ignoring unnatural cases where some of the data-points of the probability distribution are omitted altogether).
Proof. We prove the theorem by giving an explicit construction. The proof strategy is to first build a periodic function, f(x), from f (x); then to show that the expectation can be approximated by performing QAE on the components of a Fourier series decomposition of f(x) using Proposition 2; and finally to show that the total number of uses of P can be distributed amongst the Fourier series components such that the claimed overall convergence rate holds. First, we build a periodic piecewise function is also constructed such that it is continuous in value and first derivative, and has second and third derivatives that are piecewisecontinuous and bounded. For example, some f(x) of the following form, which repeats with period xũ − x l , is suitable in general: where xũ > x u andf (x) is itself sufficiently smooth (that is, continuous in value and first derivative, and has second and third derivatives that are piecewise-continuous and bounded) and It is easy to see that a cubic will always suffice to fit such ã f (x), and on occasion a simpler function may do. As f(x) is periodic, it has a Fourier series: a n cos(nωx) + b n sin(nωx) (16) where ω = 2π/T and T is the period of the periodic piecewise function. By construction, f(x) is continuous in value and first derivative and has second and third derivatives that are piecewisecontinuous and bounded, and so the Fourier series coefficients decay as 1/n 3 or faster 2 and thus we can write: where |ã n |, |b n | ∈ O(1). As f(x) = f (x) whenever p(x (i) ) = 0, we can express: To estimate µ, we can thus estimate each of the parenthesised sums in the final line of (18) individually. These can be estimated using QAE according to the procedure described in Proposition 2. To analyse the overall accuracy of the consequent estimate of µ, we first define: and further defineμ n respectively such that: where ν (a) n and ν n are biases (noting that the possibility of biased estimation has not been precluded) and n are zero-mean random variables. The inaccuracy of estimating µ in this way does not solely arise because of the inaccuracy of the estimates of each term in (18), we must also address the fact that the Fourier series is an infinite sum, and so obviously we cannot estimate all of the terms in practice. Therefore we truncate at some n max . That is, our estimate, µ, of µ is:μ = c + and to keep track of the error incurred by this truncation, we define the value of the truncated part of the sum: and it is also convenient to define: Substituting (19) -(24) into (18) and also using the definitions (25) and (26) we can express the actual mean, µ, in terms of our estimate of the mean,μ: from which we can express the MSE: using the fact that n are zero-mean by design, and also that estimates corresponding to the various terms in the Fourier series are independent.
To prove the theorem it is now incumbent upon us to show that the number of uses of P assigned to the QAE estimation of each of the terms in (18) and n max can be chosen such that the claimed convergence holds. To do this, we first let q n be the number of uses of P allowed in the QAE estimate of each of the sine and cosine components of the n th term in the Fourier series such that: q n = q 0 n −κ (29) for some constants q 0 and κ the latter of which we will later give numerical value. Recalling that |ã n |, |b n | ∈ O(1), we can define a constant k 2 such that k 2 ≥ |ã n |, |b n | (in the following analysis we will need to introduce further new constants, and will henceforth do so without explicit announcement by using k with subscript incremented relative to the most recently introduced constant). We now address the terms in (28). First, we note that the magnitude of the bias, ν n , is upper-bounded by the square root of the MSE, which has been bounded in (7), so (starting with (25) and also using (29), and k 2 ≥ |ã n | as defined above): so long as κ is chosen such that κ < 4/λ and the sum therefore converges. In exactly the same way we get: Turning now to the term (of (28)), E ( (a) n ) 2 , which is the mean-squared error of the unbiased part of the estimate of µ (a) n , and is therefore upper-bounded by k 1 q −λ n (this would hold with equality for an unbiased estimator and be strictly less for an estimator with non-zero bias). So we get (again using (7)): where the previously specified condition, κ < 4/λ, suffices to ensure that this sum also converges. Likewise we get The final quantity in (28) that we must address is the value of the truncated term, η, defined in (24). The magnitudes of cos(nωx (i) ) and sin(nωx (i) ) are bounded above by 1 for all x (i) , and so and thus from (24), using k 2 ≥ |ã n |, |b n | (by definition) we have Putting this all together by substituting (30), (31), (32), (33) and (36) into (28) we get: From which we can see that in order to achievê 2 ∝ q −λ 0 it suffices to choose n max such that: We must now count up the total number of uses of P . Recall that we have assigned q n uses of P to each of the sine and cosine components of the n th term in the Fourier series, thus we have nmax n=1 2q n uses of P . However, in the above analysis we have allowed q n to take non-integer values, but in practice we must round to an integer. To ensure that the claim of convergence holds, we round up in each case. If all n max Fourier components require q n to be rounded up by an amount approaching 1 for each of the sine and cosine components then we will require an extra 2n max uses of P in the worst case. Thus we can evaluate the total number of uses of P : From (30) we require that κ < 4/λ, and we also have that in practice 1 < λ ≤ 2 (ie, at λ ≤ 1 there is no quantum advantage, and λ = 2 is the Heisenberg limit). Thus choosing any κ < 2 always suffices, and so we let κ = 2 − δ for some 0 < δ < 1. As λ ≤ 2, we have that λ/4 ≤ 1/2, and thus we get: because the sum converges when 0 < δ < 1. Thus we have shown that the total number of queries is proportional to q 0 . Putting this together witĥ 2 ∝ q −λ 0 from (37) and (38) we can thus see that 2 ∈ Θ(q −λ tot ), proving the theorem.
It is important to be careful about the evaluation of the Fourier coefficients. If these can only be computed using numerical integration, then we have essentially shifted the computational load, but not reduced the overall complexity. For many common functions of interest, such as the mean (f (x) = x), second moment (f (x) = x 2 ) and indeed any higher-order moments, then it is easy to see that the Fourier coefficients can be calculated symbolically, and so concerns on this front can be allayed. Furthermore, if there are functions whose Fourier coefficients cannot be found symbolically, but nevertheless are commonly-used, then it may be reasonable to assume that these have been pre-computed to high accuracy and stored, and thus should not count against the complexity of any individual run. However, more generally the Fourier coefficients should be symbolically computable in order that the method attains the advertised speed-up once all of the complexity is counted.

QMCI of the Product of two Correlated Random Variables
As the proposed method of QMCI applies to quantum states that encode multivariate probability distributions, a natural next step is to apply the proposed method to functions of multiple random variables. For the sake of definiteness, and for simplicity, we restrict the following analysis to computing the product of functions of two random variables. Moreover, this setting covers computation of correlation and covariance, and indeed most quantities that may be of interest in practice. However, in principle the following approach of using trigonometric identities to break the Monte Carlo integral down into terms that can estimated using QAE also applies to the product of (functions of) three or more random variables.

Corollary 4. The quantity, µ
) and f and g are functions that are continuous in value and first derivative and whose second and third derivatives are piecewise-continuous and bounded, can be estimated with MSE ∈ Θ(q −λ ), where q is the number of uses of a circuit preparing the quantum state |p (ie, as defined in (2)) and λ is the convergence rate of some QAE subroutine (ie, as defined in (7)) which operates on circuits of the form defined in Fig. 4.
Proof (sketch). The idea is to build periodic piecewise functions for both f and g, and thus expand each as a Fourier series. This results in a double sum each term of which contains products of sines and cosines. In particular, these will always be either cos(ρ) cos(σ), sin(ρ) sin(σ) or cos(ρ) sin(σ) for some ρ, σ. The next step is to use trigonometric identities to express these terms as functions of single sines and cosines (with correspondingly altered angles). Thus we can decompose the sum representing the expectation we need to compute into a double sum over sines and cosines. The analysis then follows along the same lines as that of Theorem 1, and so to avoid tedious repetition we defer the full proof to Appendix A.

Numerical Results
To illustrate the practical utility of our method of QMCI, we present some simulated results. So that we can compare to classical MCI the rescaled QMCI method of Stamatopoulos et al [23] and conventional QMCI (with quantum arithmetic, but still taking advantage of QPE-free forms of QAE) we consider the case where f (x) = x, that is we simply compute the mean of the sampled distribution without applying any other function. We compare for a 16-point univariate probability distribution supported on {−8, −7, . . . , 7}, shown in Fig. 2(a), which is encoded in a quantum state that can be prepared using the circuit shown in Fig.2(b). This choice of probability distribution and state preparation circuit has been made largely for practical reasons: there is little to be gained (in terms of proving the principle) by addressing a multivariate distribution, and this would significantly slow down the simulations; and on a similar note, it is desirable to use some probability distribution whose encoding state can be prepared in a very shallow circuit. We have avoided a distribution whose mean is the median of the support, as in this case the rescaled QMCI method can perform artificially well. The results are shown in Fig. 3, where the simulated results for "Fourier QMCI" (the method we propose in this paper), "Rescaled QMCI" (the method of Stamatopoulos et al ) and "Quantum arithmetic QMCI" (conventional QMCI using arithmetic subcircuits) were averaged over 500 runs. We can see that for sufficiently small desired root-MSE (RMSE), Fourier QMCI outperforms both Classical MCI and Rescaled QMCI, as expected. Indeed, the slope of the line-of-bestfit for Fourier QMCI is −1.02, which is very close to the theoretical value of −1 (note that this is root-MSE, hence −1 not −2 as for MSE), and for rescaled QMCI is −0.771, clearly inferior and not too far off the theoretical value of 0.667. It is also worth highlighting that the "cross-over" point between Fourier QMCI and classical MCI occurs very close to the third data-point of the former. This represents 1100 uses of P in total, but more significantly a circuit depth of only 8 sequential uses of Q in any single quantum circuit.
Finally, we can compare to conventional QMCI (with quantum arithmetic) -as expected, this performs a little better than Fourier QMCI, although most significantly each has the same rate of decay of RMSE with number of uses of P (the line of best fit for "Quantum arithmetic QMCI" is −1.05 ≈ −1). It is however, worth remarking on the discrepancy between the constant factors for each. The lines of best fit have equations: RMSE(Fourier QMCI) = 194q −1.02 RMSE(Quantum arithmetic QMCI) = 107q −1.05 from which we can see that the constant factor is about twice as large for Fourier QMCI as conventional QMCI (194/107 = 1.81), which means that to achieve the desired RMSE the deepest circuits will require twice as many Grover iterates for Fourier QMCI compared to conventional QMCI.
However, this will generally still translate into a significant circuit depth reduction for Fourier QMCI when the complexity of the quantum arithmetic is taken into account. In order to construct a suitable circuit R for the function f (x) = x first the bank of rotation gates would be needed to prepare sin(x) (for instance -or some other trigonometric function), and then arithmetic circuits for the functions arcsine and square root would be needed. Ref. [26] provides optimised circuits for each of these, and for example  Fig. 2(a).
we can see from Ref. [26, Table II] that implementing arcsine in quantum arithmetic would require more than 100 additional qubits and nearly 5000 Toffoli gates, even in the least accurate case considered. Plugging in values to the various expressions for Toffoli counts in Ref. [26] suggests that, even were the required accuracy to be relaxed further, the number of additional Toffolis would still run into the 1000s. Thus we can see that doubling the number of Grover iterates to achieve the desired RMSE (ie, for Fourier QMCI compared to conventional QMC) is a relatively modest price to pay to avoid the heavy cost of quantum arithmetic.

Discussion
In this paper we have proposed a method for performing QMCI by expanding the Monte Carlo integral as a Fourier series and estimating each component individually using amplitude estimation. Our primary motivation for this was to achieve the full quadratic quantum advantage without incurring the need to perform arithmetic operations on the quantum computer, however an added bonus is that our method of QMCI has an innate flexibility that allows it to be readily adapted to a wide variety of settings of realworld relevance. Specifically, we have shown that QMCI can be performed when any function is to be applied to the random samples (with only a mild technical restriction on the smoothness of that function). Moreover, even the expectation of functions of multiple variables can be computedallowing the quadratic advantage to be achieved when estimating, for example, the covariance and correlation of some data.
Will Fourier QMCI yield a useful quantum advantage in the NISQ era?
To answer this question, we first must clarify what, exactly, we mean by "the NISQ era". It is fast becoming clear that the dichotomy between the NISQ era and the Full-Scale Fault-Tolerant (FSFT) era is too simplistic -at least when viewing NISQs in terms of their original definition by Preskill [27]. In these strict terms, the quantum community has very much been focussed on variational algorithms as the only route to NISQ era advantage, however this disposition is somewhat myopic. If the predictions in scaling-up of qubit number (if not fidelity) that hardware manufacturers are making (see eg. Ref [28]) come to fruition, then it is reasonable to look forward to a time when some qubits will be allocated as ancillas to enable some degree of error-correction (if not full fault-tolerance) and others as ancillas to perform entanglement distribution to reduce the depth of otherwise costly qubit routing subcircuits (see eg. Ref. [29,30]), thus effectively facilitating the simulation of a quantum computer with lower qubit-number but higher quantumvolume on the hardware. In such a case, it will in turn be reasonable to think not only of running quantum algorithms with variationallyoptimised ansätze, but also those with explicitlyconstructed circuits, such as QAE, so long as they can adequately handle the inclement noise. Removing the need to perform depth-consuming arithmetic operations on the quantum computer, as we have done in this paper, is an important step towards this.
Will Fourier QMCI still be useful beyond the NISQ era?
One of the most enduring ideas in quantum algorithm design is that, because quantum computers will not offer useful computational speed-ups for all problems, the focus should therefore be on using the quantum computer only when it offers a computational advantage. Moreover, even within a single task those parts of the computation that can be handed off to a classical computer without adversely impacting overall performance should be. This essential paradigm is unlikely to be restricted solely to the NISQ era, especially as the transition from the NISQ era to FSFT era is unlikely to be particularly sharp (as outlined above). In this paper, we have shown how the task of Monte Carlo integration can be sped-up by using a quantum computer to perform amplitude estimation (where there is a quantum advantage), but not arithmetic (where there is no quantum advantage). We therefore expect that the techniques we have proposed will find useful application for the foreseeable future.
A Proof of Corollary 4 Corollary 4. The quantity, µ = E(f (X) g(Y )), where (X, Y ) ∼ p(x (i) , x (j) ) and f and g are functions that are continuous in value and first derivative and whose second and third derivatives are piecewisecontinuous and bounded, can be estimated with MSE ∈ Θ(q −λ ), where q is the number of uses of a circuit preparing the quantum state |p (ie, as defined in (2)) and λ is the convergence rate of some QAE subroutine (ie, as defined in (7)) which operates on circuits of the form defined in Fig. 4.