Quantum mean estimation for lattice field theory

TL;DR

This paper applies quantum mean estimation to lattice field theories, achieving quadratic speedup over classical Monte Carlo methods, and explores its robustness against sign problems and errors in quantum gates.

Contribution

It demonstrates the effectiveness of quantum mean estimation in lattice field theories, including models with sign problems and error analysis for fault-tolerant quantum computing.

Findings

Quadratic advantage over Monte Carlo methods.

Effective in presence of sign problems and critical slowing down.

Analyzed impact of gate synthesis errors on quantum algorithms.

Abstract

We demonstrate the quantum mean estimation algorithm on Euclidean lattice field theories. This shows a quadratic advantage over Monte Carlo methods which persists even in presence of a sign problem, and is insensitive to critical slowing down. The algorithm is used to compute $\pi$ with and without a sign problem, a toy U(1) gauge theory model, and the Ising model. The effect of $R_{Z}$-gate synthesis errors on a future fault-tolerant quantum computer is investigated.

Figures

Equations

$$⟨\mathcal{O}⟩=\frac{1}{Z}\sum _i\mathcal{O}(g_i)e^{-S(g_i)}.$$

$$U_{\phi }\equiv \sum _i{e}^{i\phi (g_i)}|g_i⟩⟨g_i|,$$

$$\Re [⟨{\psi }_0|U_{\phi }|{\psi }_0⟩]$$

$$=\sum _i|⟨{\psi }_0⟩g_i{|}^2\Re [e^{i\phi (g_i)}]$$

$$=\sum _i|⟨{\psi }_0⟩g_i{|}^2\cos (\phi (g_i)),$$

$$U_{\mathcal{O}}|g_i⟩|0⟩=|g_i⟩|\phi (g_i)⟩.$$

$$\text{mathbbm}1\otimes e^{i\hat{\phi }}|g_i⟩|\phi (g_i)⟩=e^{i\phi (g_i)}|g_i⟩|\phi (g_i)⟩,$$

$$e^{i\phi (g_i)}{U}_{\mathcal{O}}^{†}|g_i⟩|\phi (g_i)⟩=e^{i\phi (g_i)}|g_i⟩|0⟩.$$

$$|\psi ⟩$$

$$|\psi ⟩=\frac{1}{\sqrt{2M}}\left(|0⟩\sum _i{e}^{i{\varphi }_i}|g_i⟩+|1⟩\sum _i{e}^{-i{\varphi }_i}|g_i⟩\right).$$

$$|\psi ⟩=\sqrt{\frac{Z}{M}}|{\psi }_0⟩+\sqrt{\frac{M-Z}{M}}|{\psi }_1⟩$$

$$|{\psi }_0⟩$$

$$=|0⟩\sum _i\sqrt{\frac{w_i}{Z}}|g_i⟩$$

$$|{\psi }_1⟩$$

$$=i|1⟩\sum _i\sqrt{\frac{1-w_i}{M-Z}}|g_i⟩.$$

$$Z=\sum _{i}w(g_i)>\sum _i{w}_{\text{min}}=M{w}_{\text{min}}>0.$$

$$|{\psi }_0⟩=|0⟩\sum _i\sqrt{\frac{w_i}{Z}}|g_i⟩.$$

$$|{\psi }_0⟩=|+⟩\sum _i\sqrt{\frac{w_i}{Z}}|g_i⟩.$$

$$⟨\mathcal{O}⟩$$

$$=\frac{{\sum }_i\frac{\mathcal{O}(g_i)e^{-S(g_i)}}{q(g_i)}q(g_i)}{{\sum }_{i}q(g_i)}\frac{{\sum }_{i}q(g_i)}{{\sum }_i\frac{e^{-S(g_i)}}{q(g_i)}q(g_i)}$$

$$=\frac{\frac{1}{Z_q}{\sum }_i\tilde{\mathcal{O}}(g_i)q(g_i)}{\frac{1}{Z_q}{\sum }_i\tilde{\mathcal{R}}q(g_i)}=\frac{{⟨\tilde{\mathcal{O}}⟩}_q}{{⟨\tilde{\mathcal{R}}⟩}_q}$$

$$\mathcal{Q}=\sum _q{e}^{i\chi (q)}|q⟩⟨q|.$$

$$\mathcal{Q}|{\psi }_0⟩=\sum _q{e}^{i\chi (q)}⟨q⟩{\psi }_0|q⟩,$$

$$\text{mathbbm}{U}_{\phi }=|0⟩⟨0|\otimes U_{\phi }+|1⟩⟨1|\otimes U_{\phi }^{†},$$

$$\cos (\theta )$$

$$|{\psi }_0⟩=\frac{1}{\sqrt{2}}(|{\eta }_{+}⟩-|{\eta }_{-}⟩)$$

$$\frac{\pi }{4}=\frac{1}{Z}{\iint }_0^{1}dxdy\Theta (1-r^2(x,y))$$

$$\frac{\pi }{4}$$

$$(\mathfrak{add}|x⟩|0⟩)(\mathfrak{add}|y⟩|0⟩)=|x⟩|x⟩|y⟩|y⟩.$$

$$(\mathfrak{mul}|x⟩|x⟩)(\mathfrak{mul}|y⟩|y⟩)=|x⟩|x^2⟩|y⟩|y^2⟩,$$

$$(\mathfrak{add}|x^2⟩|y^2⟩)(|x⟩|y⟩)=|x^2⟩|x^2+y^2⟩|x⟩|y⟩.$$

