This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

\authorinfo

* : Send correspondence to Kwangmin Yu, E-mail : [email protected]

Practical numerical integration on NISQ devices

Kwangmin Yu Computational Science Initiative, Brookhaven National Laboratory, Upton, New York 11973, USA Hyunkyung Lim Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook, New York 11794, USA Pooja Rao Department of Mathematics, Stony Brook University, Stony Brook, New York 11794, USA
Abstract

This paper addresses the practical aspects of quantum algorithms used in numerical integration, specifically their implementation on Noisy Intermediate-Scale Quantum (NISQ) devices. Quantum algorithms for numerical integration utilize Quantum Amplitude Estimation (QAE) (Brassard et al., 2002) in conjunction with Grover’s algorithm. However, QAE is daunting to implement on NISQ devices since it typically relies on Quantum Phase Estimation (QPE), which requires many ancilla qubits and controlled operations. To mitigate these challenges, a recently published QAE algorithm (Suzuki et al., 2020), which does not rely on QPE, requires a much smaller number of controlled operations and does not require ancilla qubits. We implement this new algorithm for numerical integration on IBM quantum devices using Qiskit and optimize the circuit on each target device. We discuss the application of this algorithm on two qubits and its scalability to more than two qubits on NISQ devices.

keywords:
Quantum Algorithm, Quantum Counting, Monte Carlo Integration

1 Introduction

Since quantum computing (QC) was first introduced in the 1980s [1], it has grown into an active and diverse field of research. In recent years, significant progress has been made in building quantum computers by companies, such as IBM, Google, and Rigetti. Recently, Google, IBM, and Intel have announced 72 qubits [2], 50 qubits [3] (53 qubits [4]), and 49 qubits [5] quantum devices, respectively. More notably, IBM has made available its cloud enabled quantum computing platform to the public. These currently available NISQ devices provide a tangible quantum programming environment and are serving as a stepping stone for the large-scale universal quantum computers of the future [6]. Despite the significant breakthrough achieved in qubit numbers in the NISQ era, the quantum systems are limited by (a) the short coherence times, the amount of useful operational time in a calculation before the information loss, (b) the circuit depth, the number of sequential quantum operations that can be performed on a quantum device, and (c) the lack of error correction [6]. They are so noisy that the number of quantum gates that can be implemented is significantly impacted [7, 8]. Furthermore, due to limited physical inter-connectivity between the qubits, the multi-qubit gates (operations), such as controlled-NOT and controlled-controlled-NOT (Toffoli) are hard to implement efficiently [7, 9].

The aforementioned device limitations restrict the implementation of the quantum algorithms developed in the 1990s such as Quantum Fourier Transform (QFT) [10] and Quantum Phase Estimation (QPE)[11]. In particular, QPE, which requires many controlled operations and ancilla qubits, plays a critical role in various quantum algorithms, including Shor’s factoring algorithm[10]. This has motivated several efforts to design a more efficient version of QPE [12, 13], but the number of controlled operations is still restrictive. Therefore, the Quantum Amplitude Estimation (QAE) of Brassard et al. (canonical QAE) [14], which relies heavily on QPE, is daunting to implement on NISQ devices.

Several variants, such as maximum likelihood quantum amplitude estimation (MLQAE), iterative quantum amplitude estimation (IQAE), etc., have recently been published [15, 16, 17], that do not need QPE. Of particular importance to this study is the MLQAE method by Suzuki et al. [16], that has many desirable properties. First and foremost, since MLQAE does not use QPE, it avoids controlled operations in QPE and instead uses maximum likelihood estimation as a post-processing method. Another advantage is that it does not need ancilla qubits, which are necessary to achieve the desired accuracy in canonical QAE (CQAE). Finally, it is parallelizable because the queries can be executed simultaneously before they are post-processed using maximum likelihood estimation.

In this paper, we study the numerical implementation of quantum Monte Carlo integration (see Sec 1.1), whose building blocks are QAE and Grover’s search algorithm, [18] and its generalization, quantum amplitude amplification (QAA) [14]. The focus of the study is the investigation of methods for QAE that are appropriate for NISQ devices. To this end, we compare and contrast the performance of MLQAE with CQAE from a practical standpoint. We implement these quantum algorithms in Qiskit [19], an open source software suite for near-term quantum computing, and apply the transpiler in Qiskit to optimize the implementation on each target quantum device. This is followed by a comparative analysis of execution of these algorithms on IBMQ devices.

1.1 Monte Carlo Method for Integration

Monte Carlo (MC) methods are a wide class of computational techniques that make use of random number sampling. In this section, we review hit-or-miss Monte Carlo Integration (MCI) to estimate a definite integral numerically.

