Short-depth circuits for efficient expectation value estimation

TL;DR

This paper introduces a method combining a single quantum phase estimation step with classical post-processing to estimate expectation values more efficiently on near-term quantum devices, reducing circuit depth compared to traditional techniques.

Contribution

It proposes a novel hybrid approach that lowers circuit depth while maintaining accuracy, improving expectation value estimation for quantum algorithms.

Findings

Reduces circuit depth to rac{ ext{target accuracy}}{ ext{power}}rac{ ext{target accuracy}}{ ext{power}} with a single QPE step.

Achieves fewer measurements than standard operator averaging for certain problem instances.

Provides conditions under which the new method outperforms existing strategies.

Abstract

The evaluation of expectation values $Tr\left[\rho O\right]$ for some pure state $\rho$ and Hermitian operator $O$ is of central importance in a variety of quantum algorithms. Near optimal techniques developed in the past require a number of measurements $N$ approaching the Heisenberg limit $N=\mathcal{O}\left(1/\epsilon\right)$ as a function of target accuracy $\epsilon$. The use of Quantum Phase Estimation requires however long circuit depths $C=\mathcal{O}\left(1/\epsilon\right)$ making their implementation difficult on near term noisy devices. The more direct strategy of Operator Averaging is usually preferred as it can be performed using $N=\mathcal{O}\left(1/\epsilon^2\right)$ measurements and no additional gates besides those needed for the state preparation. In this work we use a simple but realistic model to describe the bound state of a neutron and a proton (the deuteron) and show that the latter strategy can require an overly large number of measurement in order to achieve a reasonably small relative target accuracy $\epsilon_r$. We propose to overcome this problem using a single step of QPE and classical post-processing. This approach leads to a circuit depth $C=\mathcal{O}\left(\epsilon^\mu\right)$ (with $\mu\geq0$) and to a number of measurements $N=\mathcal{O}\left(1/\epsilon^{2+\nu}\right)$ for $0<\nu\leq1$. We provide detailed descriptions of two implementations of our strategy for $\nu=1$ and $\nu\approx0.5$ and derive appropriate conditions that a particular problem instance has to satisfy in order for our method to provide an advantage.

Tables

Table 1. Table 1: Rotation (U3) and readout errors for IBM’s 5-qubit machine ’ibmqx4’ on May 8 2019. The error on the CNOT gate on the pairs [ q ​ u ​ b ​ i ​ t ​ 2 , q ​ u ​ b ​ i ​ t ​ 1 ] 𝑞 𝑢 𝑏 𝑖 𝑡 2 𝑞 𝑢 𝑏 𝑖 𝑡 1 [qubit2,qubit1] and [ q ​ u ​ b ​ i ​ t ​ 3 , q ​ u ​ b ​ i ​ t ​ 2 ] 𝑞 𝑢 𝑏 𝑖 𝑡 3 𝑞 𝑢 𝑏 𝑖 𝑡 2 [qubit3,qubit2] is 4.88 % percent 4.88 4.88\% and 6.68 % percent 6.68 6.68\% respectively.

qubit	Rotation err	Readout err
0	0.0019	0.0865
1	0.0024	0.08
2	0.0024	0.0382
3	0.0027	0.3567
4	0.0036	0.2715

Equations

$$⟨O⟩=Tr\left[O\rho \right]=⟨\Psi |O|\Psi ⟩,$$

$$\begin{array}{rl}Tr\left[e^{i\tau O}\rho \right]& =⟨cos\left(\tau O\right)⟩+i⟨sin\left(\tau O\right)⟩\\ & =1+i\tau ⟨O⟩+\cdots \end{array}$$

$$O=O_I+O_T{O}_I\equiv {\alpha }_{0}1Tr[O_T]=0$$

$$O_T=\sum _{k=1}^L{\alpha }_k{U}_k=\sum _{k=1}^L{\alpha }_k\left[e^{i{\theta }_k}{P}_k\right],a_k>0\forall k$$

$$\widehat{O_A}={\alpha }_0+\sum _{k=1}^L{\alpha }_k\widehat{U_k}={\alpha }_0+\sum _{k=1}^L{\alpha }_k{e}^{i{\theta }_k}\widehat{P_k},$$

$$Var[\widehat{O_A}]=\sum _{k=1}^L{\alpha }_k^{2}Var[\widehat{P_k}]=\sum _{k=1}^L{\alpha }_k^2\frac{1-{\widehat{P_k}}^2}{M},$$

$$N_{tot}=\frac{L}{{ϵ}^2}\sum _{k=1}^L{\alpha }_k^2\left[1-{\widehat{P_k}}^2\right]\le \frac{L\Vert \overline{O_T}{\Vert }_2^2}{{ϵ}^2},$$

$$\Vert \overline{O_T}{\Vert }_q\equiv {\left(\sum _{k=1}^L{\alpha }_k^q\right)}^{1/q}\text{for}q\ge 1.$$

$$M_k\propto \frac{{\alpha }_k}{\Vert \overline{O_T}{\Vert }_1},$$

$$N_{tot}^{\prime }=\frac{\Vert \overline{O_T}{\Vert }_1}{{ϵ}^2}\sum _{k=1}^L{\alpha }_k\left[1-{\widehat{P_k}}^2\right]$$

$$N_A(ϵ)=\frac{\Vert \overline{O_T}{\Vert }_1^2}{{ϵ}^2}\lesssim \frac{L\Vert \overline{O_T}{\Vert }_2^2}{{ϵ}^2}$$

$$R_O=\frac{|⟨O⟩|}{\Vert \overline{O_T}{\Vert }_1}\le \frac{\Vert \overline{O}{\Vert }_1}{\Vert \overline{O_T}{\Vert }_1}\equiv R_O^{max}$$

$$N_A({ϵ}_r)=\frac{1}{{ϵ}_r^2{R}_O^2}\ge \frac{1}{{ϵ}_r^2}{\left(\frac{{\lambda }_{max}-{\alpha }_0}{|⟨O⟩|}\right)}^2$$

$$⟨sin\left(\tau O\right)⟩=\tau ⟨O⟩-\frac{{\tau }^3}{6}⟨O^3⟩+\mathcal{O}\left({\tau }^5\Vert \overline{O}{\Vert }_1\right)$$

$$\Qcircuit@C=1em@R=.7em{\ket{0}\quad\quad&\gate{H}\ctrl{1}\gate{S}\gate{H}\qw\\ \ket{\Psi}\quad\quad\qw\gate{e^{i\tau O}}\qw\qw\qw\\ }$$

$$\begin{array}{rl}⟨Z_a⟩(\tau )& =p_0(\tau )-p_1(\tau )\\ & =-⟨\Psi |sin\left(\tau O\right)|\Psi ⟩\end{array}$$

$$O_K(\tau )=-\frac{1}{\tau }\left(⟨Z{⟩}_a(\tau )+\sum _{k=1}^{K-1}{\tau }^{2k+1}\frac{(-1{)}^k{m}_k}{(2k+1)!}\right).$$

$$\begin{array}{rl}{B}_K(\tau )& =O_K(\tau )-⟨O⟩=-\sum _{k=K}^{\infty }{\tau }^{2k}\frac{(-1{)}^k{m}_k}{(2k+1)!},\end{array}$$

$$|B_K(\tau )|\le {\tau }^{2K}\frac{|m_K|}{(2K+1)!}.$$

$${ϵ}_M^2(\tau ,K)=Var\left[\widehat{O_K}(\tau )\right]+B_K(\tau {)}^2,$$

$$N_{tot}=\frac{1}{{\tau }^2}\frac{1-\widehat{Z_a}(\tau {)}^2}{{ϵ}^2-B_K(\tau {)}^2},$$

$${\tau }_{\mathrm{opt}}={\left(\frac{(2K+1)!}{\sqrt{2K+1}}\frac{ϵ}{|m_K|}\right)}^{\frac{1}{2K}}$$

$$N_{sQPE}(K,ϵ)=\frac{|m_K{|}^{1/K}}{{ϵ}^{2+1/K}}f(K)$$

$$f(K)=\frac{2K+1}{2K}{\left(\frac{\sqrt{2K+1}}{(2K+1)!}\right)}^{\frac{1}{K}},$$

$$N_{sQPE}(K,{ϵ}_r)=\frac{f(K)}{{ϵ}_r^{2+1/K}}\frac{|⟨O^{2K+1}⟩{|}^{1/K}}{|⟨O⟩{|}^{2+1/K}},$$

$$\frac{|m_K|}{\Vert \overline{O}{\Vert }_1^{2K+1}}=\frac{|⟨O^{2K+1}⟩|}{\Vert \overline{O}{\Vert }_1^{2K+1}}\le {\Gamma }_K\le 1,$$

$$N_{sQPE}(K,{ϵ}_r)\le \frac{f(K)}{{ϵ}_r^{2+1/K}}\frac{{\Gamma }_K^{1/K}}{R_O^{2+1/K}}{\left(\frac{\Vert \overline{O}{\Vert }_1}{\Vert \overline{O_T}{\Vert }_1}\right)}^{2+1/K},$$

$$R_O\ge \frac{f(K{)}^K}{{ϵ}_r}\frac{|⟨O^{2K+1}⟩|}{\Vert \overline{O_T}{\Vert }_1^{2K+1}},$$

$$R_O\ge \frac{f(K{)}^K}{{ϵ}_r}{\Gamma }_K{\left(1+\frac{|{\alpha }_0|}{\Vert \overline{O_T}{\Vert }_1}\right)}^{2K+1},$$