$$\mathfrak{comp}|x^2+y^2⟩|(2^n-1{)}^2⟩=|x^2+y^2⟩|r^2\le (2^n-1{)}^2⟩,$$

$$\frac{\pi }{12}$$

Full text

††preprint: FERMILAB-PUB-23-095-QIS-T

I Introduction

Computations in lattice field theory (LFT) are typically framed within the formalism of statistical mechanics [1, 2, 3]. This is possible by analytically continuing the field theory from Minkowski to Euclidean spacetime, which transforms complex phases into probability weights. This changes the quantum path integral into a classical partition function, $Z=\sum_{i}e^{-S(g_{i})}$, which is a sum over all configurations of the degrees of freedom, $g_{i}$, weighted by the lattice action, $S\geq 0$. Quantum expectation values of observables, $\mathcal{O}$, then transform into statistical averages

[TABLE]

Summing over all configurations is, in general, impractical. Instead, Monte Carlo methods sample a finite number, $N$, of them with an error, $\sigma$, on $\langle\mathcal{O}\rangle$ scaling asymptotically like $\sigma\sim 1/\sqrt{N}$. Often, only a few thousand samples yield precision predictions for complex observables like hadronic form factors [4] or QCD contributions to the muon $g-2$ [5].

However, there are limitations. These include situations where there is a sign problem [6, 7], and critical slowing down [8]. The sign problem crops up at finite fermion density [9] and when simulating real-time dynamics [10]. Critical slowing down is a consequence of running Monte Carlo algorithms with updates that neglect long-distance correlations. Such a situation appears as one tries to approach the continuum limit or when studying topological observables [11]. In the case of the sign problem, the required $N$ scales exponentially in model parameters, and with critical slowing down it scales as a power. These obstacles have prompted new classical algorithms that attempt to address these issues including: cluster algorithms [12, 13], dual variables approaches [14, 15], tensor networks [16, 17], complexification [10], and density of states methods [18]. Despite these successes, it is unlikely that general solutions exist cf. [19].

Through quantum computers it is possible to avoid some of these limitations. This result comes fundamentally from the abilities of quantum computers to enumerate an exponential number of states via quantum superposition, and their capacity to generate entanglement. Notable quantum algorithms which harness these properties include the quantum Fourier transform [20, 21, 22], quantum phase estimation (QPE) [23, 24, 25, 26], and Hamiltonian simulation. This last class of algorithms could further allow for tremendous advances in simulating real-time dynamics for LFT [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 44, 68, 74, 75, 76, 77, 41, 78, 79, 80, 81, 82, 83]. Investigation of quantum algorithms have also lead to faster classical algorithms [84, 85].

Another quantum algorithm relevant specifically to this work is quantum mean estimation (QME) [86, 87, 88]. Quantum mean estimation is capable of quadratically reducing the asymptotic scaling of $\sigma$ to $\sim 1/N$. It achieves this using QPE, and by using superposition to incorporate a full probability distribution into calculations of $\braket{\mathcal{O}}$. It has been developed for boolean variables, positive bounded real variables [89, 90], bounded real variables [87], and unbounded real variables [91, 86, 88].

In this article, we will demonstrate how and when QME can be used to improve classical LFT calculations that use Monte Carlo sampling. In Sec. II, we detail how to use QME on classical statistical models, and provide circuits to construct the QME algorithm. In Sec. III, we provide numerical results including circumstances with sign problems, and investigate the effects of noise. In Sec. IV, we compare traditional sampling methods to the QME algorithm. We conclude with a discussion and directions for future work in Sec. V.

II Quantum Mean Estimation

In this section, we will demonstrate how a quantum computer using QME can provide an estimate of $\langle\mathcal{O}\rangle\in[-1,1]$ with fixed precision and quadratically fewer resources than traditional sampling using the method of Ref. [87].

II.1 Overview of QME

Estimating $\braket{\mathcal{O}}$ on a quantum computer uses QPE as a subroutine. By encoding $\mathcal{O}(g_{i})$ into the phases of states and running QPE on a judiciously chosen unitary and starting state, a phase is returned approximating the mean. Alternatively, the mean can be stored as the amplitude of some target state. This amplitude can be approximated and returned—again using QPE—using quantum amplitude estimation [92, 90, 89, 91, 88].

Following Ref. [87], for a given unitary—or phase oracle—in diagonal form,

[TABLE]

QME returns

[TABLE]

where $\ket{\psi_{0}}$ is a generic initial state.111One can also estimate $\bra{\psi_{0}}U_{\varphi}\ket{\psi_{0}}$, but we focus on the real part. For appropriate $\varphi$s, one can calculate averages of real numbers $\in[-1,1]$. We provide here a method to construct $U_{\varphi}$ through qubit arithmetic, by computing the appropriate $\varphi_{i}$ into a register, applying a phase based on that register’s value, and then uncomputing the phase register back to $|0\rangle$. The computation of the $\varphi_{i}$ is done by a new oracle, $U_{\mathcal{O}}$ which for a given configuration, $g_{i}$, places the state $\ket{\varphi}$ into a separate register. In terms of states and operators,

[TABLE]

The new “working” register stores the value of the appropriate phases to apply to each configuration. From that state, the proper phase can be applied via

[TABLE]

and the working register can be uncomputed,

[TABLE]