Given a domain 𝒟n\mathscr{D}\subset\mathbb{R}^{n} with two real valued functions, f,g:𝒟nf,g:\mathscr{D}\subset\mathbb{R}^{n}\rightarrow\mathbb{R} (typically g0g\equiv 0) and f>gf>g, consider the definite integral of fgf-g. Hit-or-miss MCI starts by generating samples uniformly from 𝒟×[an+1,bn+1]\mathscr{D}\times[a_{n+1},b_{n+1}] where an+1min(fg)a_{n+1}\leq\textrm{min}(f-g) and bn+1max(fg)b_{n+1}\geq\textrm{max}(f-g). A sample X=(x1,x2,,xn,xn+1)𝒟×[an+1,bn+1]X=(x_{1},x_{2},\cdots,x_{n},x_{n+1})\in\mathscr{D}\times[a_{n+1},b_{n+1}] is considered a hit if it satisfies xn+1<f(x1,x2,,xn)x_{n+1}<f(x_{1},x_{2},\cdots,x_{n}) and xn+1>g(x1,x2,,xn)x_{n+1}>g(x_{1},x_{2},\cdots,x_{n}). This method can be formulated as follows:

∫⋯∫𝒟(fg)𝑑xnVolume(𝒟×[an+1,bn+1]))#of hits#of samples.\frac{\idotsint_{\mathscr{D}}(f-g)d\textbf{x}^{n}}{\textrm{Volume}(\mathscr{D}\times[a_{n+1},b_{n+1}]))}\approx\frac{\#\textrm{of hits}}{\#\textrm{of samples}}.

The main advantage of MCI is that the error is solely a function of the number of the samples taken and so the basic process does not depend on the dimension of the function and thus can be scaled up easily, both in terms of effort and complexity, to higher dimensions. Deterministic numerical integration methods, such as the trapezoidal rule and Simpson’s rule, have error bounds that increase exponentially as the dimension of the integral increases. All these factors make MCI an appealing method in higher dimension (3\gg 3). Figure 1 shows a pictorial representation of MCI of 01sin(πx)𝑑x=2π0.63662\int_{0}^{1}\sin(\pi x)\,dx=\frac{2}{\pi}\approx 0.63662. As shown in Fig. 1, we have 0.570.57 from 100100 samples (left figure) and 0.640.64 from 10001000 samples.

Refer to caption
(a)
Refer to caption
(b)
Figure 1: A pictorial representation of Monte Carlo integration for sin(πx)\sin(\pi x) on the interval [0,1][0,1] using 100100 (5757 sample points fall inside the region bounded by the curves, i.e, 5757 hits) and 10001000 (640640 hits) uniformly sampled points, on the left and right hand side plots.

1.2 Quantum Monte Carlo Integration

Equation (1.1) shows that MCI is evaluated by counting the number of hits when the volume of 𝒟=𝒟×[an+1,bn+1]\mathscr{D^{\prime}}=\mathscr{D}\times[a_{n+1},b_{n+1}], is easy to compute. If 𝒟\mathscr{D^{\prime}} is complicated, then it can be divided into finite sub-domains that are simpler and the same procedure can be followed on each of the subdomains. MCI basically counts the number of good states ({x|χ(x)=1}\{\textbf{x}~{}|~{}\chi(\textbf{x})=1\}) of the boolean function, χ:𝒟n+1{0,1}\chi:\mathscr{D^{\prime}}\subset\mathbb{R}^{n+1}\rightarrow\{0,1\}, which depends on f,gf,g and domain 𝒟\mathscr{D^{\prime}} in Eq. (1.1).

One way of implementing MCI on a quantum computer is by counting the solutions of Grover’s search algorithm. In fact, CQAE implements MCI by using the generalized Grover’s algorithm (Quantum Amplitude Amplification) and QPE [14].

1.3 Notations

We follow Dirac’s Bra-Ket notation for qubit representation and arithmetic, such as the inner product and the tensor product. For a multiple qubit system, we use the consecutive binary number strings with the most significant qubit located on the left and the least significant qubit on the right. For example, we have |0|0|0=|0|0|0=|000\ket{0}\otimes\ket{0}\otimes\ket{0}=\ket{0}\ket{0}\ket{0}=\ket{000} in the binary representation and |011=|3\ket{011}=\ket{3} and |110=|6\ket{110}=\ket{6} in the decimal representation of the computational basis. In the decimal representation, the number of qubits nn is denoted by a subscript as in |0n\ket{0}_{n}. The dimension of a square matrix is also denoted by a subscript. For example, 𝕀n\mathbb{I}_{n} denotes the n×nn\times n identity matrix. In quantum circuit diagrams, the top and the bottom qubits represent the least and the most significant qubits, respectively. The three Pauli matrices (Pauli gates when they are used in quantum circuits) XX, YY and ZZ are as follows:

X=(0110),Y=(0ii0),Z=(1001).X=\begin{pmatrix}0&1\\ 1&0\end{pmatrix},~{}~{}Y=\begin{pmatrix}0&-i\\ i&0\end{pmatrix},~{}~{}Z=\begin{pmatrix}1&0\\ 0&-1\end{pmatrix}.

The Hadamard matrix, H = 12(1111)\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\ 1&-1\end{pmatrix}, is a unitary matrix that is also Hermitian, so it is its own inverse.

2 Quantum Amplitude Estimation Algorithms

First, we briefly review quantum amplitude amplification and amplitude estimation algorithms [14]. The amplitude amplification is a generalization of Grover’s search algorithm [18], without the loss of the quadratic quantum speedup over its classical counterpart.