$$\frac{|{\lambda }_{\Psi }|}{\Vert \overline{O_T}{\Vert }_1}\le \frac{{ϵ}_r^{1/2K}}{\sqrt{f(K)}}.$$

Full text

In this work we use a simple but realistic model to describe the bound state of a neutron and a proton (the deuteron) and show that the latter strategy can require an overly large number of measurement in order to achieve a reasonably small relative target accuracy $\epsilon_{r}$. We propose to overcome this problem using a single step of QPE and classical post-processing. This approach leads to a circuit depth $C=\mathcal{O}\left(\epsilon^{\mu}\right)$ (with $\mu\geq 0$) and to a number of measurements $N=\mathcal{O}\left(1/\epsilon^{2+\nu}\right)$ for $0<\nu\leq 1$. We provide detailed descriptions of two implementations of our strategy for $\nu=1$ and $\nu\approx 0.5$ and derive appropriate conditions that a particular problem instance has to satisfy in order for our method to provide an advantage.

††preprint: INT-PUB-19-016

As we are approaching the era of noisy intermediate scale quantum devices (NISQ Preskill (2018)) the growing interest in practical applications of quantum computing techniques has led to an increased interest in algorithms that are robust against errors and require a limited number of qubits and gates. A central component of many quantum computing algorithms, in particular for applications to quantum simulation Feynman (1982); Lloyd (1996), is the estimation with error $\epsilon$ of the expected value of an hermitian operator $O$ on some quantum state described by the density matrix $\rho$ of a system of $n$ qubits:

[TABLE]

where here and in the following we assume the qubits to be in a pure state described by $\rho=\ket{\Psi}\bra{\Psi}$. A common example is the calculation of the expectation value of the Hamiltonian of a many-body system which can be used, together with the variational principle, to guide the preparation of low energy states in algorithms such as the Variational Quantum Eigensolver (VQE) Peruzzo et al. (2014); McClean et al. (2016); Kandala et al. (2017); Shen et al. (2017). For unitary operators we can evaluate the expectation value in Eq. (1) by simply measuring the state of the device in the appropriate basis. If the change of basis is expensive one can employ the Hadamard Test algorithm from Ortiz et al. (2001) which only requires a single application of $O$ conditional on an ancilla prepared in a superposition state.

In general however we are interested in evaluating $\langle O\rangle$ for hermitian operators $O$ (such as the Hamiltonian of many-body system) and thus alternative strategies have to be devised. Optimal quantum algorithms have been discovered in the past (see eg. Knill et al. (2007)), with efficiencies approaching the Heisenberg limit $N=\mathcal{O}(1/\epsilon)$ for the number of times the experiment needs to be repeated in order to achieve some target additive error $\epsilon$. This remarkable result is achieved by making effective use of both the Quantum Phase Estimation Cleve et al. (1998); Abrams and Lloyd (1999) and Amplitude Amplification Brassard et al. (2000) algorithms, along with the observation that we can estimate Eq. (1) from the linear part of

[TABLE]

up to additive errors that vanish in the limit $\tau\to 0$.

The scheme requires $\mathcal{O}(1/\epsilon)$ applications of a controlled version of the unitary operator $U_{\tau}=e^{i\tau O}$ and it’s inverse (for the QPE part) and $\mathcal{O}(1/\epsilon)$ application of the state preparation unitary $W$ (ie. $W\ket{0}=\ket{\Psi}$). Near term quantum devices will be characterized by a substantial level of noise Preskill (2018) and this will in general prevent us from employing these algorithms even in situations where performing $U_{\tau}$ is simple (eg. for operators $O$ diagonal in the computational basis) as the need to apply $W$ multiple times will lead to an excessively large gate count.

In general, schemes that are robust to noise are therefore preferable for early applications and this has led to a proliferation of hybrid quantum-classical algorithms Bravyi et al. (2016) for a variety of purposes ranging from quantum simulation Li and Benjamin (2017), to approximate optimization Farhi et al. (2014) and quantum compiling Khatri et al. (2019). Our work follows the same philosophy in that we will delegate a substantial computational effort to classical computing resources while at the same time leveraging the available capabilities of current quantum hardware.

The paper is organized as follows. In Sec. I we review the standard method commonly used to evaluate expectation values of Hermitian operators known as Operator Averaging McClean et al. (2014) and in Sec. II we present a detailed discussion of our proposed method based on a single step of phase estimation for performing this task. We then proceed in Sec. III to apply our approach to a simple but challenging nuclear physics problem: the calculation of the deuteron binding energy with a realistic interaction. We have devoted Sec. IV to a more thorough exploration of the effect of noise on both methods and discuss more in depth why an ancilla-based scheme could provide potential benefits on near term noisy devices. As our method is not expected to be competitive in the asymptotically small error limit, in Sec. V we describe how to overcome some of the major implementation challenges in assessing if our scheme could prove advantageous for a particular problem instance. Finally in Sec. VI we summarize the results and possible future directions of the present work.

I Operator Averaging

We discuss here the standard method commonly used in the literature to evaluate expectation values of Hermitian operators Peruzzo et al. (2014); McClean et al. (2016) and start setting up the notations we will use throughout this work. As a first step, and without loss of generality, let us separate out the traceless component $O_{T}$ from the observable

[TABLE]

for some $\alpha_{0}\in\mathbb{R}$. The idea behind the Operator Averaging (OA) approach first proposed in Peruzzo et al. (2014) is to exploit a decomposition of $O_{T}$ of the following form

[TABLE]

where $P_{k}\in\{\mathbb{1},X,Y,Z\}^{\otimes n}$ are tensor products of Pauli matrices and the angles $\theta_{k}$ are introduced in order to keep the coefficient $\alpha_{k}$ positive. The expectation value of each of the terms in the expansion can be estimated with $O(1)$ circuit depth by directly measuring the corresponding Pauli operator on the quantum hardware to get the finite sample estimators $\widehat{P_{k}}$ (here and in the following we will use a wide hat to indicate sample estimators) and combine them to form

[TABLE]

which converges to $\langle O\rangle$ in the infinite measurement limit. In particular, if we perform $M$ measurement for every one of the $L$ terms in Eq. (5) the variance of this estimator is given simply by

[TABLE]

due to the independence of the separate measurements. Note that this implies that in general $Var[\widehat{O_{A}}]\neq Var[O]$ and in fact Eq. (6) may be large even for eigenstates of $O$. Using Eq. (6) we can estimate the total number of measurement $N_{tot}=LM$ required to evaluate $\widehat{O_{A}}$ with precision $\epsilon$ to be

[TABLE]

where we defined

[TABLE]

Different strategies have been proposed in the literature to reduce the scaling of operator averaging shown in Eq. (7). For instance, in situations like quantum chemistry where a large number of coefficients $\alpha_{k}$ have possibly a small magnitude, one can use efficient truncation schemes McClean et al. (2014, 2016) to improve the performance considerably. A complementary approach is to group the $L$ terms into $G$ groups of operators which can be measured together in a single experiment McClean et al. (2016); Kandala et al. (2017); Izmaylov et al. (2019), if the newly introduced correlations are not too large this approach can allow again a great reduction in the number of measurements. In general however we still expect $G$ to scale with system size (for quantum chemistry applications see eg. McClean et al. (2014)). We could also choose the number of measurements to be performed for $k$-th term to depend on the magnitude of the expansion coefficients $\alpha_{k}$ as proposed in Romero et al. (2018):

[TABLE]

this leads in turn to the estimate

[TABLE]

and in this way obtaining a better bound for $N_{tot}$ by a factor $\|\overline{O_{T}}\|_{1}^{2}/(L\|\overline{O_{T}}\|_{2}^{2})\leq 1$. For our nuclear physics problem this strategy produced actually a small increase in the cost since the magnitude of the variance of $k$-th term does not necessarily correlate with the magnitude of coefficient $\alpha_{k}$. By also using the asymptotic improvement from this adaptive variant, we will consider

[TABLE]

as an estimate for the number of repetition needed to obtain a precision $\epsilon$ with operator averaging.

As expected we find that the estimator of Eq. (5) shows in all cases the usual shot noise behavior for small errors $N_{tot}=\mathcal{O}(1/\epsilon^{2})$ but its explicit dependence on the operator norm of $O_{T}$ (which is a lower bound for $\|\overline{O_{T}}\|_{1}$) can be unfavorable when the expectation value $\langle O\rangle$ becomes too small. Introducing the ratio

[TABLE]

we can express the number of shots $N_{tot}$ in terms on the relative error $\epsilon_{r}=\epsilon/\lvert\langle O\rangle\rvert$ as

[TABLE]

which makes explicit the quadratic dependence of the classical effort (the number of repetitions) with the ratio between the largest eigenvalue $\lambda_{max}$ of the target operator $O$ and its expectation value in the state $\ket{\Psi}$.

The scheme we propose in the next section can be advantageous whenever the ratio $R_{O}$ becomes excessively small by providing a scheme with $N_{tot}$ independent of $R_{O}$ in the special case of eigenvalue estimation and possibly well performing in general (see condition Eq. (28)). This could be important in a large system when $\|O_{T}\|$ is extensive (indeed for applications in quantum simulation we expect $\|O_{T}\|=\mathcal{O}\left(poly(n)\right)$ in the number of qubits, see e.g. McClean et al. (2014)) or simply because the expectation value we are after is much smaller than the largest eigenvalue of $O_{T}$ like for the ground state of the deuteron with hard-core potentials Wiringa et al. (1995); Carlson et al. (2015) studied in Sec. III.

II Expectation values from single step phase estimation