This entire circuit can be captured by a single unitary $U_{\varphi}=U_{\mathcal{O}}^{\dagger}(\mathbbm{1}\otimes e^{i\hat{\varphi}})U_{\mathcal{O}}$, as seen in Fig. 1.

Using QME corresponds to recasting Eq. (1) in the form of Eq. (II.1). The relevant identification is to see that the product $|\braket{\psi_{0}}{g_{i}}|^{2}\cos(\varphi(g_{i}))\sim(e^{-S(g_{i})}/Z)\mathcal{O}(g_{i})$. This identification in turn allows for flexibility in defining $\ket{\psi_{0}}$ and the phases $\varphi(g_{i})$. Since we can identify $|\braket{\psi_{0}}{g_{i}}|^{2}$ in Eq. (II.1) as the probability for configuration $g_{i}$ used in the average, a natural correspondence is given by $|\braket{\psi_{0}}{g_{i}}|^{2}\sim e^{-S(g_{i})}/Z$ as well as $\cos(\varphi(g_{i}))\sim\mathcal{O}(g_{i})$, which can be used to define $\varphi(g_{i})$. However, one is allowed the freedom to reweigh parts of the observables into the weights, and vice versa, to optimize. By rescaling, this method can compute averages in any finite range.

To actually compute Eq. (II.1) requires preparing a state $|\psi_{0}\rangle$, and the use of QPE with $U_{\varphi}$. We will present two separate implementations of QME using two different $|\psi_{0}\rangle$. One directly prepares the Boltzmann weights, and the other shifts the difficulties of preparing such states into the computation of observables through reweighting.

II.2 Using state preparation

We can use any state-preparation method, whether it be through efficient classical methods [93], black-box methods [94, 95, 96], using the quantum singular value transform [97], or with state search [98] to construct an appropriate $\ket{\psi_{0}}$. Here, we prepare $\ket{\psi_{0}}$ using arithmetical oracles such that $|\braket{\psi_{0}}{g_{i}}|^{2}=e^{-S(g_{i})}/Z$. We will denote this state preparation method as $\mathcal{A}^{\text{SP}}$.

To prepare a quantum state with the desired probability distribution, let $0<w(g_{i})\leq 1$ be the weights of the system, and $Z=\sum_{i}^{M}w(g_{i})$ with $M$ the total number of configurations. We can prepare a maximal superposition over all configurations, along with an ancilla in the state $|+\rangle$:

[TABLE]

Next, we apply a controlled-$U_{\phi}$ with $\phi_{i}=\arccos(\sqrt{w_{i}})$ using the procedure from Sec. II.1,

[TABLE]

The phases used in state preparation are denoted with $\phi$, and those used in encoding observables with $\varphi$, since both use the same phase oracle. After applying a Hadamard to the ancilla we can rewrite the whole state as,

[TABLE]

where $\ket{\psi_{0}}$ is the desired state,

[TABLE]

and $\ket{\psi_{1}}$ is an orthogonal state,

[TABLE]

$\ket{\psi_{0}}$ contains the target probability distribution; however, $\ket{\psi_{1}}$ needs to be removed. We use fixed-point oblivious amplitude amplification [99, 100, 101], which requires a lower-bound on the amplitude, to isolate the zero-ancilla state. A crude lower bound is given by assuming every $g_{i}$ comes with the smallest possible weight, then,

[TABLE]

We find that $\sqrt{Z/M}>\sqrt{w_{\text{min}}}$, and we can use $\sqrt{w_{\text{min}}}$ as a lower bound for amplitude amplification. Asymptotically the algorithm will take $O(\log(1/\epsilon)/\sqrt{w_{\text{min}}})$ amount of time to achieve a desired $\epsilon$ accuracy [99]. Deriving a tighter lower bound would reduce the algorithmic time; however, since the actual amplitude is $\sqrt{Z/M}$, at best we can expect $O(\sqrt{M}\log(1/\epsilon)/\sqrt{Z})$ time. The algorithm requires no conditional measurements, or “repeat-until-success” steps, and is captured by a unitary matrix. After amplitude amplification we end with the desired state,

[TABLE]

The QME algorithm requires an ancilla in the state $|+\rangle$, so after applying a Hadamard gate, the initial state is

[TABLE]

Having prepared $|\psi_{0}\rangle$, we can use QPE to perform the mean estimation, but this state preparation is potentially expensive; therefore, we discuss an alternative approach.

II.3 Using reweighting

Previously, we discussed how there exists freedom in associating the product $\mathcal{O}(g_{i})e^{-S(g_{i})}/Z$ to $|\braket{\psi_{0}}{g_{i}}|^{2}$ and $\cos(\varphi(g_{i}))$. The choices made determine the state preparation $\mathcal{A}$ and the specific means computed. This reweighting [102] can be understood by introducing into Eq. (1) a second probability distribution $q(g_{i})$:

[TABLE]

where $Z_{q}=\sum_{i}q(g_{i})$ and $\braket{\tilde{\mathcal{R}}}$ is a necessary reweighting factor. So instead of $\ket{\psi_{0}}$ with the distribution $e^{-S(g_{i})}$, we can construct a different one by shuffling part of the distribution into the observable, and then computing the ratio of two means. One nice choice is the uniform distribution $q(g_{i})=1/M$ which can be prepared by $\mathcal{A}^{\text{RW}}=H^{\otimes}$. With two different state preparation methods described, we now discuss the role of QPE.

II.4 Quantum phase estimation

