Iterative phase estimation

Note

Click here to view this tutorial as an interactive Jupyter notebook in IBM Quantum Lab (requires sign-in).

The goal of this exercise is to understand how the Iterative Phase Estimation (IPE) algorithm uses dynamic circuits, and why would we use the IPE algorithm instead of the QPE (Quantum Phase Estimation) algorithm. We demonstrate how to perform this algorithm with Qiskit using circuits containing reset in order to apply gates conditioned by the values resulting from previous measurements and stored in a classical register.

If you want to learn more about iterative phase estimation, visit the references below.

References

Preamble

import os
from typing import Any, List, Dict, Union

import numpy as np
import matplotlib.pyplot as plt

from qiskit import IBMQ, Aer, QuantumCircuit, QuantumRegister, ClassicalRegister, transpile, execute
from qiskit.tools.visualization import plot_histogram
from qiskit.result import marginal_counts
import qiskit.tools.jupyter

from run_openqasm3 import run_openqasm3

%matplotlib inline

import warnings
warnings.filterwarnings("ignore")

pi = np.pi
# Note: This can be any hub/group/project that has access to the required device and the Qiskit runtime.
# Verify that ``qasm3`` is present in ``backend.configuration().supported_features``.
# hub = "<hub>"
# group = "<group>"
# project = "<project>"
# backend_name = "<your backend>"
IBMQ.load_account()
provider = IBMQ.get_provider(hub=hub, group=group, project=project)
backend_real = provider.get_backend(backend_name)
backend_sim = Aer.get_backend('aer_simulator')
config = backend_real.configuration()
defaults = backend_real.defaults()

basis_gates = config.basis_gates
# To get more information about your backend
# import qiskit.tools.jupyter
backend_real
from qiskit import qasm3

def dump_qasm3(circuit, backend=backend_real):
    return qasm3.Exporter(includes=[], basis_gates=basis_gates+["reset"], disable_constants=True).dumps(circuit)

Conditioned gates: the c_if method

Before starting the IPE algorithm, let’s understand the Qiskit conditional method c_if, as it goes into building the IPE circuit.

c_if is a function (actually a method of the gate class) to perform conditioned operations based on the value stored previously in a classical register. With this feature you can apply gates after a measurement in the same circuit conditioned by the measurement outcome.

For example, the following code will execute the X gate if the value of the classical register is 0 .

q = QuantumRegister(1,'q')
c = ClassicalRegister(1,'c')
qc = QuantumCircuit(q, c)
qc.h(0)
qc.measure(0,0)
qc.x(0).c_if(c, 0)
qc.draw(output='mpl')
../../../../_images/output_12_03.png

We highlight that the method c_if expects as the first argument a whole classical register, not a single classical bit (or a list of classical bits), and as the second argument a value in decimal representation (a non-negative integer), not the value of a single bit, 0 or 1 (or a list/string of binary digits).

Let’s make another example. Consider that we want to perform a bit flip on the third qubit after the measurements in the following circuit, when the results of the measurement of q_0 and q_1 are both 1 .

q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
qc = QuantumCircuit(q, c)
qc.h(q[0])
qc.h(q[1])
qc.h(q[2])
qc.barrier()
qc.measure(q,c)
qc.draw('mpl')
../../../../_images/output_14_03.png

We want to apply the X gate, only if both the results of the measurement of q_1 and q_2 are 1 . We can do this using the c_if method, conditioning the application of X depending on the value passed as argument to c_if.

We will have to encode the value to pass to the c_if method such that it will check the values 011 and 111 (in binary representation), since it does not matter what is in the rightmost position.

The two-integer values in decimal representation:

../../../../_images/binary3.png

We can check the solutions using the bin() method in Python (the prefix 0b indicates the binary format).

print(bin(3))
print(bin(7))
0b11
0b111

We apply X to q_2 using c_if two times, one for each value corresponding to 011 and 111 .

q = QuantumRegister(3,'q')
c = ClassicalRegister(3,'c')
qc = QuantumCircuit(q, c)
qc.h(0)
qc.h(1)
qc.h(2)
qc.barrier()
qc.measure(q,c)

qc.x(2).c_if(c, 3)  # for the 011 case
qc.x(2).c_if(c, 7)  # for the 111 case

qc.draw(output='mpl')
../../../../_images/output_20_03.png