The widespread use of the direct algorithm described in the previous section comes from the appealing property of minimizing the required quantum resources by not requiring any additional quantum operation on the qubits. This is especially important for NISQ era devices where coherence time and noise will limit the attainable circuit depth Preskill (2018). As mentioned in the introduction, algorithms based on full fledged quantum phase estimation, like the one described in Knill et al. (2007), will possibly allow to approach Heisenberg limited scaling $\epsilon=\mathcal{O}(1/N)$ as a function of the number of measurements $N$. The price for this is a circuit depth that scales as $C_{D}=\mathcal{O}(1/\epsilon)$ making its implementation challenging on noisy devices with circuit depth $C\lesssim 100$ Preskill (2018). In this section we show how to effectively use a single step of time evolution to obtain a large efficiency gain in measurement number with respect to the operator averaging method, while enjoying short circuit depths $C_{sQPE}=\mathcal{O}(\epsilon^{\alpha})$ with $\alpha\geq 0$.

Similarly to the technique of Knill et al. (2007), we can use the small time expansion of the imaginary part of the expectation value of the time evolution unitary $U_{\tau}=e^{i\tau O}$ on the state $\ket{\Psi}$ (cf. Eq. (2))

[TABLE]

to extract the expectation value $\langle O\rangle$. In particular for any Hermitian observable $O$ we can consider the following standard circuit (cf. Cleve et al. (1998); Ortiz et al. (2001)):

[TABLE]

The imaginary part of $Tr\left[e^{i\tau O}\rho\right]$ can be extracted by measuring the ancilla along the z axis

[TABLE]

where $p_{0}$ ($p_{1}$) are the probabilities of measuring the ancilla in the $\ket{0}$ ($\ket{1}$) state. In practice we need to estimate the expectation value Eq. (16) by performing $N$ independent measurements and computing their average $\widehat{Z_{a}}(\tau)$. For ease of use, from here on we will call the process of extracting $\langle O\rangle$ from a polynomial fit of a set of estimates of Eq. (16) for different times $\tau$ as sQPE for single step quantum phase estimation.

It is important to note now that, due to the presence of a bias for finite values of $\tau$ (coming from the necessity to truncate the series expansion in Eq. (14)), sQPE has a worse asymptotic scaling $N=\mathcal{O}\left(\epsilon^{-(2+\kappa)}\right)$ with $\kappa>0$. Despite this efficiency loss in the asymptotic $\epsilon\to 0$ limit, we will show in the next subsections that the constant factors can be very small when the ratio $R_{O}\ll 1$. provided certain conditions are met (see Eq. (28) and Sec. V.1). For instance in the situation when we prepare an eigenstate of $O$ and we want to estimate it’s eigenvalue $\lambda_{\Psi}$, a bound $R_{O}\leq 1\%$ is sufficient to guarantee an advantage for the simplest possible version of sQPE (see Sec. II.2) up to a respectably small relative error $\epsilon_{r}=\mathcal{O}(10^{-5})$.

We next describe the general high-order sQPE in an idealized setting and in the following two subsection we propose practical implementation for the two lowest order sQPE algorithms where the truncation of Eq. (14) occurs at either the linear (Sec. II.2) or the cubic term (Sec. II.3). For now we will focus our attention on the classical resources (the number of measurements that needs to be performed) and postpone the discussion of the cost associated with the implementation of the (controlled) time evolution unitary $U_{\tau}=e^{i\tau O}$ in Eq. (15) to section Sec. V.2. We can already anticipate that the gate count will be low since the total evolution time $\tau$ has to be kept small enough to minimize the effect of the bias.

II.1 General scaling

Let’s start by considering the idealized case where we know the coefficients $m_{k}=\langle O^{2k+1}\rangle$ for $k={1,\dots,K}$ (in practice we will need to estimate these resulting in a sub-optimal algorithm). We can now use the Taylor expansion in $\tau$ of Eq. (16) to construct the following biased estimator for the expectation value:

[TABLE]

Here the bias comes from neglecting higher order terms with $k\geq K$ in the expansion

[TABLE]

and its magnitude can bound by using Lagrange’s Remainder theorem:

[TABLE]

Due to the presence of this bias for any finite value of $\tau$ we choose to characterize the deviations of our estimator to the exact expectation value by means of the Mean Squared Error (MSE) defined as:

[TABLE]

where, as before, we denote with $\widehat{O_{K}}(\tau)$ a finite population estimator of Eq. (17). The expected total number of measurements required to achieve a final precision (meaning MSE) target $\epsilon$ can then be estimated as was done before in Sec. I. In particular for any $\epsilon>B_{K}(\tau)$ we have the following

[TABLE]

where we used $Var[\widehat{O_{K}}(\tau)]=(1-\widehat{Z_{a}}(\tau)^{2})/(N\tau^{2})$ with $N$ the size of the population used to estimate $\widehat{Z_{a}}(\tau)$. This estimate is minimized with the choice $\tau=\tau_{\rm opt}$ with

[TABLE]

which leads to the following estimate for the total number of measurements needed for sQPE at order $K$:

[TABLE]

where we have defined

[TABLE]

and as promised the shot noise regime $N_{sQPE}=\mathcal{O}(1/\epsilon^{2})$ is recovered only asymptotically for large values of $K$. As anticipated then, for a sufficiently small error we expect the Operator Averaging method of Sec. I to outperform the scheme presented here. The advantage of sQPE is in the possibility of having a much weaker dependence of the expectation value to norm ratio $R_{O}$. To see this let us first rewrite Eq. (23) using the relative error instead

[TABLE]

where we simply used the definition of $m_{k}$. Consider now the special case of eigenvalue estimation where $\ket{\Psi}$ is an eigenvector of $O$ with eigenvalue $\lambda_{\Psi}$ that we want to compute. In this limit the ratio of expectation values on the right-hand side of Eq. (25) is just $1$ and the total number of measurement required by sQPE is completely independent on system considered, and in particular it does not depend on $R_{O}$. In the more general case where $\ket{\Psi}$ is not a single eigenvector we will need an additional ingredient to assess the performance of sQPE: a tight upper bound on the bias Eq. (19) or equivalently of the moment $m_{K}$. For this purpose let us first introduce $\Gamma_{K}$ as an upperbound of the following ratio

[TABLE]

where we remind that $\|\overline{O}\|_{1}=\lvert\alpha_{0}\rvert+\|\overline{O_{T}}\|_{1}$. We can now use $\Gamma_{K}$ to bound the number of measurements as

[TABLE]

and this could be smaller than Eq. (13) when $\Gamma_{K}\ll 1$. To be more quantitative, the sQPE estimator Eq. (17) will become efficient when $N_{sQPE}(K,\epsilon_{r})\leq N_{A}(\epsilon_{r})$ at the desired relative accuracy $\epsilon_{r}$. Using Eq. (13) and Eq. (25), this condition can be written equivalently as

[TABLE]

which can be turned in the following sufficient condition

[TABLE]

where we reintroduced the upperbound $\Gamma_{K}$ described above and wrote explicitly the dependence on the coefficient $\alpha_{0}$ in Eq. (4). We will discuss in detail how to check if the condition Eq. (28) is satisfied in practical application in Sec. V.1 while for now we focus on the simpler situation of eigenvalue estimation mentioned above. In this case the right hand side of Eq. (28) takes a simple form and the full inequality can be written as:

[TABLE]

We now need only a reasonably tight upperbound $\lambda_{u}$ of $\lvert\lambda_{\Psi}\rvert$ to judge when the condition Eq. (30) is satisfied. This requirement is rather loose in practice since the left-hand side is bounded from above by $\lambda_{max}/(\lambda_{max}-\alpha_{0})\approx\mathcal{O}(1)$ while the right hand side becomes large quickly as a function of $K$.

In order to visualize this effect, we plot in Fig. 1 the minimum relative error $\epsilon_{r}$ that can be achieved for a fixed value of the eigenvalue ratio $R_{O}$ using the condition from Eq. (30). As this will depend on the chosen order $K$, in the figure we report the boundaries starting from $K=1$ in the top left corner up to $K=8$ in the bottom right. In particular this means that to achieve lower error rates than the first solid curve we will need to use sQPE with $K>1$ in order to ensure the condition Eq. (30) is satisfied (ie. sQPE could provide an advantage over the Operator Averaging method). Conversely the darkest and innermost region in Fig. 1 is accessible only for $K>8$.

We can deduce a number of interesting conclusions from this figure. For instance, we can see that if we can place an upperbound of $R_{O}$ smaller than $\approx 10^{-2}$, then sQPE with $K=2$ will be more efficient than Operator Averaging down to extremely small relative errors $\epsilon_{r}\approx 10^{-9}$. This is a major improvement from the limit $\epsilon_{r}\approx 10^{-4}$ achievable with the linear method of Sec. II.2. Furthermore, as we can clearly see in the inset, for target relative error at the $1\%$ level, the linear method can be used effectively up to $R_{O}\approx 0.1$ while by increasing the order to cubic (i.e. $K=2$) we can push this up to $R_{O}\approx 0.8$. Given these observations and the increasing difficulty in implementing sQPE efficiently for large $K$, it’s likely that for many practical situations sQPE with $K=1,2$ will be sufficient to achieve a substantial speedup in terms of the number of measurements required to estimate the eigenvalue $\lambda_{\Psi}$. We will provide a similar scaling analysis for the more general problem of expectation value estimation in Sec. V.1 (see Fig. 11).