The QPE algorithm returns an estimate of a phase $\chi$ of an eigenvector $\ket{q}$ for a unitary $\mathcal{Q}$ such that $\mathcal{Q}\ket{q}=e^{i\chi}\ket{q}$. Thus, QPE expects an input state, or a state preparation method $\mathcal{A}$, a unitary matrix $\mathcal{Q}$ whose eigenvalues are of interest, and a result register of $N_{r}$ qubits prepared in a maximal superposition. It then returns the best, approximate phase with a probability $\geq 4/\pi^{2}$ [20, 26].

We can see the effect of the input state not being an eigenstate of $\mathcal{Q}$ by inspecting its spectral decomposition,

[TABLE]

For a generic state $\ket{\psi_{0}}$,

[TABLE]

and QPE operating on $\ket{\psi_{0}}$ returns angle $\chi(q)$ for each eigenstate into the result register in superposition with probability $|\braket{q}{\psi_{0}}|^{2}$. The accuracy of the returned phase and how likely that phase is to be measured, is specified by $N_{r}$. The possible phase values are restricted to the possible binary fractions expressible by $N_{r}$ qubits.

II.5 Synthesis of parts

Here we relate $\mathcal{A}$ to $\mathcal{Q}$ to see how QME emerges from QPE. In Secs. II.2 and II.3 we discussed two $\mathcal{A}$ methods. With those in mind, we will now construct $\mathcal{Q}$ used in QME. The unitary $\mathcal{Q}$ is given by the product $\mathcal{Q}\equiv S_{\mathcal{A}}\mathbbm{U}_{\varphi}$ with a phase oracle, $\mathbbm{U}_{\varphi}$, and the operator, $S_{\mathcal{A}}$. The form of $S_{\mathcal{A}}$ is: $S_{\mathcal{A}}\equiv\mathcal{A}S_{0}\mathcal{A}^{\dagger}$, with $S_{0}=2\ket{0}\bra{0}-\mathbbm{1}$. The circuit for $-S_{0}$ can be seen in Fig. 2, and consists of applying a minus sign to the “all zeros” state. The phase oracle of Ref. [87] has two terms,

[TABLE]

and is responsible for applying the phases encoding $\mathcal{O}$ onto the states. The oracle $U_{\varphi}$ is described in Sec. II.1. The circuit for $\mathbbm{U}_{\varphi}$ can be seen in Fig. 3, and with it $\mathcal{Q}$ can be completely defined. The circuit for $\mathcal{Q}$ can be seen in Fig. 4. The entire circuit for QPE—and through it QME—can be seen in Fig. 5.

With the relationship between $\mathcal{A}$ and $\mathcal{Q}$ established the desired mean must come about from an eigenvalue of $\mathcal{Q}$. To see this, let us revisit Eq. (II.1), and define the angle $\theta$ along with $\cos(\theta)$ through,

[TABLE]

Written in this form we see that Eq. (19) computes the average value of cosine. Now, consider the initial state $\ket{\psi_{0}}\equiv\mathcal{A}\ket{0}$. While not an eigenstate of $\mathcal{Q}$, it is a simple linear combination of eigenstates of $\mathcal{Q}$ (see Ref. [87]),

[TABLE]

with $\mathcal{Q}\ket{\eta_{\pm}}=e^{\pm i\theta}\ket{\eta_{\pm}}$. Since $\ket{\psi_{0}}$ is constructed from eigenvectors of $\mathcal{Q}$, applying QPE to $\ket{\psi_{0}}$ returns either of the two phases associated with those states, $\theta$ or $2\pi-\theta$, with equal probability. Then by taking the cosine of the output angle—which is insensitive to $\theta$ or $2\pi-\theta$—we recover the average in Eq. (19).

To summarize,

Using a qubit encoding of $g_{i}$, construct $\mathcal{A}$. This can either be done by making a state with the correct probability distribution as in Sec. II.2, or by reweighting as in Sec. II.3.

• 2)

Construct $\mathbbm{U}_{\varphi}$:

• –

Construct $U_{\mathcal{O}}$. This is model- and observable-dependent but consists of qubit arithmetic. Define the phases, $\varphi(g_{i})$, as needed.

• –

Construct $U_{\varphi}$ from $U_{\mathcal{O}}$ and $e^{i\hat{\varphi}}$ using Fig. 1.

• –

Construct $\mathbbm{U}_{\varphi}$ from $U_{\varphi}$ using Fig. 3.

• 3)

Construct $\mathcal{Q}$ from $\mathcal{A}$, $\mathbbm{U}_{\varphi}$, and $S_{0}$ using Fig. 4.

• 4)

Perform QPE with $\mathcal{A}$ and $\mathcal{Q}$ as in Fig. 5. The output is an estimated angle, $\omega$, where $\cos(\omega)\approx\sum_{i}|\braket{\psi_{0}}{g_{i}}|^{2}\cos(\varphi(g_{i}))$.

With the method now explained, we proceed with some examples to help solidify each step above.

III Applications

We consider several numerical examples: computing $\pi$ with and without negative weights, a toy U(1) gauge theory, and the two-dimensional Ising model. Then, a general QME is described for lattice gauge theories. Finally, how $R_{z}$-gate synthesis errors affect QME is discussed. All of the numerical results in this section were performed on the IBM Qiskit QASM simulator using 1024 shots.

III.1 $\pi$ with and without a sign problem