Iterative phase estimation versus quantum phase estimation

The motivation for using the IPE algorithm is that while the Quantum Phase Estimation algorithm works fine for short-depth circuits, it doesn’t work properly once the circuit starts to grow, due to gate noise and decoherence times.

IPE example with a 1-qubit gate for U

We want to apply the IPE algorithm to estimate the phase for a 1-qubit operator U . For example, here we use the S -gate.

Let’s apply the IPE algorithm to estimate the phase for the S -gate. Its matrix is

\begin{split}S =\begin{bmatrix} 1 & 0 \\ 0 & e^\frac{i\pi}{2} \end{bmatrix}\end{split}

That is, the S -gate adds a phase \pi/2 to the state |1\rangle , leaving unchanged the phase of the state |0\rangle

S|1\rangle = e^\frac{i\pi}{2}|1\rangle

In the following, we will use the notation and terms used in Section 2 of lab 4 of the Qiskit Textbook.

To estimate the phase \phi=\frac{\pi}{2} for the eigenstate |1\rangle , we should find \varphi=\frac{1}{4} (where \phi = 2 \pi \varphi ). Therefore, to estimate the phase we need exactly two phase bits, i.e.,  m=2 , since 1/2^2=1/4 . Therefore, \varphi=0.\varphi_1\varphi_2 .

Remember from the theory that for the IPE algorithm, m is also the number of iterations, so we need only 2 iterations or steps.

First, we initialize the circuit. IPE works with only one auxiliary qubit, instead of m counting qubits of the QPE algorithm. Therefore, we need two qubits: one auxiliary qubit, and one for the eigenstate of U -gate, as well as a classical register of two bits, for the phase bits \varphi_1 , \varphi_2 .

IPE - First step

Now we build the quantum circuit for the first step (that is, the first iteration of the algorithm), to estimate the least significant phase bit \varphi_m , which in this case is \varphi_2 . For the first step we have four sub-steps:

  • qubit initialization

  • initialization

  • application of the Controlled- U gates

  • measure of the auxiliary qubit in the x-basis

qubits = [0, 1]
shots = 1000

nq = 2
m = 2
qr = QuantumRegister(2,'q')
c0 = ClassicalRegister(1,'c0')
c1 = ClassicalRegister(1, 'c1')
qc_S = QuantumCircuit(qr, c0, c1)

Eigenstate initialization

The initialization consists of applying the Hadamard gate to the auxiliary qubit, and the preparation of the eigenstate |1\rangle .

def initialize_eigenstate(qc):
    """Initialize the eigenstate and prepare our ancilla qubit"""
    qc.h(0)
    qc.x(1)

initialize_eigenstate(qc_S)
qc_S.draw('mpl')
../../../../_images/output_26_03.png

Application of the Controlled- U gates

We now have to apply 2^t times the Controlled- U operators (see also in the Qiskit documentation Two-qubit gates), that, in this example, is the Controlled- S gate ( CS for short).

To implement CS in the circuit, since S is a phase gate, we can use the controlled phase gate \text{CP}(\theta) , with \theta=\pi/2 .

theta = 1 * np.pi / 2
cu_circ = QuantumCircuit(2)
cu_circ.cp(theta,0,1)
cu_circ.draw('mpl')
../../../../_images/output_29_03.png

Let’s apply 2^t times \text{CP}(\pi/2) . Since for the first step t=m-1 , and m=2 , we have 2^t=2 .

for _ in range(2**(m-1)):
    qc_S.cp(theta,0,1)
qc_S.draw('mpl')
../../../../_images/output_31_03.png

Measure in the x-basis

Finally, we perform the measurement of the auxiliary qubit in the x-basis.

In this way we obtain the phase bit \varphi_2 and store it in the classical bit c_0 .

Note

The eigenvectors of the Hadamard gate are the basis vectors of the x-basis.

def x_measurement(qc, qubit, cbit):
    """Measure 'qubit' in the X-basis, and store the result in 'cbit'"""
    qc.h(qubit)
    qc.measure(qubit, cbit)

x_measurement(qc_S, qr[0], c0)
qc_S.draw('mpl')
../../../../_images/output_33_03.png

Subsequent steps