It is now time to come back to the problem of performing the polynomial fit needed in Eq. (17) in the realistic case where we do not know the high-order coefficients $m_{k}$. For the last coefficient with $k=K$ we can easily use a reasonably tight upperbound $\Gamma_{K}$ to manage the influence of the bias in Eq. (19) as we did before. This allows for a complete algorithm in the simplest case $K=1$ achieving the scaling of Eq. (27). In higher order algorithms we need to estimate the higher order contributions for $k<K$ by collecting data at different values of the time-step $\tau$ and performing a non-linear fit. If no other information is available, we will need at least $K$ values of $\tau$ to properly perform the reconstruction (eg. using the expansion proposed in Knill et al. (2007)).

In the following sections we describe an implementation of the simple $K=1$ (linear) algorithm and a more efficient scheme that adaptively finds the optimal pair $(\tau_{a},\tau_{b})$ of time parameters for the $K=2$ (cubic) case.

II.2 Linear Algorithm

The simplest case is where we neglect the cubic terms in the expansion of the $sin$ so that our estimator Eq. (17) becomes

[TABLE]

In the linear case the optimal time-step Eq. (22) is

[TABLE]

and now, using the estimate of Eq. (25), we can estimate the number of measurements needed as

[TABLE]

where the limit on the right hand side holds when we prepare $\ket{\Psi}$ in any eigenvector of $O$. Note that in this limit, the resources required for a given target relative accuracy $\epsilon_{r}$ are completely independent on the chosen operator $O$ or even eigenvalue $\ket{\Psi}$. For instance at the $1\%$ level we have

[TABLE]

However, in order to get an advantage with sQPE for this situation, the inequality of Eq. (30) has to be satisfied, and that one does depend on the system details.

The problem now is that in general situations we will not be able to calculate $\tau_{opt}$ without at least an approximate estimate for the wanted expectation value $\langle O\rangle$, a good bound $\Gamma_{K}$ on the bias is in fact not sufficient. Even in the simpler case of eigenvalue estimation we can only use an upperbound $\lambda_{u}$ for the absolute value of eigenvalue to compute the approximation

[TABLE]

and this will cause the total number of measurement required to increase:

[TABLE]

This is somewhat acceptable in the eigenvalue estimation case as an upperbound for the eigenvalue is used initially to check if sQPE is at all convenient using condition Eq. (30). The relatively strong dependence of $\widetilde{N}_{T}(1,\epsilon_{r})$ on the value of $\lambda_{u}$ and the difficulty to obtain the optimal time-step in the general case are two of the main problems of the linear algorithm just described and strong motivations for developing the self-consistent algorithm described in the next section.

Before moving to the cubic algorithm, we want to further illustrate the sensitivity of the linear method on the particular choice for the time step $\tau$ by anticipating some results from the deuteron model described in detail in Sec. III. In Fig. 2 we show with a black solid line the analytical estimate (cf. Eq. (7)) for the error scaling of the operator averaging method of Sec. I

[TABLE]

as a function of the total number of measurements $N_{tot}$. As mentioned in Sec. I above the simple adaptive scheme of Eq. (10) produces a slightly worse performance ($\approx 10\%$ larger constant factor, not shown) than the naive operator averaging method.

For sQPE with $K=1$ (or linear method) we show instead the mean squared error

[TABLE]

achievable with different time-steps. The green line corresponds to sQPE with $K=1$ using the optimal choice for the time step Eq. (32) while the red dashed curve shows the detrimental effect of using a worse upperbound $\lambda_{u}=2\lvert\lambda_{\Psi}\rvert$ resulting in $\tau=\tau_{opt}/2$. In both cases sQPE provides an important speedup over operator averaging but this advantage is very fragile. The purple dot-dashed curve shows results obtained using the smallest time-step which could reasonably provide an advantage: the time step for which the variance in Eq. (20) equals an upperbound of the variance of the operator averaging estimator: $1/\tau=\|\overline{O_{T}}\|_{1}$. The performance of operator averaging can be no worse than that, and indeed we see in Fig. 2 that linear sQPE with this time-step requires $\approx 3$ times more measurements than the original scheme. The grey band spans the whole region covered by linear order algorithms with varying time-steps, in particular the upperbound (shown as dotted black line in Fig. 2) corresponds to the worse possible choice for the eigenvalue upperbound $\lambda_{u}=\lambda_{max}$.

The linear algorithm is simple to implement and can provide already important efficiency gains over Operator Averaging whenever we have the ability to make a good choice for the time-step parameter $\tau$. We will now show how we propose to tackle this issue by using the next order sQPE algorithm corresponding to $K=2$.

II.3 Cubic algorithm

In order to use the estimator Eq. (17) for $K=2$ we need to be able to estimate $m_{1}=\lvert\langle O^{3}\rangle\rvert$. In our implementation we achieve this by computing $\langle Z\rangle_{a}(\tau)$ for two different values of the time step and use these to extract both $\langle O\rangle$ and $m_{1}$ using a cubic fit. Given a pair of time steps $(\tau_{a},\tau_{b})$ the outcome of $M$ independent measurements over of the projector $\Pi_{a}=\rvert 0\rangle\langle 0\lvert$ is described by a pair of binomial random variables $X_{a}\sim B(M,P_{a})$ and $X_{b}\sim B(M,P_{b})$ with probabilities given by

[TABLE]

For small values of the time steps we can approximate these distributions with

[TABLE]

Estimators for the two parameters $\mu$ and $\eta$ can be obtained by extremizing the likelihood $L(X_{a},X_{b}|\mu,\eta,\tau_{a},\tau_{b})$ to observe a particular realization $(X_{A},X_{B})$ given the distribution parameters $(\mu,\eta,\tau_{a},\tau_{b})$:

[TABLE]

The resulting maximum likelihood estimators are:

[TABLE]

[TABLE]

where the time-step dependent coefficients are

[TABLE]

We can estimate the variance of these estimators by computing the inverse of the Fisher information matrix

[TABLE]

where the derivatives are taken over $\{i,j\}=\{\mu,\eta\}$ and the parametric dependence of $I(\mu,\eta)$ on the two time-steps $\tau_{a/b}$ has been suppressed for clarity. The results are

[TABLE]

[TABLE]

where the probabilities $P_{a/b}$ will need to be estimated using only a finite sample. In this work we used the Bayesian estimators

[TABLE]

obtained using a slightly informative Beta prior with $\alpha=\beta=1$, but in general any accurate sample estimator $\widehat{P_{a/b}}$ of the true probabilities $P_{a/b}$ will do.

Note that estimators of the variance obtained in this way are in principle accurate only in the limit of large statistics $M\gg 1$, but in practice we found their use to be reasonable for the application studied in this work. Since the adaptive algorithm we describe below relies on their quality, further work on the construction of better (ie. more robust to noise or more rapidly converging) estimators of the fluctuations of $\mu$ and $\eta$ may prove useful.

In addition to these statistical sources of error we also have a bias coming from the approximation $P_{a/b}\to\widetilde{P}_{a/b}$. This can be estimated to be

[TABLE]

A useful upperbound can be obtained by noticing that $B_{\mu}(\tau_{a},\tau_{b})$ is obtained from the remainder of the Taylor expansion in Eq. (39) as

[TABLE]

where the remainder operator $R_{5}$ is defined as

[TABLE]

as follows directly from the integral representation of the remainder of a Taylor series. We can now bound the bias in the cubic algorithm using for instance

[TABLE]

Achieving a tight bound for the bias is generally important as it controls the final efficiency of the method. As for the linear method above, for now we will focus on the special case of eigenvalue estimation while leaving the discussion on how to obtain practical upperbounds in more general situations in Sec. V.1. In the calculations performed in our work we found a weak dependence of the computational effort with the particular choice of estimator for the bias (see Fig. 7) and we will discuss the different options we used in Sec. III.

In the next subsection we present our strategy to determine the time steps for the cubic algorithm.

II.3.1 Optimal determination of the times steps

As was pointed out at the end of Sec. II.2 one of the major drawbacks of the linear algorithm is its sensitivity to the choice of the time-step $\tau$. For the cubic algorithm with solve this issue by using ideas form Optimal Design Chernoff (1972); Loredo (2004). Optimal Design (OD) techniques have been used in a variety of applications to quantum computing ranging from quantum tomography Kosut et al. (2004); Nunn et al. (2010); Huszár and Houlsby (2012), to parameter estimation Ballo´and Hangos (2011); Granade et al. (2012), to quantum-gate synthesis Schulte-Herbrüggen et al. (2005). The general underlying idea in OD for parameter estimation is to try to optimize some, possibly unconstrained, hyper-parameters of an experiment (eg. the pair of time steps ($\tau_{a},\tau_{b})$ to be used in the cubic sQPE) in order to minimize an estimator for the error in the parameter we want to estimate. In many situations this minimization procedure is translated into the maximization of some measure of the “size” of the Fisher information matrix $I(\mu,\eta|\tau_{a},\tau_{b})$ of Eq. (45) (these can be eg. one of its norms or its determinant). This procedure can be seen effectively to be a minimization of the Cramer-Rao bound for an unbiased estimator of the target parameter.

Since in our application the maximum likelihood estimator Eq. (42) has a bias, we will minimize the mean squared error of $\mu_{mle}$ instead:

[TABLE]

One possible adaptive algorithm works by choosing, for any given iteration $i$, a new pair of time steps $(\tau^{i+1}_{a},\tau^{i+1}_{b})$ for the next rounds of $M$ measurements by minimizing Eq. (53) using the estimators $\widehat{\mu}_{mle}$ and $\widehat{\mu}_{mle}$ available at the current iteration $i$. The initial pair can be chosen randomly provided both time steps are small, in the results shown below we sample one of the two from a uniform distribution $U(0,0.1)$111The upperbound on $\tau$ here is somewhat arbitrary since the time step will get readjusted anyway. while the second is chosen to minimize the following upper bound for the variance