Suppose we have a boolean function ff and a unitary operator 𝒜\mathcal{A}, where ff maps the good states to 11 and the bad states to 0 on domain 𝒟\mathscr{D} such that |𝒟|=N=2n|\mathscr{D}|=N=2^{n}, and 𝒜\mathcal{A} acts on n+1n+1 qubits such that

|Ψ=𝒜|0n|0=1a|ψ0n|0+a|ψ1n|1,\ket{\Psi}=\mathcal{A}\ket{0}_{n}\ket{0}=\sqrt{1-a}\ket{\psi_{0}}_{n}\ket{0}+\sqrt{a}\ket{\psi_{1}}_{n}\ket{1}, (1)

where the good state is |ψ1n\ket{\psi_{1}}_{n} with ||ψ1n|=k|\ket{\psi_{1}}_{n}|=k, the bad state is |ψ0n\ket{\psi_{0}}_{n}, and a=k/N[0,1]a=k/N\in[0,1] is unknown. The job of the quantum amplitude estimation algorithm is to find aa approximately. The algorithm complexity is measured by the number of quantum queries to the operator 𝒜\mathcal{A}.

To achieve the quantum speedup, instead of measuring the last qubit of |Ψ=𝒜|0n|0\ket{\Psi}=\mathcal{A}\ket{0}_{n}\ket{0} directly, |Ψ\ket{\Psi} is first amplified by the following unitary operator Q:

Q=𝒜S0𝒜1Sχ,\textbf{Q}=\mathcal{A}\textbf{S}_{0}\mathcal{A}^{-1}\textbf{S}_{\chi}, (2)

where S0=𝕀n+12|0n+10|n+1\textbf{S}_{0}=\mathbb{I}_{n+1}-2\ket{0}_{n+1}\bra{0}_{n+1} and Sχ=(n𝕀2)Z\textbf{S}_{\chi}=(\bigotimes\limits^{n}\mathbb{I}_{2})\otimes Z. Sχ\textbf{S}_{\chi} puts a negative sign to the good state |ψ1n|1\ket{\psi_{1}}_{n}\ket{1} and does nothing to the bad state |ψ0n|0\ket{\psi_{0}}_{n}\ket{0}. Let us define a parameter θ[0,π/2]\theta\in[0,\pi/2] so that sin2θ=a\sin^{2}\theta=a. With this, we can rewrite Eq. (1) as:

|Ψ=𝒜|0n|0=cosθ|ψ0n|0+sinθ|ψ1n|1.\ket{\Psi}=\mathcal{A}\ket{0}_{n}\ket{0}=\cos\theta\ket{\psi_{0}}_{n}\ket{0}+\sin\theta\ket{\psi_{1}}_{n}\ket{1}. (3)

By applying Q (amplitude amplification operator) repeatedly mm times on |Ψ\ket{\Psi}, we get

Qm|Ψ=cos((2m+1)θ)|ψ0n|0+sin((2m+1)θ)|ψ1n|1.\textbf{Q}^{m}\ket{\Psi}=\cos((2m+1)\theta)\ket{\psi_{0}}_{n}\ket{0}+\sin((2m+1)\theta)\ket{\psi_{1}}_{n}\ket{1}. (4)

From Eqs. (1) and (3), it can be observed that the measurement after applying Qm\textbf{Q}^{m} on 𝒜|0n|0\mathcal{A}\ket{0}_{n}\ket{0} shows a quadratically larger probability of obtaining the good state (provided θ\theta is sufficiently small so that (2m+1)θ<π2(2m+1)\theta<\frac{\pi}{2}) than measuring 𝒜|0n|0\mathcal{A}\ket{0}_{n}\ket{0} directly [14].

CQAE [14] estimates θ\theta in Eq. (3) by QPE which includes the inverse QFT. QPE is implemented by the controlled operation on the oracle queries. Hence, it needs a number of multi-controlled operations, which are further decomposed into many basis gates. When this algorithm is implemented on NISQ devices, the accuracy is strongly limited by the connectivity between qubits because all the ancilla registers need connectivity with the target register and a sufficiently large number of ancillae are needed to ensure desired accuracy of the estimation. On the other hand, MLQAE [16] is implemented by post-processing the result of the quantum computation (and measuring the process without QPE) using maximum likelihood estimation. Since different levels of amplification (depending on power of QQ) can be executed independently, the algorithm is parallelizable. The implementation aspects of both CQAE and MLQAE will be discussed in the following sections.

3 Implementation

For our study, a two-qubit domain is used to implement CQAE and MLQAE algorithms, and the solution is fixed to be |01\ket{01} out of the four possible values, |00\ket{00}, |01\ket{01}, |10\ket{10}, and |11\ket{11}. The algorithms are implemented in Qiskit [19], and optimized for IBM quantum devices using the transpiler in Qiskit.

3.1 Implementation

The quantum circuit implementation of Eqs. (1) and (2) with a two-qubit domain in Qiskit are described in Figs. 2 and 3, respectively. Figure 2 also includes the quantum circuit of 𝒜1\mathcal{A}^{-1} as it is used in Eq. (2).