Now we build the quantum circuit for the other remaining steps (in this example, only the second one). In these steps we have four sub-steps: the three sub-steps as in the first step and, in the middle, the additional step of the phase correction.

  • initialization with reset

  • phase correction

  • application of the Control- U gates

  • measurement of the auxiliary qubit in x-basis

Reset the auxiliary qubit

As we want to perform an iterative algorithm in the same circuit, we need to reset the auxiliary qubit q0 after the measurement gate and initialize it again as before, to recycle the qubit.

We use the built-in qc.reset operation to perform a reset of qubit 0, followed by the Hadamard operation.

def reset_auxiliary(qc):
    qc.reset(0)
    qc.h(0)

reset_auxiliary(qc_S)
qc_S.draw('mpl')
../../../../_images/output_36_03.png

Phase correction (for step 2)

As seen in the theory, in order to extract the phase bit \varphi_{1} , we perform a phase correction of -\pi\varphi_2/2 . Of course, we need to apply the phase correction in the circuit only if the phase bit \varphi_2=1 , i.e., we have to apply the phase correction of -\pi/2 only if the classical bit c_0 is 1.

Therefore, after the reset we apply the phase gate P(\theta) with phase \theta=-\pi/2 conditioned by the classical bit c_0 ( =\varphi_2 ) using the c_if method. As we saw in the first part of this tutorial, we have to use the c_if method with a value of 1, as 1_{10} = 001_{2} (the subscripts _{10} and _2 indicate the decimal and binary representations).

qc_S.p(-np.pi/2,0).c_if(c0, 1)
qc_S.draw('mpl')
../../../../_images/output_38_03.png

Apply the controlled- U gates and x-measurement (for step 2)

We apply the CU operations as we did in the first step. For the second step we have t=m-2 , hence 2^t=1 . We therefore apply \text{CP}(\pi/2) once, and then we perform the x-measurement of the qubit q_0 , storing the result, the phase bit \varphi_1 , in bit c_1 of the classical register.

## 2^t c-U operations (with t=m-2)
for _ in range(2**(m-2)):
    qc_S.cp(theta,0,1)

x_measurement(qc_S, qr[0], c1)

We now have our final circuit!

qc_S.draw('mpl')
../../../../_images/output_42_03.png

Let’s execute the circuit with the qasm_simulator, the simulator without noise that runs locally.

count = execute(qc_S, backend_sim).result().get_counts()
count
{'0 1': 1024}
count0 = execute(qc_S, backend_sim).result().get_counts()
plot_histogram(count0)
../../../../_images/output_45_03.png
count0
{'0 1': 1024}

In the picture we have the same histograms, but on the left we have on the x-axis the string with phase bits \varphi_1 , \varphi_2 , and on the right the actual phase \varphi in decimal representation.

As we expected, we have found \varphi=\frac{1}{4}=0.25 with a 100\% probability.

Run on real hardware

qc_S_backend = transpile(qc_S, backend_real, initial_layout=qubits, optimization_level=3)
qc_S_backend.draw(output="mpl")
../../../../_images/output_49_03.png
real_ipe_job = run_openqasm3(qc_S_backend, backend_real, verbose=False, shots=shots)
print(f"IPE job id: {real_ipe_job.job_id()}")
IPE job id: caosq7nhrqeo3qj6m3og
real_counts0 = real_ipe_job.result().get_counts()
plot_histogram(real_counts0)
../../../../_images/output_51_03.png

Exercise 1 - Try estimating the phase of another gate

Now that we have estimated the phase of the S gate, let’s try finding the phase of another gate. Recall that we implemented IPE for S using the P(\theta) gate. Choose a new instance of the P gate that can be estimated with two bits of precision and rerun.

Click here for the solution

Modify \theta and rerun

Exercise 2 - Adding an extra bit of precision

Now that we have mastered IPE for two bits of precision, let’s try adding a third bit. Modify the code of this exercise to add another round of IPE.

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit SoftwareVersion
qiskit-terra0.21.0
qiskit-aer0.10.4
qiskit-ignis0.7.0
qiskit-ibmq-provider0.20.0.dev0+4f1f8c6
qiskit0.36.1
System information
Python version3.9.5
Python compilerClang 10.0.0
Python builddefault, May 18 2021 12:31:01
OSDarwin
CPUs8
Memory (Gb)32.0
Tue Jun 21 09:52:47 2022 EDT

This code is a part of Qiskit

© Copyright IBM 2017, 2022.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.