[TABLE]

keeping the first fixed. This procedure ensures that the estimator for $\mu$ obtained from the new set of measurements have the smallest MSE possible and is thus rather efficient early on. As we collect more data and the variance of our estimator $\widehat{\mu}_{mle}$ gets reduced we should however reduce the contribution of the bias by reducing in magnitude the new pair of time steps. In order to incorporate this effect we obtain a new set $(\tau^{i+1}_{a},\tau^{i+1}_{b})$ of time steps by minimizing the expected variance of $\widehat{\mu}_{mle}$ after the new block of data is collected

[TABLE]

which has the correct shot-noise scaling coming from $Var[\widehat{\mu}]\approx Var[\mu]/i$ as a function of the number of data blocks $i$ collected so far. In practice we cannot evaluate the variance exactly and for the results presented in this work we found it sufficient to use the approximation

[TABLE]

where we have used

[TABLE]

to estimate the two probabilities $P_{a/b}$ (cf. Eq. (39)). Introducing an upperbound $B^{u}_{\mu}$ for the bias in Eq. (52), the final cost function we use to find the new set of time steps is therefore

[TABLE]

The approximation of the variance Eq. (56) is good in the limit of small time steps $\tau_{a/b}$ but this is not a problem since the presence of the bias term forces the optimal solutions to be numerically small automatically. In order to prevent numerical instabilities early on in the optimization we use a cost function that becomes extremely large whenever $\widetilde{P}_{a/b}\notin[0,1]$.

III The deuteron ground state

The deuteron is the simplest nucleus present in nature. It is a bound state of a neutron and a proton in a state having total isospin $T=0$, spin $S=1$ and angular momentum-parity $J^{\pi}=1^{+}$. It has a small binding energy of approximately $2.2$ MeV. The ground state of the deuteron has a non-zero quadrupole moment, originated by the mixing between s- and d-waves generated by pion-exchanges (see eg. Garcon and W. Van Orden (2001) for a pedagogical introduction). A simple model for the system is to consider the 2-level system built from an s-wave orbital $\ket{\phi_{S}}$ and a d-wave one $\ket{\phi_{D}}$:

[TABLE]

Using the Argonne Av6’ potential Wiringa and Pieper (2002) we obtain (approximately 222The exact matrix elements with the Av6’ potential are real, the integer representation used here is however a very good approximation and simplifies the discussion.) the following Hamiltonian matrix:

[TABLE]

Large cancellations among different contributions produce a ground state energy orders of magnitude smaller than the norm:

[TABLE]

where, in analogy to the previous sections, we defined $H_{T}$ top be the traceless part of $H$ while $R_{O}$ is the ratio defined above in Eq. (12).

This large cancellation is a direct consequence of the hard-core nuclear repulsion that we can see in the central component of the interaction in Fig. 3 which introduces states with vary large energies in the many-body Hilbert space. This is a notorious problem which causes calculations on a finite basis to converge extremely slowly. General strategies to alleviate the issue have been developed in the past (eg. the Similarity Renormalization Group approach Bogner et al. (2007); Hergert and Roth (2007)) but they are usually accompanied by an increase in the degree on nonlocality of the Hamiltonian (see eg. Anderson et al. (2010)) which in general will require a (possibly large) increase in the number of terms needed in expansions of the form Eq. (4).

Note however that, even when the detrimental effects of hard-core interactions are mitigated trough an effective theory like the one mentioned above, the requirement of ensuring basis-size convergence by performing multiple calculations with progressively larger basis sets will still lead to a possibly large mismatch between the ground state energy and the Hamiltonian norm. In that case this is due to fact that, as the basis size increases, the ground state energy will decrease at a much slower rate than the maximum eigenvalue (indeed $E_{gs}$ will reach a plateau for large basis while the highest eigenvalue will grow indefinitely). A general strategy to reduce the importance of this problem (like the sQPE scheme presented in this work) is thus welcome more generally.

Let’s now start to discuss the performance of the Operator Averaging method of Sec. I on our model deuteron problem Eq. (60). Using the estimate from Eq. (13) we find that the number of measurement required for target relative accuracy $\epsilon_{r}$ is given by

[TABLE]

where the last limit holds for $\epsilon_{r}=1\%$. Even tough this estimate might not be very tight since we neglected the variances in Eq. (7) in order to derive Eq. (13), the fact that we are dealing with a simple one-qubit system and that we haven’t considered yet the effect of errors, makes this requirement already alarming. In order to put this number in perspective, the IBM group Kandala et al. (2017) estimated that $N_{A}\approx 10^{6}$ measurements would be sufficient to reach chemical accuracy for a $6$ qubit model of $BeH_{2}$ with hundreds of Pauli terms in the Hamiltonian expansion of Eq. (4).

In Fig. 4 we show results for an ideal implementation (no noise apart from statistical fluctuations) of our one qubit model. As we can see the upperbound of Eq. (61) is only a factor of a few larger and we find that $\approx 9.3\times 10^{6}$ measurements are needed in this ideal noiseless case. Before moving on to discuss the results we have obtained using sQPE, we want to point out that if our goal was only to optimize a variational state $\ket{\Phi(\vec{\theta})}$ using the energy expectation value (ie. we want to run VQE Peruzzo et al. (2014)) then low accuracy results for the energy could be sufficient to get close to the optimum $\vec{\theta}_{min}$. To illustrate this we show in the inset the deuteron energy expectation value obtained using $N_{tot}=10^{3}$ as a function on the iteration in the minimization procedure (for these results we used a simple Nelder-Mead optimizer). Even though the error in the energy is $\gtrsim 200\%$ the angle $\theta$ converges towards $\theta_{min}$ to within a few percent error in only a small number of iterations (note that for this simple model a single angle is sufficient to prepare the ground state). This striking difference is probably peculiar to simple models like our one-qubit deuteron, since the lack of excited states with low energy in the spectrum of the Hamiltonian gives rise to large gradients in the the variational energy $E(\theta)=\langle\Phi(\theta)\lvert H\rvert\Phi(\theta)\rangle$ and therefore to a relatively easy optimization. Where sQPE could be most useful in this case is for the final estimation of the energy, but in general for more complex systems low order sQPE could be advantageous also in the last stages of optimization where large statistical fluctuations could prevent to reach the minimum.

Let’s now turn to discuss the sQPE method starting with the implementation of the controlled time evolution from Eq. (15). In our simple two qubit situation the circuit for the controlled unitary can be constructed using only 2 CNOT gates (cf. Barenco et al. (1995)):

[TABLE]

where $\Phi(\alpha)=diag\left(1,e^{i\alpha}\right)$ is a phase gate and the $R_{A},R_{B},R_{C}$ blocks are formed by appropriate single qubit rotations (see Appendix A for more details).

In Fig. 5 we present the results obtained with both sQPE and operator averaging: the solid curves correspond to empirical results while the dotted lines correspond to the analytical estimates discussed in Sec. II.2. As for Fig. 2, the red square marks the location of the upperbound $N_{A}$ from Eq. (61). For the linear method, the optimal choice Eq. (32) is shown in blue, while in green we present the results obtained using a more conservative value $\tau=\tau_{opt}/2$ (cf. Eq. (35) and the discussion following it). In both cases we see that, in agreement with the results presented in Fig. 2, the linear algorithm requires about an order of magnitude less measurements than operator averaging. In addition to this we see that the estimated number of measurements $N_{sQPE}$ are in very good agreement with the empirical results: the maroon diamond indicates the upperbound $N_{sQPE}(1)$ for linear sQPE Eq. (34) while the cyan circle corresponds to the worse bound obtained trough Eq. (36).

As we discussed in Sec. II.2 the speedup offered by the linear method is very sensitive to the particular choice of time step used and by employing the adaptive strategy described in Sec. II.3 to find the time steps $(\tau_{a},\tau_{b})$ the cubic algorithm partially overcomes this problem. In order to implement the algorithm we use the same circuit described for the linear method: we just run it twice for the two time steps separately. The results obtained from $6$ different runs are presented in Fig. 5 as purple lines, in all cases we update the time step pair every block of $N_{b}=40$ measurements. In addition, the blue triangle is twice the expected number of measurements $N_{sQPE}(2)$ obtained from the general expression Eq. (23)

[TABLE]

where the factor of $2$ is introduced to account for the fact that in our adaptive scheme we are actually estimating two expectation values: $\langle O\rangle$ and $m_{1}$. The optimal pair of time steps $(\tau_{a},\tau_{b})$ obtained during the execution of the algorithm fluctuate around the value $(0.15,0.3)$ which is not very far from the optimal time step $\tau_{opt}\approx 0.4$ found from Eq. (22). The spread of results at large measurement count is possibly a signature that the optimization of Eq. (58) gets stuck in local minima, we plan to investigate this further in future work.

As explained in Sec. II.3, the cubic algorithm needs a good approximation of the bias term $B_{\mu}(\tau_{a},\tau_{b})$ in Eq. (52) as this enters directly the cost function Eq. (58) used to determine the optimal time steps. The results presented above where obtained using the following ansatz

[TABLE]

where $\widehat{\mu}_{mle}$ and $\widehat{\eta}_{mle}$ are the current best estimators for the distribution parameters $(\mu,\eta)$. This form reduces to the correct one in the eigenvalue estimation limit relevant here where $\ket{\Psi}$ is an eigenstate of $O$. We want now to present results showing the weak sensitivity of the cubic algorithm to the specific choice of the estimator for the bias, in particular we will use two additional estimators: the exact one from Eq. (49)

[TABLE]