In the Monte Carlo estimate of $\pi$, one randomly chooses points, $x_{g},y_{g}$ in the $x$-$y$ plane—which we assume to be a grid (see Fig. 6)—such that $x_{g},y_{g}\in[0,1)$. A quarter circle in that plane has an area of $\pi/4$, while the square has area one. Then we expect that on average the fraction of points that fall within the boundary of the circle is $\approx\pi/4$. We can decide this for each $x_{g}$-$y_{g}$ pair if $x_{g}^{2}+y_{g}^{2}\leq 1$, and by counting the number of points that fall inside the circle we can compute $\pi$.

To compute $\pi$ using a quantum computer, parameterize $x$ and $y$ using $N_{g}=2n$ qubits: $\ket{x}=H^{\otimes n}\ket{0}$, $\ket{y}=H^{\otimes n}\ket{0}$, which creates an equal superposition of all possible $x$-$y$ pairs. The $x$ ($y$) coordinate is defined as the integer representation of the bit-string associated with $\ket{x}$ ($\ket{y}$) normalized by $2^{n}-1$. For example, the state $\ket{111}\ket{010}=\ket{7}\ket{2}$ when normalized corresponds to $\ket{7/7}\ket{2/7}=\ket{1}\ket{0.2857..}$. With the addition of a single ancilla qubit prepared in the $\ket{+}$ state, this preparation of the $x$-$y$ coordinates constitutes $\mathcal{A}=\mathcal{A}^{\text{SP}}=\mathcal{A}^{\text{RW}}$.

To construct $\mathcal{Q}$, we need to define the phase oracle, $\mathbbm{U}_{\varphi}$. Let us consider a statistical mechanics average:

[TABLE]

where $\Theta$ is the Heaviside function, $r^{2}(x,y)=x^{2}+y^{2}$, and $Z\equiv\int_{0}^{1}dx\int_{0}^{1}dy$. Going forward we will give the $x$ and $y$ dependence once, but omit it afterward. We approximate this integral using a grid of points,

[TABLE]

where $Z=2^{N_{g}}$, i.e. the number of points in the grid.

The phases for $U_{\mathcal{O}}$ are then defined as $\varphi(x,y)=\arccos(\Theta(1-r^{2}))$, with the two cases: $\arccos(0)=\pi/2$ when $r^{2}>1$, and $\arccos(1)=0$ when $r^{2}\leq 1$. Since the $x$-$y$ coordinates are in superposition we can demonstrate the action of $U_{\mathcal{O}}$ on a single pair without loss of generality. We will use quantum arithmetic gates [103, 104] specifically an addition gate, $\mathfrak{add}$, and a multiplication gate, $\mathfrak{mul}$. Further we require a comparator $\mathfrak{comp}$, that returns a boolean variable [105]. First we add $x$ and $y$ into two registers,

[TABLE]

Then, we multiply the registers together,

[TABLE]

and add the products together

[TABLE]

Finally, we compare to a register$\ket{(2^{n}-1)^{2}}$,

[TABLE]

where the comparison register is either zero or one. This defines $U_{\mathcal{O}}$. Using this register, a phase can be applied which gives $e^{i\pi/2}$ if the comparison bit is zero, or $1$ if it is one. The working registers can be uncomputed, defining the entire unitary $U_{\varphi}$. With $U_{\varphi}$ defined, $\mathbbm{U}_{\varphi}$ is defined, and with $\mathcal{A}$, $\mathcal{Q}$ is defined, and QPE can be performed.

Histograms showing counts of bit-strings using different $N_{g}$ and $N_{r}$ are in Figs. 7, and 8, respectively. Along with the histograms are lines for the finite-grid ratios of inner-circle points to total points, as well the exact value of $\pi/4$. In Fig. 7, we see an increase in $N_{g}$ at fixed $N_{r}=8$, leading towards $\pi/4$, with a constant bin width at no less than an exponential rate. In practice, one could extrapolate the measured value as a function of $N_{g}$ with potentially fewer resources than a single fine grid might require. In Fig. 8, we see varying $N_{r}$ while keeping $N_{g}=8$ yields a most likely value at the finite-grid ratio, and the bin width diminishes as $2^{-N_{r}}$. Thus, larger $N_{g}$ improves accuracy, while larger $N_{r}$ improves precision.

We can include a sign problem by weighing locations in the quarter circle using positive and negative weights separated at $\theta=\pi/3$ which corresponds to $(x/r)^{2}=1/4$ (See Fig. 6 (d)). When $(x/r)^{2}\leq 1/4$, the weight is $-1$, and when $(x/r)^{2}\geq 1/4$, the weight is $+1$. This corresponds to the integral

[TABLE]

To estimate $\pi$ from Eq. (27) we must modify $\mathbbm{U}_{\varphi}$ to also check if $x^{2}/r^{2}\geq 1/4$. We compute $\ket{x^{2}}$ and $\ket{r^{2}}$ as before, and using a division oracle calculate $\ket{x^{2}/r^{2}}$ which is then compared using $\mathfrak{comp}$ with $\ket{1/4}$. We then apply a phase based on the compare register. Therefore, with a sign problem there are two deciding qubits, $\ket{q_{\text{circle}}}$ and $\ket{q_{\text{angle}}}$, which indicate if the point is inside the circle and $(x/r)^{2}\geq 1/4$, respectively. From their values a phase $e^{i\varphi}$ is imparted, with $\varphi=\pi(1+q_{\text{circle}}-2q_{\text{angle}})/2$. With the working register arithmetic outlined above, $U_{\mathcal{O}}$, $U_{\varphi}$, and $\mathbbm{U}_{\varphi}$ are defined.