Refer to caption
(a) 𝒜\mathcal{A}
Refer to caption
(b) 𝒜1\mathcal{A}^{-1}
Figure 2: Quantum circuit implementation of 𝒜\mathcal{A} from Eq. (1) in Qiskit.
Refer to caption
Figure 3: Quantum circuit implementation of 𝒜\mathcal{A} and Q from Eqs. (1) and (2), in Qiskit.

3.1.1 Canonical Quantum Amplitude Estimation

Refer to caption
Figure 4: A CQAE quantum circuit implementation of the two-qubit domain with m=1m=1. The circuit depth here is 7878.
Refer to caption
Figure 5: The optimization of the two-qubit canonical QAE circuit using the transpiler in Qiskit on IBMQ VIGO. The circuit depth of this optimized circuit is 195195.

CQAE needs additional ancilla qubits to read out the amplitude. The number of ancilla qubits is represented by mm in this paper. Figures 4 and 11 show the quantum circuits with mm = 11 and 33, respectively. Since the possible measured angle with mm ancillae are 2πi/2m2\pi i/2^{m} for i=0i=0 to 2m12^{m}-1, having more ancilla qubits increases the accuracy of the readout. For example, the possible measurements for m=1m=1 are 0 and π\pi. The quantum circuit in Fig. 4 is a variation of the circuit in Fig. 3. The former adds control on the Q operator followed by the inverse QFT.

To run a circuit on a quantum device, each quantum gate in the circuit should be decomposed into the basis gates supported by the device. Figure 5 shows the decomposed and optimized circuit for the IBMQ VIGO device. For optimization, the transpiler in Qiskit is used with the optimization level set at 33, the highest possible value allowed in the package. In spite of applying the optimization, the circuit depth increases by more than double because of decomposition of the gates into the basis gates, and the addition of SWAP gates to the circuit to supplement disconnection between the qubits.

For m=3m=3, the possible measurements of θ\theta are 0,π/4,π/2,,7π/40,\pi/4,\pi/2,\cdots,7\pi/4. The quantum circuit depth increases from 7878 (Fig. 4) to 514514 (Fig. 11) when mm is increased from 11 to 33 for more accuracy. In Fig. 11, a0a_{0}, a1a_{1} and a2a_{2} represent the ancilla qubits and the circuit is composed of 𝒜\mathcal{A}, followed by a0a_{0} controlled Q, followed by a1a_{1} controlled Q2\textbf{Q}^{2} (two consecutive a1a_{1} controlled Q), followed by a2a_{2} controlled Q4\textbf{Q}^{4} (four consecutive a2a_{2} controlled Q), and the inverse QFT at the end. The main reasons that lead to a sharp increase in the circuit depth are (a) exponential increase of controlled Q operations, which increase to 2m12^{m}-1 if mm ancillae are used, and (b) controlled operations imposed on Toffoli operations in 𝒜\mathcal{A}, Q, and 𝒜1,\mathcal{A}^{-1}, as shown in Fig. 3. The controlled Toffoli (ccc-NOT) are composed of a number of basic quantum gates as shown in Figs. 3 and 4. From Fig. 4 to Fig. 5 (m=1m=1 case), there is a twofold increases in circuit depth upon its decomposition into the basis gates for the quantum device. For m=3m=3, the circuit in Fig. 11, with a depth of 514514, should be decomposed into the basis quantum gates. As a consequence, its depth will become more than 10001000! And for that reason, we do not show the decomposed circuit in this paper. The implementation of such a circuit is not feasible on current NISQ devices because of decoherence and device error.

3.1.2 Maximum Likelihood Quantum Amplitude Estimation

The two main benefits of MLQAE over CQAE are (a) MLQAE does not need ancilla qubits to read out the amplitude a\sqrt{a} of Eq. (1), and (b) parallel execution. The analysis of error bounds and the number of queries (application of 𝒜\mathcal{A} and 𝒜1\mathcal{A}^{-1}) are discussed by Suzuki et al. [16] and Aaronson et al. [15]. Both CQAE and MLQAE need O(1ϵ(KN))O\bigl{(}\frac{1}{\epsilon}\bigl{(}\sqrt{\frac{K}{N}}\bigr{)}\bigr{)} oracle queries to estimate KK good states in NN samples (domain) with error ϵ\epsilon [15].

Figure 6 shows the quantum circuit implementation of MLQAE on a two-qubit domain. When only one circuit (only 𝒜\mathcal{A}) is used, then there is no quantum speed up [16]. More circuits increase the accuracy of the estimation of the amplitude. Suzuki et al. [16] discuss two options for circuit sequencing in MLQAE, linearly incremental sequence (LIS) and exponential incremental sequence (EIS), and suggest EIS as the asymptotically optimal choice. Thus, we adapt the EIS which has exponential power of Q from the second circuit. For example, the nn-th circuit has Q2n2\textbf{Q}^{2^{n-2}} operator after 𝒜\mathcal{A} when we have n>1n>1. Figure 6 shows four circuits of the MLQAE implementation and the last circuit (Fig. 6 (d)) has Q4\textbf{Q}^{4} operator.