which in practical situations we won’t have access to, and a different variant of the estimator $B_{A1}$ above defined as

[TABLE]

This estimator is a tighter bound that can be obtained from Eq. (51) using the additional condition

[TABLE]

In Fig. 6 we show how these estimators evolve as the algorithm proceeds in three different situations. The top panel shows the ideal situation where the optimal pair $(\tau_{a},\tau_{b})$ for the next step is obtained using the exact scaled mean squared error (cf. Eq. (55))

[TABLE]

together with the exact bias $B_{E}$. Even though in practical situation we won’t be able to run the cubic algorithm this way, these results provide a ceiling for the performance of approximate algorithms while at the same time show clearly the source of the advantage that is achieved with sQPE: initially the time steps are raised to relatively large values in order to reduce the shot noise limited variance term in the equation above at the expense of a larger bias term. This allows to quickly reduce the error in the expectation value early on when the dominant contribution are statistical fluctuations, as the accuracy increases the importance of the bias term grows and the adaptive algorithm starts to reduce the magnitude of the time steps in order to keep $B_{\mu}$ under control. Furthermore we see that the three bias estimators follow each other rather closely.

The other two panels instead correspond to results obtained using the approximate cost function $\Delta_{i}$ from Eq. (58) with either the $B_{A1}$ ansatz used also in Fig. 5 (central panel) or the $B_{A2}$ ansatz described above (bottom panel). In both situations we recover the same qualitative behavior seen in the ideal case: the bias gets initially increased and then reduced gradually as the accuracy improves. The main difference with the results of the top panel is the lower efficiency obtained in the first stage of this procedure where the actual exact bias remains much smaller than it could have been for the first few hundred measurements (note that as above we update the time-steps every $N_{b}=40$ measurements) but the discrepancy quickly vanishes later when the bias becomes the limiting factor.

The see the impact of these approximations on the final convergence of the expectation value we show in Fig. 7 the results obtained by using the exact pair of time steps $\tau_{a/b}$ obtained as before using Eq. (55) but different approximations to the final mean squared error. The reason this is important is that we need to estimate the bias in order to provide an estimate for the final error of our estimated expectation value. The results in the main panel show the apparent reduced efficiency that is a result of using a bigger bias than the exact one (in this case we used $B_{A1}$ and $B_{A2}$ as above). The cyan up triangle is the same as the blue triangle in Fig. 5 and indicates an optimistic expectation on the efficiency of the cubic algorithm. Somewhat not surprisingly we achieve that estimate only in the ideal exact case shown as black lines in the main panel of Fig. 7 while in the worst case (corresponding to the red curves) we need to perform as much as $\approx 50\%$ more measurements. Interestingly when statistical fluctuations in the estimation of the variance are included these differences mostly disappear and there is no clear preference for different choices of the bias, we can see this from the results shown in the inset of Fig. 7 where the time steps were estimated using the approximate cost function $\Delta^{i}$ from Eq. (58).

These results are encouraging as they show that, even if a tight upperbound for the bias is helpful for the algorithm, an approximate expression can work very well at least for the special case of eigenvalue estimation. Evaluating upperbounds becomes more important in the general case and we leave the discussion for Sec. V.1 while we now turn to the problem of accounting for the presence of noise in the quantum device.

IV Effect of noise

As we have mentioned in the introduction, noise will be unavoidable for near term quantum devices and it is therefore critical for algorithms to provide robustness against noise if we want to deploy them on a non fault-tolerant quantum computer. Since the methodology we propose goes against this trend in that we are trading classical resources (the number of experimental trials) with quantum ones (one more qubit for the ancilla and more gates), we need to provide supporting evidence that our method shows advantages even in the presence of noise and is thus practical.

The importance of this assessment is critical as there are known cases where the advantage of a quantum algorithm can be drastically reduced by the presence of even small noise sources (see eg. Shenvi et al. (2003); Kołodyński and Demkowicz-Dobrzański (2010); Demkowicz-Dobrzański et al. (2012); Fröwis et al. (2018)).

We start by showing how, in situations where condition Eq. (28) is valid by a large margin (as in the deuteron model discussed here), a substantial increase in classical resources is required to minimize the effect of measurement noise for the Operator Averaging method of Sec. I while with sQPE this problem can be mitigated. In the last part we provide a more general argument in support of our measurement strategy in situations where a clean ancilla is available in the spirit of the D1QC model Knill and Laflamme (1998).

IV.1 Measurement noise

Assignment errors in the measurement device used for qubit read out are an important source of bias that needs to be accounted for properly in order to obtain meaningful results. In order to illustrate the problem we show in Fig. 8 results obtained by executing the deuteron problem described in the previous section on a emulated version of the IBM 5-qubit machine ’ibmqx4’ using the Qiskit software package Aleksandrowicz et al. (2019) (see Tab. 1 for details).

Despite the fact that we are not correcting for any source of errors in these results, the linear sQPE (with optimal time $\tau$) presented on the right seems to provide considerably higher quality results despite the much increased circuit depth. In the following we provide an argument to explain the observed results.

Here we will use an extremely simplified model for these errors that nevertheless captures their essential features, this is achieved by replacing the projectors $\Pi_{0}$, $\Pi_{1}$ on the states $\ket{0}$, $\ket{1}$ of a qubit with the following ones

[TABLE]

where $0<p<1$. This model can be justified in the limit where assignment errors are both qubit independent and symmetric with respect to the interchange $\ket{0}\leftrightarrow\ket{1}$ and is sufficient for our purpose (see eg. Supplemental Material of  Kandala et al. (2017) for details on a more accurate model). Using Eq. (69) we find that the noisy expectation value of some one-qubit Pauli operator $\langle\widetilde{P_{\sigma}}\rangle$ is related to the noise-free value by the relation:

[TABLE]

which can be easily inverted to estimate $\langle P_{\sigma}\rangle$ from $\langle\widetilde{P_{\sigma}}\rangle$. Despite it’s simplicity this model is sufficient to completely account for the error afflicting the OA results of Fig. 8 as we can see from the error mitigated results presented in Fig. 9. Note that the linear sQPE energies in the right panel are still biased due mostly to the noise introduced by using the CNOT gates, the mitigation of which is beyond the scope of our discussion here (note however that mitigation techniques Temme et al. (2017); Endo et al. (2018) will be required also for OA for larger target systems).

Generalizations of Eq. (70) to expectation values of multi-qubit Pauli operators can also be obtained Kandala et al. (2017), but we will limit our discussion to the one-qubit case relevant to our deuteron calculations and to sQPE more generally. In fact, an important feature of sQPE is that the output of the algorithm is obtained through measurements on a single qubit thus avoiding the problem of exponentially reduction of signal to noise ratio as a function of the number of qubits involved in the measurement of individual terms in the expansion Eq. (4) (see eg. Kandala et al. (2017)).

For a generic one-qubit observable

[TABLE]

we can now use Eq. (70) to estimate the noise free expectation value

[TABLE]

where $\widehat{p}$ is a finite sample estimator estimator, with variance $\delta p$, of the error probability $p$. The variance of this estimator can be approximated as

[TABLE]

with

[TABLE]

where we used a linear expansion to propagate the error (ie. we used $Var[f(x)]\approx f^{\prime}(x)^{2}Var[x]$). The second error term in Eq. (73) comes from the uncertainty in the determination of the error parameter $p$ and provides a noise floor that we need to minimize in order to achieve good accuracies. Assuming the sample estimator $\widehat{p}$ was obtained from $N_{C}$ calibration measurements, we can bound the contribution of this background as

[TABLE]

where in the first step we used the decomposition Eq. (4).

The case of the sQPE algorithm is simpler because we have always to deal with a single qubit to be measured. Using the same correction scheme employed above, assuming again that the higher order coefficients $m_{k}$ are known, we find

[TABLE]

with variance

[TABLE]

and

[TABLE]

It is then clear that the sQPE algorithm will reduce the importance of measurement noise in the same situations where it provides an advantage in the noise free case, namely whenever $\|\overline{O_{T}}\|\gg\frac{1}{\tau^{2}}$ with $\tau$ the estimated optimal time step for a particular problem.

In order to assess the practical impact of this error term for the deuteron calculation of Sec. III we have numerically minimized the total number of measurements $N_{tot}=2N+N_{C}$ needed to achieve $\epsilon_{r}=1\%$ ($N$ for each of the Pauli terms and $N_{C}$ to estimate $\widehat{p}$) as a function of the error parameter $p$ using directly Eq. (73). In Fig. 10 we show the results of this study: the full black line is the minimal value of $N_{tot}$ needed to reach a target relative error $\epsilon_{r}=1\%$, while the green dashed line corresponds to the situation where we have performed a calibration using $N_{C}^{0}=10^{7}$ initial measurements to estimate the error rate $p$ of the machine.

The inset shows the ratio of tuning measurements needed to achieve the target accuracy as a function of the error rate $p$ in both situations. We see that small values of the error rate $p\lesssim 5\%$ can be dealt with with relative ease using data obtained from previous calibrations, while more substantial efforts are needed for more noisy qubits. Using $qubit3$ with its large error $p\gtrsim 35\%$ would require two orders of magnitude more measurements than in the noiseless case (shown as the red square in Fig. 10) out of which more than $\approx 80\%$ would be needed for characterizing the noise. Fortunately noise levels so high are not common on modern machine, with typical values of a few percent in the superconducting circuit case, but the fact that this very simple source of error is capable of completely swarming the results of a simple single qubit calculation provides another motivation to explore the use of the sQPE scheme on near term architectures.

IV.2 General advantage of ancilla-based schemes