Analogous histograms to the case of $\pi$ without a sign problem can be seen in Figs. 9 and 10. We see in Fig. 9 as $N_{g}$ is varied with fixed $N_{r}$, the most likely value moves towards $\pi/12$ for larger values of $N_{g}$. Likewise, in Fig. 10 as $N_{r}$ is increased we see the distribution narrows around the finite-grid ratio. For QME this problem is as computationally difficult as without a sign problem. One should note that the nonlinear transform from the measured $\omega$ to $\cos(\omega)$ results in non-constant bin width. This is acutely important for sign problems and reweighting when means are often near zero, where the bins are largest. Thus rescaling the result register to optimize bin locations could be valuable for reducing costs.

III.2 U(1) gauge theory toy model

We introduce a toy model resembling the Villian [106, 107] approximation of a two-dimensional, Euclidean U(1) gauge theory, and the calculation of the Wilson loop. This has been used to investigate classical methods for ameliorating the sign problem [108]. We aim to calculate

[TABLE]

Here, $n$ controls the degree of oscillations in the integrand and can be any real number. The value of $\mathcal{P}(n)$ is

[TABLE]

We again approximate $x$ and $y$ by a grid and $Z=2^{N_{g}}$.

State preparation is given by $\mathcal{A}=\mathcal{A}^{\text{SP}}=\mathcal{A}^{\text{RW}}$. We will encapsulate $\mathcal{O}(g_{i})$ inside the phases and define $\varphi=\arccos\left(\Theta(1-r^{2})\cos(n\theta)e^{-\theta^{2}}\right)$. Designing the circuit for computing the phases requires floating-point arithmetic, which we will not elaborate upon, but will assume exist. With the ability to create $\ket{\varphi}$, we can apply $\exp(i\hat{\varphi})$ to the phase register and extract the phase. After uncomputing the phase register, this procedure defines $U_{\varphi}$ and $\mathbbm{U}_{\varphi}$.

Results of these circuits for $\mathcal{P}(3)$ are shown in Figs. 11 and 12. In Fig. 11 we see that as $N_{r}$ is increased the distribution narrows around the exact value. Figure 12 shows the drift of the most likely value of $\mathcal{P}(3)$ towards the exact value of $\approx 0.041$ as $N_{g}$ increases at fixed $N_{r}$.

Having given circuits and explained QME for a variety of toy and demonstrative models, we now turn to an archetypal LFT, the Ising model.

III.3 The Ising model

We consider a two-dimensional Ising model on a $V=N_{s}\times N_{t}$ Euclidean lattice with periodic boundary conditions and spatial and temporal extents of $N_{s}$ and $N_{t}$. The action is

[TABLE]

where $\sum_{\langle ij\rangle}$ sums over nearest-neighbor pairs, $s_{i}$ is a lattice spin taking values $\pm 1$, and $\beta$ is a coupling constant commonly identified as an inverse temperature. We can rewrite this as a sum over binary values, $n_{i}=0,1$, using $s_{i}=2n_{i}-1$. Then

[TABLE]

As an observable, we will compute the square of the magnetization density

[TABLE]

which is bounded between zero and one.

On the quantum computer each qubit in the system is associated with a spin on the lattice. To prepare the system, we put every qubit into the $\ket{+}$ state, $\ket{\{n\}}\equiv H^{\otimes V}\ket{0}$. This generates all spin configurations in superposition. By appending one ancilla prepared in the $\ket{+}$ state we have the state preparation circuit $\mathcal{A}^{\text{RW}}$, but must calculate a ratio of observables (see Sec. II.3).

Just as in the $\pi$ example, we require a phase oracle. We will focus on a single spin configuration without loss of generality. To compute the magnetization, we must sum the bit-string. We do this by first summing every even-odd spin together in parallel. We can then add the results of the first step again in even-odd pairs. Repeating this we can compute the sum in $\sim\log_{2}(V)$ steps. The result of these sums is $\ket{\sum_{i}n_{i}}$. By multiplying by two, dividing by the volume, and shifting by one we arrive at the magnetization density. We can add this number to zero, and multiply the two $m$ registers to obtain $m^{2}$.

To compute the action, we already have the term containing $\sum_{i}n_{i}$. To compute the nearest-neighbor interaction we can again compute the product of disjoint pairs of spins in each of the two directions in parallel, and sum the resulting collection of products. This can be done in $\sim\log_{2}(2V)$ steps, resulting in a register $\ket{\sum_{\langle ij\rangle}n_{i}n_{j}}$. Multiplying by the appropriate factors and subtracting we can form the register $\ket{S}$. Then using exponentiation, multiplication, and taking the arccosine we arrive at $\ket{\arccos(m^{2}e^{-S})}$. This defines $U_{\mathcal{O}}$ for the numerator. $U_{\mathcal{O}}$ for the denominator can be calculated similarly without the factor of $m^{2}$. Now by applying $\mathbbm{1}\otimes e^{i\hat{\varphi}}$ we can extract the phase, then uncompute the phase register. These steps define the unitary $U_{\varphi}$, and hence $\mathcal{Q}$, in QPE.

In Figs. 13 and 14 we see histograms of $\langle m^{2}\rangle$ when varying $N_{g}$ and $N_{r}$ respectively, for the one-dimensional, classical Ising chain at $\beta=0.4$. As $N_{g}$ is increased for fixed $N_{r}$, the mean tends towards the infinite-volume value. Likewise, with $N_{g}$ fixed and $N_{r}$ increased, estimates converge to the exact finite-volume value.