The number of Q queries of MLQAE of nn circuits is equivalent to that of CQAE of n1n-1 ancilla qubits. For example, when we have four circuits for MLQAE, the number of Q queries is equivalent to that of CQAE having three ancilla qubits (m=3m=3 case). Since both CQAE and MLQAE have O(1ϵ(KN))O\bigl{(}\frac{1}{\epsilon}\bigl{(}\sqrt{\frac{K}{N}}\bigr{)}\bigr{)} oracle queries, we can compare the two circuits of Figs. 11 and 6 with the same accuracy. In this comparison, CQAE has a circuit depth of 514514 and MLQAE has a circuit depth of 9393 in total (sum of all circuit depths in Fig. 6). The advantage offered by MLQAE with respect to decoherence due to its smaller circuit depth is even more pronounced because the circuits that make up the MLQAE run independently, and the measurements are post-processed. For MLQAE, the noise is dominated by the last (longest) circuit due to decoherence and gate noise present in that circuit. The circuits (a), (b), (c), and (d) in Fig. 6 are named MLQAE[0]\left[0\right], MLQAE[1]\left[1\right], MLQAE[2]\left[2\right], MLQAE[3]\left[3\right], respectively.

Refer to caption
(a) 𝒜\mathcal{A}, circuit depth 44
Refer to caption
(b) Q𝒜\textbf{Q}\mathcal{A}, circuit depth 1515
Refer to caption
(c) Q2𝒜\textbf{Q}^{2}\mathcal{A}, circuit depth 2626
Refer to caption
(d) Q4𝒜\textbf{Q}^{4}\mathcal{A}, circuit depth 4848
Figure 6: An MLQAE quantum circuit implementation of two-qubit domain in Qiskit. We name circuits MLQAE[0]\left[0\right], MLQAE[1]\left[1\right], MLQAE[2]\left[2\right], MLQAE[3]\left[3\right] from (a) to (d), respectively.

4 Results

In this section, we discuss the results obtained by running MLQAE on IBMQ simulator and IBMQ devices. The simulator does not have quantum decoherence and gate errors. Thus, we can regard the result from the simulator as perfect, meaning devoid of any quantum errors, although it does have errors and inaccuracies which arise from the algorithm, such as the discretization error. Therefore, we can observe quantum noise or error from the implementation of the algorithm by comparing the results between the simulator and the devices.

4.1 Maximum Likelihood Estimation

The key idea of MLQAE is the post-processing of the measurements from each of the quantum circuits. The likelihood functions and the resultant maximum likelihood function are defined as following in the domain [0,π/2][0,\pi/2] for θ\theta:

Lk(hk;θ)={sin2((2mk+1)θ)}hk{cos2((2mk+1)θ)}Nhk,L_{k}(h_{k};\theta)=\{\sin^{2}((2m_{k}+1)\theta)\}^{h_{k}}\{\cos^{2}((2m_{k}+1)\theta)\}^{N-h_{k}}, (5)
L(h;θ)=k=0mLk(hk;θ),L(\vec{\textbf{h}};\theta)=\prod\limits_{k=0}^{m}L_{k}(h_{k};\theta), (6)

where mkm_{k}, hkh_{k}, and NN are the power of Q, hit count of 11, and number of shots for the kk-th circuit, respectively, and h=(h0,h1,,hm)\vec{\textbf{h}}=(h_{0},h_{1},\cdots,h_{m}). The maximum likelihood estimation estimates θ^\hat{\theta}, which maximizes L(h;θ^)L(\vec{\textbf{h}};\hat{\theta}) in the domain. But instead of L(h;θ)L(\vec{\textbf{h}};\theta), logL(h;θ)\log L(\vec{\textbf{h}};\theta) is used to estimate θ^\hat{\theta} since the log\log function is monotonically increasing. The Python code implementation is given in Appendix A.

4.2 Simulator

The results from the IBMQ simulator are shown in Fig. 7. Each histogram represents the result of its corresponding quantum circuit given in Fig. 6. Since aa is 1/41/4 and θ\theta is π/6\pi/6 in Eqs. (1) and (3), in this setting, 3θ3\theta, 5θ5\theta, and 9θ9\theta are π/2\pi/2, 5π/65\pi/6, and 3π/23\pi/2. Therefore, circuits MLQAE[0]\left[0\right], MLQAE[1]\left[1\right], MLQAE[2]\left[2\right], and MLQAE[3]\left[3\right] will collapse to |1\ket{1} with a probability of 0.250.25, 1.001.00, 0.250.25, 1.001.00, respectively, when they are measured in the states of Eq. (4). The results in Fig. 7 are consistent with the analytical computation.

To estimate θ\theta from the measured probability of each circuit, MLQAE uses the maximum likelihood estimation method [16] (see Appendix A). In this case, the estimated θ\theta is 0.5240.524 (see Eqs. (1) and (3)). Since the probability aa, of measuring |1\ket{1}, is sin2(θ)\sin^{2}(\theta), aa is 0.25040.250.2504\simeq 0.25, as expected.