The purpose of this section is to show that the ancilla-based construction of sQPE (cf. Eq. (15)) can be advantageous in general when we want to estimate the value of an expectation value $\langle O\rangle$ in presence of depolarizing noise (see eg. Nielsen and Chuang (2010)) in the quantum device. For a realistic advantage to be found we will assume we have either a clean ancilla (meaning completely error free) or at least a qubit subject to a well characterized noise channel and with high measurement fidelity.

Let’s start by considering a slight generalization of the sQPE circuit Eq. (15) where we leave unspecified the state of the ancilla before we apply the controlled unitary

[TABLE]

and in the second circuit we made explicit the presence of a new rotation matrix $W$ which prepares $\ket{\phi}$ starting from $\ket{0}$. The action of the time evolution unitary $U_{\tau}$ on the target state $\ket{\Psi}$ can be conveniently expressed as

[TABLE]

with $\kappa,\nu\in\mathbb{C}$ with

[TABLE]

In other words, the total Hilbert space explored with circuits of type Eq. (79) is only $4$-dimensional: we have $\mathbb{C}^{2}$ for the ancilla and the linear span of $\ket{\Psi}$ and $\ket{\Psi^{\perp}}$ for the target system. At this point, a measurement of the $y$-polarization of the ancilla after the circuit in Eq. (79) will reveal the wanted quantity (cf. Sec. II)

[TABLE]

from which we can extract the expectation value as discussed above.

If we trace out the system qubits, the circuit in Eq. (79) can be represented as a quantum channel $\Lambda_{\tau}$ acting on the ancilla:

[TABLE]

where the output state of the ancilla is indicated here as a density matrix $\rho_{f}$. We want to show now how by performing quantum process tomography Chuang and Nielsen (1997); Altepeter et al. (2003); Mohseni and Lidar (2006); Mohseni et al. (2008) on the ancilla we can extract $\kappa_{I}$.

We can represent the quantum channel $\Lambda_{\tau}$ using the following Kraus decomposition (see eg. Nielsen and Chuang (2010))

[TABLE]

with

[TABLE]

but this choice is not unique. A better parametrization of the channels that overcomes this difficulty is to use the Pauli Transfer Matrix Chow et al. (2012) defined as

[TABLE]

where the $P_{i}$’s are the Pauli operators $\{\mathbb{1},X,Y,Z\}$. For our quantum channel this matrix takes the form

[TABLE]

where $\kappa=\kappa_{R}+i\kappa_{I}$. This form makes it apparent that the channel $\Lambda_{\tau}$ is a composition of a dephasing (or phase damping Nielsen and Chuang (2010)) channel

[TABLE]

with error probability $p_{z}=\lvert\nu\rvert^{2}=1-\lvert\kappa\rvert^{2}$ and a rotation around the z-axis with angle $\theta=tan^{-1}\left(\kappa_{I}/\kappa_{R}\right)$ described by

[TABLE]

As shown by Wiebe et al. in Wiebe et al. (2014) a Bayesian reconstruction strategy can be effectively employed to learn $\kappa_{R}$ and $\kappa_{I}$ from measurement of the device even in the presence of substantial depolarizing noise described by the channel

[TABLE]

where $d=2^{n}$ and $n$ is the number of qubits used to encode $\ket{\Psi}$. According to the results in Wiebe et al. (2014), large values $p_{D}=50\%$ could be handled with relative ease.

Furthermore, the reconstruction is simplified in our case since the structure of our channel is known beforehand, and one can tailor strategies aimed at estimating matrices of the form Eq. (87). The extent to which these could be leveraged to minimize the negative effect of more realistic noise channels in the system qubits is left for future explorations. Before concluding, we want to point out that the possibility of removing read-out lines from all but the ancilla qubit could also help more generally in reducing the overall noise in the device.

V Implementation Challanges

During our exposition of our methodology in Sec. II we have only briefely touched upon the practical cost of implementing the core parts of the algorithm on near term quantum devices. This section is dedicated to address those issues. In particular we first present a discussion on how to estimate the the potential gain of using sQPE using condition Eq. (28) in practical situations where only partial information on the high order coefficients $m_{k}$ is available. We then provide a description of the resources needed to implement the time evolution needed for sQPE (cf. Eq. (15)) using different strategies.

V.1 Practical bound estimation

Due to the presence of the the expectation value $\langle O^{2K+1}\rangle$, it is difficult in most situations to asses directly if the condition in Eq. (28) holds. In order to obtain a more manageable condition we can rewrite Eq. (28) as

[TABLE]

Due to fact that $h(x)=x^{K}$ with $K\geq 1$ has a bounded first derivative on a finite interval $\Omega$ we have

[TABLE]

which in turn implies

[TABLE]

with $\lambda_{max}$ the largest singular value of the operator $O$. Using the bound

[TABLE]

which can be obtained by using Jensen’s inequality, we arrive at the following, looser, condition

[TABLE]

Reasonably tight ubberbounds on $Var[O]$ and $\langle O^{2K}\rangle$ for small $K=\mathcal{O}(1)$ can be obtained in many situations of interest. An important example are many-body calculations of ground-state properties where a variational calculation with a classically simulatable trial states can provide such bounds with reasonable efficiency (eg. one could use Quantum Monte Carlo methods Foulkes et al. (2001); Carlson et al. (2015)). Another situation is when we have some control on the fidelity of the prepared state. For instance consider the case where we are preparing an initial state $\ket{\Psi}$ with large overlap with some eigenvector $\ket{\phi}$ of $O$

[TABLE]

with $\langle\phi|\phi^{\perp}\rangle=0$ and $\alpha^{2}+\beta^{2}=1$. We can obtain a bound on $\langle O^{2K}\rangle$ if we have a upperbound for the state fidelity

[TABLE]

by using

[TABLE]

When only a bound on the variance is available instead we can use the inequality

[TABLE]

to obtain the even looser condition

[TABLE]

Unfortunately this condition can bee too loose to be of practical value as we can see by looking at the limit $Var[O]\to 0$: the approximate expression Eq. (95) recovers the correct limit Eq (30) as does the exact expression Eq. (91), while the right hand side of Eq. (100) is always larger than $f(k)^{K}/\epsilon_{r}$ and this can produce an overly pessimistic assessment of the efficiency gain achievable with the sQPE scheme of Sec. II.

To get more insight on this problem we plot in Fig. 11 the regions in a two-dimensional $\left(Var[O],\epsilon_{r}\right)$ space where the inequalities above predict an advantage of our proposed scheme for the estimation of the deuteron’s ground state energy. The solid black lines define the maximum relative error $\epsilon_{r}$ achievable for a trial state $\ket{\Psi}$ with a given variance obtained using the condition of Eq. (95). The simple linear algorithm of Sec. II.2 is predicted to be efficient for situations that stay in the top left corner of parameter space bounded by the first black line. For lower target $\epsilon_{r}$ or larger variance $Var[O]$ we then progressively need to increase the order $K$ of the algorithm to ensure Eq. (95) is satisfied. From Fig. 11 we also see that for a $1\%$ target error the simple linear scheme with $K=1$ is predicted to be more efficient up to $Var[O]\approx E^{2}_{GS}$ before the growth of the bias term forces us to increase the order in $K$. This is encouraging since this condition is not necessarily tight in the sense that, due to the use of the upperbound Eq. (94), the inequality of Eq. (95) is a sufficient but not necessary condition for Eq. (91) to hold true. The extent to witch one can still find an efficiency gain by using the linear algorithm past this condition will most likely depend on the particular problem instance and is therefore difficult to predict without some prior knowledge on $Cov[O^{2K},O]$ that would allow to use Eq. (91).

The red contours in Fig. 11 are obtained instead by using the looser condition Eq. (100) and as we can see they are overly pessimistic: a relatively high order algorithm with $K>4$ is judged to be needed even in the limiting case when the trial state has zero variance. This observation reinforces the importance of being able to use a tighter condition like Eq. (91) in order to assess meaningfully the possibility of a gain in using the strategy proposed in this work.

V.2 Time evolution

We finally turn our attention to the problem of estimating the circuit depth required for the implementation of the time-evolution unitary $U_{\tau}=e^{i\tau O}$ needed for the sQPE method. In particular we consider the more realistic situation where we only have an approximation $\widetilde{U_{\tau}}$ of $U_{\tau}$ with error bounded by $\delta_{\tau}$ is available:

[TABLE]

In this situation the induced error on the sQPE estimator Eq. (17) is then $\epsilon_{\tau}=\delta_{\tau}/\tau$. A simple way to control the total error $\epsilon$ of the calculation is to require that $\epsilon_{\tau}<\epsilon/2$ and similarly for the MSE $\epsilon_{M}$ in Eq. (20). The latter modification will increase the bound reported in Eq. (23) only by a factor $2^{1+1/K}$. This requirement can be relaxed if one employs optimal algorithms like Berry et al. (2015); Hao Low and Chuang (2016); Low and Chuang (2017) that allow to implement $\widetilde{U_{\tau}}$ for $\epsilon_{\tau}\ll\epsilon$ with only a small increase in gate count.

Due to the relatively short propagation time required by the sQPE algorithm, we can also use simpler strategies based on the Trotter–Suzuki decomposition Suzuki (1991) while still maintaining a short circuit depth. To see how this works let us start by considering a simple first order scheme obtained by dividing the propagation time $\tau$ into $r$ segments:

[TABLE]

where we used the decomposition Eq. (5) for the operator $O$. As shown in Childs et al. (2018) the error for this approximation can be bounded by

[TABLE]