Fig. 14 demonstrates a desirable feature that results from 1) the bin-width providing a 100% confidence interval of error, and 2) the maximal bin changing location. That is that the ideal result must lie in the overlap between distributions. Therefore, the overlap between bins of different $N_{r}$ can estimate the mean more precisely than any single $N_{r}$ alone. This could be used to engineer sets of $N_{r}$ that maximize information gained. In Fig. 14 we see a sliver of overlap between the $N_{r}=5$ and $6$ maximal bins, which the exact answer must lie within. Fig. 11 provides another example. Having now seen QME used with simple examples, as well as the formulation and execution on a bona fide LFT, we now give the formulation of QME for a generic lattice gauge theory.

III.4 Generic lattice gauge theory

Take a $D$ dimensional lattice with $V\equiv N_{s}^{D-1}N_{t}$ and periodic boundary conditions. We consider the action

[TABLE]

where $g_{\mu\nu}(x)=g_{\mu}(x)g_{\nu}(x+\hat{\mu})g^{\dagger}_{\mu}(x+\hat{\nu})g^{\dagger}_{\nu}(x)$ is a plaquette, $g_{\mu}(x)$ is a gauge group element, and the sums are over lattice sites and the two directions $\mu$ and $\nu$ respectively .

Prepare a qubit register for each gauge link, $\ket{g_{\mu}(x)}$. Let the operator $\mathcal{H}$ prepare a maximal superposition of all possible gauge link states, analogous to the Hadamard gate,

[TABLE]

where $|G|$ is the size of the local state space. Then for state preparation we can apply $\mathcal{H}$ for every link register,

[TABLE]

This creates a superposition of all field configurations with equal probability. We then append an ancilla qubit $\ket{+}$. This constructs $\mathcal{A}^{\text{RW}}$. For $\mathcal{A}^{\text{SP}}$ we follow Sec. II.2 with $\phi(g_{i})=\arccos(e^{-S(g_{i})/2})$. When encoding observables one sets $\varphi(g_{i})=\arccos(\mathcal{O}(g_{i}))$ in $U_{\varphi}$. With these prescriptions of the phases, the operators for QME of lattice gauge theory are defined. In the next section we will study the effects of noise on QME.

III.5 The effects of noise

On a fault-tolerant quantum computer the $R_{Z}$ gates are approximated with an infidelity $\varepsilon$ by interleaved $T$ and $H$ gates [109, 110, 111, 112, 113, 26, 114, 115]. While one can derive a string of $T$ and $H$ gates that approximates any $R_{Z}$ gate, implementing this string and classically simulating it drastically extends the length of a quantum circuit, and becomes a computationally intensive problem. Instead, to study the effects of nonzero $\varepsilon$, we approximate each $R_{Z}$ gate,

[TABLE]

The parameter $\varepsilon$ approximates the infidelity in quantum computations, and drives coherent angle-dependent drift in the $R_{Z}$ gates. We tested variations on Eq. (36) by sending $\theta\varepsilon\rightarrow-\theta\varepsilon$ in $R_{X}$ and found negligible effects.

For a test case we consider the model from Sec. III.2. Using $N_{r}=N_{g}=6$ we run the QME circuit for various $\varepsilon$. This calculation requires synthesizing $N_{Z}\sim O(10^{4})$ $R_{Z}$ gates. Fig. 15 shows the histogram of the average returned from QME. For $\varepsilon\approx 10^{-4}$ a statistically clear signal for the most likely value consistent with the ideal value appears with $N_{\text{shots}}=1024$.

In context, this indicates there is a threshold $\varepsilon$ for a given set of shots below which better synthesis of $R_{Z}$ provides no benefit at fixed $N_{r}$ and $N_{g}$. Exact rotation gates are unnecessary, and imprecise rotation gates can give sufficient accuracy and precision. This is analogous to precision LFT calculations where half-precision floating point numbers are used to accelerate calculations [116, 117, 118, 119]. Heuristically, the threshold is $\varepsilon\sim\sqrt{N_{\text{shots}}}/N_{Z}$.

IV Quantum advantage

To appreciate the advantage provided by the QME algorithm, consider $\mathcal{A}^{\text{SP}}\ket{0}=\ket{\psi_{0}}$. At this stage, traditional sampling methods could be applied, and measuring $\ket{g_{i}}$ provides a single configuration. By iterating this procedure one collects an ensemble of configurations from which averages can be taken. We expect $\sigma\sim 1/\sqrt{N}$, but $N$ is precisely the number of times $\mathcal{A}^{\text{SP}}$ is run. In this sense we can see $\mathcal{A}^{\text{SP}}$ as the quantum mechanical analog to having access to the partition function through, say, brute-force calculation.

Now, without loss of generality, take $N$ to be a power of two, $N=2^{k}$, implying $\sigma\sim 1/2^{k/2}$. For the QME algorithm, we want the same precision of $1/2^{k/2}$. Using QME we see this requires $N_{r}=k/2$ qubits. Looking at Fig. 5, QPE executes $\mathcal{Q}$ a total of $\sum_{j=0}^{k/2-1}2^{j}=2^{k/2}-1$ times. Since $\mathcal{A}^{\text{SP}}$ appears in $\mathcal{Q}$ ${O}(1)$ times, $\mathcal{A}^{\text{SP}}$ is similarly called ${O}(2^{k/2})$ times. Thus we see the number of times $\mathcal{A}^{\text{SP}}$ is needed in the QME algorithm ($2^{k/2}$) is quadratically less compared to the traditional sampling method ($2^{k}$). Therefore rather than using the prepared state to sample and calculate expectation values, it is advantageous to pass the quantum state into the QME algorithm to avoid greater calls to the state preparation algorithm overall, while maintaining the same precision.