Refer to caption
(a) 𝒜\mathcal{A}, 248248 hits
Refer to caption
(b) Q𝒜\textbf{Q}\mathcal{A}, 10241024 hits
Refer to caption
(c) Q2𝒜\textbf{Q}^{2}\mathcal{A}, 249249 hits
Refer to caption
(d) Q4𝒜\textbf{Q}^{4}\mathcal{A}, 10241024 hits
Figure 7: The results of the MLQAE (Fig. 6) on the IBMQ simulator with 10241024 shots.

4.3 IBM Quantum Devices

The optimizer (transpiler) maps logical qubits to device qubits optimally by considering qubit connectivity, and adds a SWAP gate if there is no physical connection between qubits. Figure 8 shows the layout of how qubits, such as q0q_{0}, q1q_{1}, and q2q_{2} in Fig. 6(d), are mapped to IBM quantum devices. For example, on IBMQX2, q0q_{0}, q1q_{1}, and q2q_{2} are mapped to 0, 11, and 22 as can be seen in Fig. 8.

Refer to caption
(a) IBMQX2 U2 error (5.452e-4, 1.007e-3) CNOT error (1.434e-2, 2.689e-2)
Refer to caption
(b) IBMQ VIGO U2 error (3.632e-4, 6.552e-4) CNOT error (7.252e-3, 1.108e-2)
Figure 8: The optimized layout of qubits for the circuits transpiled for IBMQX2 and IBMQ VIGO backends from Fig. 6(d), MLQAE[3]\left[3\right]. Errors shown are for a single qubit from the device calibration done on March 28, 2020.

The optimized circuits of the MLQAE implementation (see Fig. 6) for IBMQX2 and IBMQ VIGO are shown in Figs. 12 and 13, respectively. The main difference between the two devices is the qubit connectivity. The three qubits on IBMQX2 have inter-connectivity with one another while the qubits 11 and 22 on IBMQ VIGO are not connected. Thus, the optimized circuit for IBMQ VIGO, Fig 12 (b) in Appendix B, includes SWAP gates which are highlighted by red boxes. Since SWAP gates are decomposed into three CNOT gates, optimized circuits on IBMQ VIGO have longer circuit depth than the circuits for IBMQX2, as shown in Figs. 12 and 13 in Appendix B.

Refer to caption
(a) 𝒜\mathcal{A}, 468468 hits
Refer to caption
(b) Q𝒜\textbf{Q}\mathcal{A}, 738738 hits
Refer to caption
(c) Q2𝒜\textbf{Q}^{2}\mathcal{A}, 595595 hits
Refer to caption
(d) Q4𝒜\textbf{Q}^{4}\mathcal{A}, 667667 hits
Figure 9: The results of the MLQAE (Fig. 12) on IBMQX2 quantum device with 10241024 shots.
Refer to caption
(a) 𝒜\mathcal{A}, 274274 hits
Refer to caption
(b) Q𝒜\textbf{Q}\mathcal{A}, 712712 hits
Refer to caption
(c) Q2𝒜\textbf{Q}^{2}\mathcal{A}, 401401 hits
Refer to caption
(d) Q4𝒜\textbf{Q}^{4}\mathcal{A}, 589589 hits
Figure 10: The results of the MLQAE (Fig. 13) on VIGO quantum device with 10241024 shots.

Figures 9 and 10 are histograms of measurements on quantum devices. The results from IBMQ VIGO are more accurate than those from IBMQX2. The result from the MLQAE[0] circuit on IBMQX2 (Fig. 9 (a)) shows a relative error of 83.8%83.8\%, whereas that value on IBMQ VIGO (Fig. 10 (a)) is only 7.2%7.2\%. This is because IBMQ VIGO has more accurate quantum device operations than IBMQX2, as depicted in Fig. 8. For the θ\theta-estimation problem (see Appendix A), IBMQX2 and IBMQ VIGO predict θ=0.795\theta=0.795 (a=0.509a=0.509) and θ=0.780\theta=0.780 (a=0.494a=0.494), respectively. The former has 51.8%51.8\% and 103.6%103.6\% as the relative errors of θ\theta and aa and those values for the latter are 49.0%49.0\% and 97.6%97.6\%. In contrast with the result from the MLQAE[0] circuit, there is not a significant difference between the overall accuracy of the two devices. We hypothesize that IBMQ VIGO loses its accuracy because of the long circuit depths in MLQAE[2] and MLQAE[3].

5 Conclusion

In this paper, we have implemented and discussed two quantum amplitude estimation algorithms on IBM quantum devices, in the broader context of Monte Carlo integration. Even though CQAE [14] is a monumental algorithm and a brilliant extension of Grover’s algorithm [18], it is extremely challenging and in some cases infeasible to implement on NISQ devices, because of the required number of controlled operations acting on Q operators in QPE. Recently developed quantum amplitude estimation algorithms do not use QPE. We have implemented one such algorithm, MLQAE, along with CQAE on IBM quantum devices and discussed their applicability on NISQ devices. As discussed in Sec. 4, MLQAE is practical even on NISQ devices and has several advantages, including smaller circuit depth and parallel execution. The advancements in QAE for quantum Monte Carlo integration rely on quantum amplitude amplification (QAA).