with the the norm $\|\overline{O}\|_{1}=\lvert\alpha_{0}\rvert+\|\overline{O_{T}}\|_{1}$ defined as in Eq. (8). Note that this is slightly tighter than the result Proposition F.3 obtained in Childs et al. (2018) and can be found following the same proof. We can now provide the following analytic bound for the number of intervals $r$ required to guarantee that $\epsilon_{\tau}<\epsilon/2$

[TABLE]

where we used the same derivation as in Childs et al. (2018). For convenience we rewrite the optimal time step Eq. (22) as

[TABLE]

where we defined

[TABLE]

Note the additional factor of $2$ in the denominator coming from the choice $\epsilon_{M}=\epsilon/2$. We can now express the bound on $r$ for the simple linear product formula as

[TABLE]

where we have defined

[TABLE]

and for reasonably small errors the bound Eq. (107) is maximized with the rightmost expression giving the algorithm an overall depth scaling at best as $r_{1}=\mathcal{O}\left(1/\sqrt{\epsilon}\right)$. It is instructive to express these bounds in terms of the relative error $\epsilon_{r}$ and the expectation value ratio $R_{O}$ from Eq. (12), for instance in the case of eigenvalue estimation we find

[TABLE]

which leads to a generic scaling given by

[TABLE]

which in the linear case of Sec. II.2 simplifies to

[TABLE]

This provides only a minor advantage over the $\mathcal{O}\left(1/\epsilon\right)$ scaling associated with full fledged QPE algorithms Knill et al. (2007) which can be easily spoiled with a sufficiently small $R_{O}$ ratio. It is therefore important to use higher order expansions that are able to achieve a more favourable scaling and in the following we will consider higher order product formulas as an example. If we denote the ($2j$)-th order Trotter-Suzuki formula with $r$ intervals Suzuki (1991); Childs et al. (2018) as $S_{2j}(\tau,r)$ we can generalize the error bound Eq. (103) obtained above to

[TABLE]

and bound the number of intervals as

[TABLE]

with $\rho_{j}\equiv 2\rho_{1}5^{j-1}$. Again this is slightly tighter than the result obtained by Childs et al. in the Supplemental Material of Childs et al. (2018). As for the linear decomposition the right term dominates for reasonably small errors $\epsilon$ and we find the overall scaling

[TABLE]

As above we can express this in terms of relative quantities in a compact way for the special case of eigenvalue estimation

[TABLE]

and due to the fast growth of the first term in the above expression we might want to keep the order $j$ as small as possible. For instance using $j=K$ will already guarantee a gate count independent on the target precision $\epsilon$ and scaling as $\mathcal{O}\left(1/R_{O}^{3/2}\right)$ in terms of the eigenvalue ratio. Notably by simply choosing $j=K+1$, for the price of a fixed increase in cost of less than a factor of $5$ we can achieve a circuit depth that decreases as a function of the target relative error.

Note that these estimates are based on the possibly very pessimistic bounds Eq. (103) and Eq. (112) which means that these circuit depths could possibly be greatly reduced in practice (see eg. Childs et al. (2018)). Before finishing this section we want to point out that even though the estimates provided above are for the implementation of $U_{\tau}$ a complete implementation of it’s controlled version needed in Eq. (15) can be obtained with an overall linear increase in depth as a function of the number of qubits in the system register used to represent $\ket{\Psi}$. Tighter bounds will require further knowledge of the particular operator $O$ whose time evolution we want to simulate, but the results presented here give us good reasons to believe practical implementations could be achievable for interesting systems on near term devices.

VI Summary and Conclusions

In this work we reviewed the standard methodology of Operator Averaging Peruzzo et al. (2014); McClean et al. (2014, 2016) to evaluate expectation values $\langle O\rangle$ of general Hermitian operators $O$ efficiently on quantum computers by minimizing the number of quantum operations needed while maintaining a shot-noise limited number of measurement $N_{tot}=\mathcal{O}\left(1/\epsilon^{2}\right)$ to attain precision $\epsilon$ in the expectation value estimate. This provides a great advantage on current generation noisy devices where the asymptotically optimal behaviour $N_{tot}=\mathcal{O}\left(1/\epsilon\right)$ of methods that employ Quantum Phase Estimation Abrams and Lloyd (1999); Knill et al. (2007) cannot be attained in practice due to the large circuit depths $C_{QPE}=\mathcal{O}\left(1/\epsilon\right)$ involved. As we explain in Sec. I however, the Operator Average strategy has a major drawback in that in terms of relative error $\epsilon_{r}$ the total measurement count grows as $N_{tot}=\mathcal{O}\left(1/(\epsilon_{r}R_{O})^{2}\right)$ where $R_{O}$ defined in Eq. (12) is approximately the ratio between the wanted expectation value and the largest eigenvalue of $O$.

In this work we propose to use a single step of phase estimation as in the well known Hadamard Test to learn the expectation value $\langle O\rangle$ by looking at the short time behaviour of $\langle sin(\tau O)\rangle$ instead. This strategy was already discussed in the context of full QPE calculations in Knill et al. (2007) and it remains the method of choice for fully error-corrected devices capable of executing accurately very long gate sequences. Our contribution is in showing how, by using circuits implementing only a single Hadamard Test with appropriately chosen time-steps, one can greatly reduce the classical cost (the number of experimental measurements $N_{tot}$) using much shorter circuits than those needed for QPE. For instance in the important case of eigenvalue estimation we can achieve $N_{tot}=\mathcal{O}(1/\epsilon_{r}^{2+1/K})$ for $K=\mathcal{O}(1)$ independent on $R_{O}$ while keeping the gate count bounded by $C_{sQPE}\mathcal{O}\left(5^{K}/R_{O}^{1+1/2K}\right)$ using a very simple general purpose strategy employing the high-order Trotter-Suzuki decomposition Suzuki (1991); Childs et al. (2018). As we argue in Sec. V.2 the latter requirement can possibly be further reduced by using more advanced simulation strategies Berry et al. (2015); Hao Low and Chuang (2016); Low and Chuang (2017) and we plan to further this possibility along the lines of the study presented in Childs et al. (2018) in a future work.

We presented a complete analysis of the first two lowest order sQPE algorithms with $K=1$ and $K=2$ in Sec. II together with a self consistent procedure aimed at finding the optimal time-steps to be used in the calculation. As our approach could be extremely helpful in some situations but it is not efficient in an asymptotic scaling sense in general, we provide both strict and easy to estimate conditions to help determine if the use of sQPE con provide a speedup for a particular problem instance (see Sec. II and Sec. V.1). As these conditions require the availability of bounds on the expectation value to be computed and some control over the operator spectrum (like bounds on the $n$-th cumulant $\langle O^{n}\rangle$) further work on classically efficient strategies to estimate them (using for instance ideas from Baumgratz and Plenio (2012); Harrow and Montanaro (2017)) could have a possible big impact on the practicality of our approach. As discussed in Sec. V.1 classical Quantum Monte Carlo simulations could be employed efficiently in the meantime. Finally in Sec. IV we have shown some evidence on the robustness of our proposed methodology to readout noise on the quantum device and provided arguments to justify the expectation that ancilla based algorithms like sQPE provide in general a much more robust layout to deploy and execute non trivial quantum algorithms on NISQ devices. It will be very interesting in the future to see the impact of adaptive machine-learning techniques as those presented in Wiebe et al. (2014) on the practical feasibility of scaling up quantum computations in the near term.

Acknowledgements.

We wish to thank N. Klco and R. Schiavilla for helpful discussions regarding the subject of this work. We also thank A. Matsuura of Intel Labs for useful conversations and for encouragements. The work of A.R. was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research (ASCR) quantum algorithm teams and testbed programs, under field work proposal number ERKJ333 and by U.S. Department of Energy grant No. DE-FG02-00ER41132. We also grateful to the T2 group at Los Alamos National Laboratory for the hospitality while working on this project thanks in part to funds from the U.S. Department of Energy, Office of Science, HEP Contract No.DE-KA2401032. A.B. acknowledges support from U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Award No. DE-SC0010300 and DE-SC0019647.

Appendix A Implementation of controlled time-evolution

We report in this section the implementation of the two-qubit controlled time evolution appearing in Eq. (15) needed for the sQPE algorithm and schematically presented in Eq. (79) in the main text. Throughout this section we assume that the initial state $\ket{\Psi}$ has been prepared with a rotation $R_{y}(\theta)$ with $\theta$ the angle of interest. For a system of two qubits a generic controlled unitary operation associated with a $2\times 2$ unitary matrix $U$ can be represented in the computational basis $\rvert 00\rangle$, $\rvert 01\rangle$,$\rvert 10\rangle$, $\rvert 11\rangle$ as

[TABLE]

with $1_{2\times 2}$ and $0_{2\times 2}$ indicating the two-by-two identity matrix null matrix respectively.

Let’s now recall the general decomposition of a $U(2)$ unitary

[TABLE]

for appropriately chosen angles. The rotation matrices here are defined as

[TABLE]

Using this decomposition, together with the definition Eq. (116), we can implement the controlled time-evolution operator with the following circuit (see eg. Barenco et al. (1995))

[TABLE]

where we have defined the phase gate as

[TABLE]

Using the following decomposition of the Hamiltonian matrix

[TABLE]

with $X,Z$ Pauli spin matrices, $\mathbb{1}$ the identity matrix and $(\alpha,\beta,\gamma)$ real numbers, we can write the exact time-propagator as

[TABLE]

with

[TABLE]

From this expression we can easily determine the needed angles $(\theta_{0},\dots,\theta_{3})$. The full circuit for sQPE is then

[TABLE]

and the last rotation can be avoided.