This increased efficiency in calls to $\mathcal{A}^{\text{SP}}$ is seen in Fig. 16, where QME is compared with the traditional sampling method to calculate $\pi$ for $N_{g}=2$ and varying $N_{r}$. Both the QME and sampling circuits are decomposed under the same basis gates. We see the overhead for QME is large, however for a fixed accuracy, the scaling of QME is superior to sampling. The ratio between the sampling and QME slopes is $\approx 2.09$, approaching the asymptotic quadratic value. The sampling circuit depth is computed using the circuit depth for $\mathcal{A}$, multiplied by the number of samples required to achieve a precision of $2^{-N_{r}}$ assuming the error converges like $2^{-N_{r}/2}$. After $N_{r}=8$, QME surpasses the traditional sampling method in efficiency.

What about in the case of a sign problem, or an exponentially small value to calculate? This is relevant in the case of $\mathcal{A}^{\text{RW}}$. Consider the case where the degrees of freedom of the model are binary variables themselves. Now the prefactor $Z/M$ in Sec. II.2 can be rewritten as

[TABLE]

where $f$ is the free energy density, $\log(Z)/V$. This same factor appears in the final ratio when computing $\braket{\mathcal{O}}$ using Sec. II.3 as well. We see that this ratio depends exponentially on the number of qubits used in the system—note that $f$ is always less than $\log(2)$ here, since the Boltzmann weights will always be less than or equal to one, and $f=\log(2)$ when all Boltzmann weights are one. Therefore, to compute a number as small as $2^{-V}$ using QME will require $O(V)$ qubits in the result register, which will entail $O(2^{V})$ calls to $\mathcal{A}$. Notice this is still quadratically faster than traditional sampling, since naively

[TABLE]

Therefore even in the case of a sign problem or an exponentially small signal, the QME algorithm is superior.

For many state preparation methods, QME depends exponentially on the system volume. This is paid either in the amplitude amplification step, or in the precision of calculating the mean during QPE. This dependence comes from producing the full probability distribution; however, one can instead approximate the distribution. This is relevant to the case of classical sampling algorithms. Using pseudo-random number generators, these algorithms sample from approximate probability distributions. Since classical algorithms can be ported onto a quantum computer [120, 26], these distributions can be realized by an $\mathcal{A}$ circuit. Even still, QME calls $\mathcal{A}$ quadratically fewer times than would be required by sampling, and so calculations of $\langle\mathcal{O}\rangle$ are accelerated—at least asymptotically. Moreover, should classical algorithms improve in the future, QME can take those state preparation methods and use them, along with a further advantage.

In light of this discussion, we can study the resource costs for QME. We do this for the example of computing $\pi$ without a sign problem. In Figs. 17 and 18 we see the total $R_{Z}$ and CNOT gate counts as a function of $N_{r}$, and $N_{g}$, respectively. In both cases we see an exponential dependence. This can be understood as follows: an increase in $N_{r}$ exponentially increases the precision of the result in QPE. Likewise, as we increase $N_{g}$ the grid fineness increases, and the accuracy improves exponentially. The growth in gate count is therefore a reflection of the improvement in the accuracy and precision of the result with growth in both $N_{g}$ and $N_{r}$.

V Conclusion

We have shown how QME can estimate Euclidean LFT observables. The method relies on QPE and provides a quadratic speed-up over traditional methods. The possible improved scaling in the calculation of expectation values is tantalizing. While the gate costs seem to require fault-tolerant quantum computers, further optimizations can reduce those costs. Possible directions include improving QPE [121, 25, 122] over the version used here, or advancements in encoding phases into the result register in other formats, e.g. floating point. Crucially, quantum state preparation of the probability distribution dominates the resource estimate. However when using amplitude amplification, it is still quadratically faster than classically constructing the distribution. Further reductions may occur by fitting the measured phase distribution instead of taking only the most likely value.

Quantum mean estimation also provides an alternative perspective on sampling. When calculating expectation values, the quantum system must be prepared, measured, and re-prepared. Similarly, classical sampling algorithms follow this pattern, albeit obfuscated. Quantum mean estimation instead embeds the preparation and measurement actions into the quantum algorithm itself via quantum superposition. The average is returned with high probability, meaning, the only sampling necessary is used to distinguish the correct bit-string from others. Moreover this speed-up persists even for sign problems and near criticality. This advantage is appealing for LFT practitioners, where these issues are prohibitive. Conversely, the speed-up provided by QME is only polynomial, and hence, the prefactors will dictate its utility. Applying the algorithm on real quantum hardware with simple examples will elucidate its long-term scaling and practicality.

Acknowledgements.

The authors thank Prasanth Shyamsundar for reading drafts of the manuscript, and Michael Wagman for stimulating discussions. This work is supported by the Department of Energy through the Fermilab QuantiSED program in the area of “Intersections of QIS and Theoretical Particle Physics”. Fermilab is operated by Fermi Research Alliance, LLC under contract number DE-AC02-07CH11359 with the United States Department of Energy.