As shown in Figs. 3 and 6, one of the crucial components of QAA is the S0\textbf{S}_{0} operator, and it needs a multi-controlled NOT operator. In general, when we have an nn-qubit problem domain (2n2^{n}), we need nn controlled NOT gates for each Q operator. For example, when we have a four-qubit domain (1616), the S0\textbf{S}_{0} gate needs the following quantum gate,

\Qcircuit@C=1.0em@R=0.2em@!R\lstickq0&\targ\qw\qw\lstickq1\ctrl1\qw\qw\lstickq2\ctrl1\qw\qw\lstickq3\ctrl1\qw\qw\lstickq4\ctrl1\qw\qw\Qcircuit@C=1.0em@R=0.2em@!R{\lstick{{q}_{0}}&\targ\qw\qw\\ \lstick{{q}_{1}}\ctrl{-1}\qw\qw\\ \lstick{{q}_{2}}\ctrl{-1}\qw\qw\\ \lstick{{q}_{3}}\ctrl{-1}\qw\qw\\ \lstick{{q}_{4}}\ctrl{-1}\qw\qw\\ }

The implementation of the ‘advanced’ mode of the above circuit in Qiskit is shown in Fig. 14. Even though the ‘basic’ mode has shorter circuit depth, it needs two additional ancilla qubits as shown in Fig. 15. The circuit decomposition of the ‘advanced’ mode for IBMQX2 and IBMQ VIGO are shown in Fig. 16 and 17, respectively. Since IBMQ VIGO has less connectivity than IBMQX2 (see Fig. 8), the circuit for IBMQ VIGO has much longer depth (181181) than that of IBMQX2 (142142). This example shows that the multi-controlled NOT gates have mainly two difficulties on NISQ devices. The first is that a multi-controlled NOT gate should be decomposed into many basis gates as shown in Figs. 16 and 17. The second difficulty is that all qubits which are involved in the gate should be connected. As shown in Fig. 17, the lack of connectivity between qubits are complemented by SWAP gates (decomposed into three CNOT gates), increasing decoherence through a longer circuit and increasing gate error. Therefore, it is necessary to develop more efficient circuits for S0\textbf{S}_{0}, or better algorithms, to avoid the problems S0\textbf{S}_{0} has in order to run scalable quantum Monte Carlo integration on NISQ devices.

References

  • [1] Feynman, R. P., “Simulating physics with computers,” International journal of theoretical physics 21(6), 467–488 (1982).
  • [2] Kelly, J., “A preview of bristlecone, google’s new quantum processor,” Google AI Blog (March 2018).
  • [3] Knight, W., “Ibm raises the bar with a 50-qubit quantum computer,” (November 2017).
  • [4] Shankland, S., “Ibm’s new 53-qubit quantum computer is its biggest yet,” CNET (September 2019).
  • [5] Hsu, J., “Ces 2018: Intel’s 49-qubit chip shoots for quantum supremacy,” IEEE Spectrum (January 2018).
  • [6] Preskill, J., “Quantum Computing in the NISQ era and beyond,” Quantum 2, 79 (Aug. 2018).
  • [7] Tannu, S. S. and Qureshi, M. K., “Not all qubits are created equal: a case for variability-aware policies for NISQ-era quantum computers,” in [Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems ], 987–999 (2019).
  • [8] Córcoles, A. D., Kandala, A., Javadi-Abhari, A., McClure, D. T., Cross, A. W., Temme, K., Nation, P. D., Steffen, M., and Gambetta, J., “Challenges and opportunities of near-term quantum computing systems,” arXiv preprint arXiv:1910.02894 (2019).
  • [9] Adbo, B. et al., “Ibm q 5 yorktown v1.x.x,” (January 2019).
  • [10] Shor, P. W., “Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer,” SIAM J. Comput. 26, 1484–1509 (Oct. 1997).
  • [11] Kitaev, A. Y., “Quantum measurements and the abelian stabilizer problem,” arXiv preprint quant-ph/9511026 (1995).
  • [12] Svore, K. M., Hastings, M. B., and Freedman, M., “Faster phase estimation,” Quantum Info. Comput. 14, 306–328 (Mar. 2014).
  • [13] van den Berg, E., “Practical sampling schemes for quantum phase estimation,” (2019).
  • [14] Brassard, G., Hoyer, P., Mosca, M., and Tapp, A., “Quantum amplitude amplification and estimation,” Contemporary Mathematics 305, 53–74 (2002).
  • [15] Aaronson, S. and Rall, P., “Quantum approximate counting, simplified,” Symposium on Simplicity in Algorithms , 24–32 (Jan 2020).
  • [16] Suzuki, Y., Uno, S., Raymond, R., Tanaka, T., Onodera, T., and Yamamoto, N., “Amplitude estimation without phase estimation,” Quantum Information Processing 19(2), 75 (2020).
  • [17] Grinko, D., Gacon, J., Zoufal, C., and Woerner, S., “Iterative quantum amplitude estimation,” arXiv preprint arXiv:1912.05559 (2019).
  • [18] Grover, L. K., “A fast quantum mechanical algorithm for database search,” Proceedings of the twenty-eighth annual ACM symposium on Theory of computing - STOC ’96 (1996).
  • [19] Aleksandrowicz, G., Alexander, T., Barkoutsos, P., Bello, L., Ben-Haim, Y., Bucher, D., Cabrera-Hernández, F. J., Carballo-Franquis, J., Chen, A., Chen, C.-F., and et al., “”qiskit: An open-source framework for quantum computing”,” (Jan 2019).

Appendix A Maximum Likelihood Estimation Implementation

  

1
2import numpy as np
3from scipy.optimize import brute
4
5def MaximumLikelihoodEstmator(circuit_length, ones, zeros):
6
7 grid = 20000
8 epsilon = 1/grid
9 domain = [0.0 + epsilon, np.pi/2 - epsilon] # to avoid zero
10
11 def logL(theta):
12 fval = 0
13 for i in range(circuit_length):
14 if i==0:
15 fval += 2 * ones[i] * np.log( np.absolute(np.sin(theta)) )
16 fval += 2 * zeros[i] * np.log( np.absolute(np.cos(theta)) )
17 else:
18 fval += 2*ones[i]*np.log(np.absolute(np.sin((2*(2**(i-1))+1)*theta)))
19 fval += 2*zeros[i]*np.log(np.absolute(np.cos((2*(2**(i-1))+1)*theta)))
20 return -fval # to compute maximum
21
22 return brute(logL, [domain], Ns=grid)[0]
23
24
25
26ones_sim = [248, 1024, 249, 1024]
27zeros_sim = [776, 0, 775, 0]
28
29ones_ibmqx2 = [468, 738, 595, 667]
30zeros_ibmqx2= [556, 286, 429, 357]
31
32ones_vigo = [274, 712, 401, 589]
33zeros_vigo = [750, 312, 623, 435]
34
35est_theta_sim = MaximumLikelihoodEstmator(4, ones_sim, zeros_sim)
36est_prob_sim = np.sin(est_theta_sim)**2
37
38print(’Estimated theta Simulator : , est_theta_sim)
39print(’Estimated probability Simulator : , est_prob_sim)
40
41est_theta_ibmqx2 = MaximumLikelihoodEstmator(4, ones_ibmqx2, zeros_ibmqx2)
42est_prob_ibmqx2 = np.sin(est_theta_ibmqx2)**2
43
44print(’Estimated theta IBMQX2 : , est_theta_ibmqx2)
45print(’Estimated probability IBMQX2 : , est_prob_ibmqx2)
46
47est_theta_vigo = MaximumLikelihoodEstmator(4, ones_vigo, zeros_vigo)
48est_prob_vigo = np.sin(est_theta_vigo)**2
49
50print(’Estimated theta IBMQ VIGO : , est_theta_vigo)
51print(’Estimated probability IBMQ VIGO : , est_prob_vigo)
Listing 1: The Maximum Likelihood Estimator Implementation

Appendix B Quantum Circuit Diagrams

Refer to caption
Figure 11: A CQAE quantum circuit implementation of a two-qubit domain for m=3m=3. The circuit depth is 514514.
Refer to caption
(a) 𝒜\mathcal{A}, circuit depth 1010 before the measurement
Refer to caption
(b) Q𝒜\textbf{Q}\mathcal{A}, circuit depth 4242 before the measurement
Refer to caption
(c) Q2𝒜\textbf{Q}^{2}\mathcal{A}, circuit depth 7474 before the measurement
Figure 12: The optimization of the two-qubit MLQAE circuit using the transpiler in Qiskit on IBMQX2 device. The optimization of Fig. 6 (d) is omitted because the circuit length is too long.
Refer to caption
(a) 𝒜\mathcal{A}, circuit depth 1313 before the measurement
Refer to caption
(b) Q𝒜\textbf{Q}\mathcal{A}, circuit depth 7272 before the measurement
Refer to caption
(c) Q2𝒜\textbf{Q}^{2}\mathcal{A}, circuit depth 121121 before the measurement
Figure 13: The optimization of the two-qubit MLQAE circuit using the transpiler in Qiskit on IMBQ VIGO device. The optimization of Fig. 6 (d) is omitted because the circuit length is too long. In (b), SWAP gates are highlighted by red boxes.
Refer to caption
Figure 14: The cccc-NOT implementation in ‘advanced’ mode in Qiskit. The circuit depth is 6262.
Refer to caption
Figure 15: The cccc-NOT implementation in ‘basic’ mode in Qiskit. This mode needs two additional ancilla qubits. The circuit depth is 3333.
Refer to caption
Figure 16: The cccc-NOT implementation on IBMQX2 in ‘advanced’ mode. The circuit depth is 142142.
Refer to caption
Figure 17: The cccc-NOT implementation on IBMQ VIGO in ‘advanced’ mode. The circuit depth is 181181.