Lab06_IterativePhaseEstimation.ipynb

Kernel: Python 3

Lab 6 Iterative Phase Estimation Algorithm

Prerequisite

• Ch.3.5 Quantum Fourier Transform

• Ch.3.6 Quantum Phase Estimation

• Ch.1.4 Single Qubit Gates

• Summary of Quantum Operations

Other relevant materials

• Device backend noise model simulations

• Hellinger fidelity

from qiskit import *
import numpy as np
from qiskit.visualization import plot_histogram
import qiskit.tools.jupyter
import matplotlib.pyplot as plt

sim = Aer.get_backend('qasm_simulator')

Part 1: Implementation of Iterative Phase Estimation algorithm

Goal

Estimate a phase value on a system of two qubits through Iterative Phase Estimation (IPE) algorithm.

Having gone through previous labs, you should have noticed that the "length" of a quantum circuit is the primary factor when determining the magnitude of the errors in the resulting output distribution; quantum circuits with greater depth have decreased fidelity. Therefore, implementing algorithms based on shallow depth circuits is of the great importance in near-term quantum computing. In Lab 5, we learn one such algorithm for estimating quantum phase called the Iterative Phase Estimation (IPE) algorithm which requires a system comprised of only a single auxiliary qubit and evaluate the phase through a repetative process.

1. Understand a circuit with non-unitary operations.

Before we learn how the IPE algorithm works, lets review reset and conditional operations in Qiskit that go into building a IPE circuit. Read the Qiskit tutorial here ( go to Non-unitary operations section ) to understand how to build a circuit that performs conditional operations and reset.

📓Step A. Run the following cell and predict the outcome of the circuit.

q = QuantumRegister(1)
c = ClassicalRegister(2)

qc0 = QuantumCircuit(q, c)
qc0.h(q[0])
qc0.measure(q[0], c[0])
qc0.reset(q[0])
qc0.h(q[0])
qc0.p(np.pi/3, q[0]).c_if(c,1)
qc0.h(q[0])
qc0.measure(q[0],c[1])

qc0.draw()

Execute the cell below to see the result.

count0 = execute(qc0, sim).result().get_counts()
plot_histogram(count0)

📓Step B. Complete the rest of the circuit so that the auxiliary qubit ( top qubit ) after the reset would be in the state $\frac{1}{\sqrt2}(|0\rangle + e^{-i\frac{\pi}{2}}|1\rangle)$ if the value of the classical bit is one or remains zero state otherwise.

q = QuantumRegister(2)
c = ClassicalRegister(2)

qc1 = QuantumCircuit(q,c)
qc1.h(q[0])
qc1.x(q[1])
qc1.cp(np.pi/5, q[0], q[1])
qc1.measure(q[0], c[0])
qc1.reset(q[0])

###### your code goes here #####





##########################
qc1.h(q[0])
qc1.measure(q[0],c[1])

qc1.draw()

Running the following cell to display the result.

count1 = execute(qc1, sim).result().get_counts()
plot_histogram(count1)

2. Iterative Phase Estimation (IPE) Algorithm.

The Quantum Phase Estimation (QPE) circuit that we have learned and used previously is limited by the number of qubits necessary for the algorithm’s precision. Every additional qubit has added costs in terms of noise and hardware requirements; noisy results that we obtained from executing the QPE circuit on a real quantum device in Lab 4 would get worse as the number of the qubits on the circuit increases. The IPE algorithm implements quantum phase estimation with only a single auxiliary qubit, and the accuracy of the algorithm is restricted by the number of iterations rather than the number of counting qubits. Therefore, IPE circuits are of practical interest and are of foremost importance for near-term quantum computing as QPE is an essential element in many quantum algorithms.

Consider the problem of finding $\varphi$ given $|\Psi\rangle$ and $U$ such that $U |\Psi\rangle = e^{i \phi} | \Psi \rangle$, with $\phi = 2 \pi \varphi$. Let's assume for now that $\varphi$ can be written as $\varphi = \varphi_1/2 + \varphi_2/4 + ... + \varphi_m/2^m = 0.\varphi_1 \varphi_2 ... \varphi_m$, where we have previously defined the notation $0.\varphi_1 \varphi_2 ... \varphi_m$.

Assume that $U$ is a unitary operator acting on one qubit. We therefore need a system of two qubits, $q_0$ and $q_1$, where $q_0$ is auxiliary qubit and the qubit $q_1$ represents the physical system on which $U$ operates. Having them initialized as $q_0 \rightarrow |+\rangle$ and $q_1 \rightarrow |\Psi \rangle$, application of control-U between $q_0$ and $q_1$ $2^t$ times would change the state of $q_0$ to $|0\rangle + e^{i 2 \pi 2^{t} \varphi} | 1 \rangle$. That is, the phase of $U$ has been kicked back into $q_0$ as many times as the control operation has been performed.

Therefore,

for $t=0$, the phase encoded into $q_0$ would be $e^{i 2 \pi 2^{0} \varphi} = e^{i 2 \pi \varphi} = e^{i 2 \pi 0.\varphi_1 \varphi_2 ... \varphi_m}$ and

for $t=1$, the phase would be $e^{i 2 \pi 2^{1} \varphi} = e^{i 2 \pi \varphi_1} e^{i 2 \pi 0.\varphi_2 \varphi_3 ... \varphi_m}$ and

for $t=2$, $e^{i 2 \pi 2^{2} \varphi} = e^{i 2 \pi 2 \varphi_1} e^{i 2 \pi \varphi_2} e^{i 2 \pi 0.\varphi_3 \varphi_4 ... \varphi_m}$ and

for $t=m-1$, $e^{i 2 \pi 2^{m-1} \varphi} = e^{i 2 \pi 2^{m-2} \varphi_1} e^{i 2 \pi 2^{m-3} \varphi_2} ... e^{i 2 \pi 2^{-1} \varphi_m} = e^{i 2 \pi 0.\varphi_m}$.

Note that for the last case with $t=m-1$, the state of $q_0$ is $|0\rangle + e^{i 2 \pi 0.\varphi_m}|1\rangle$; $|+\rangle$ if $\varphi_m = 0$ and $|-\rangle$ if $\varphi_m = 1$ which would produce outcomes $|0\rangle$ and $|1\rangle$ respectively when it gets measured in x-basis.

In the first step of the IPE algorithm, we directly measure the least significant bit of the phase $\varphi$, $\varphi_m$, by initializing the 2-qubit registers as described above ( $q_0 \rightarrow |+\rangle$ and $q_1 \rightarrow |\Psi \rangle$ ), performing $2^{m-1}$ control-$U$ operations between the qubits, and measuring $q_0$ in the x-basis.

For the second step, we initialize the systems in the same way and apply $2^{m-2}$ control-$U$ operations. The relative phase in $q_0$ after these operations is now $e^{i 2 \pi 0.\varphi_{m-1}\varphi_{m}}= e^{i 2 \pi 0.\varphi_{m-1}} e^{i 2 \pi \varphi_m/4}$. To extract the phase bit $\varphi_{m-1}$, first perform a phase correction by rotating around the $Z-$axis of angle $-2 \pi \varphi_m/4=-\pi \varphi_m/2$, which results in the state of $q_0$ to be $|0\rangle + e^{i 2 \pi 0.\varphi_{m-1}} | 1 \rangle$. Perform a measurement on $q_0$ in x-basis to obtain the phase bit $\varphi_{m-1}$.

Therefore, the $k$th step of the IPE, getting $\varphi_{m-k+1}$, consists of the register initialization ($q_0$ in $|+\rangle$, $q_1$ in $|\Psi\rangle$), the application of control-$U$ $2^{m-k}$ times, a rotation around $Z$ of angle $\omega_k = -2 \pi 0.0\varphi_{m-k+2} ... \varphi_m$, and a measurement of $q_0$ in x-basis: a Hadamard transform to $q_0$, and a measurement of $q_0$ in the standard basis. Note that $q_1$ remains in the state $|\Psi\rangle$ throughout the algorithm.

3. Estimate the phase of the $T$-gate implementing IPE algorithm.

Review the section 2. Example: T-gate in Ch.3.8 Quantum Phase Estimation and the section 4. Measuring in Different Bases in Ch.1.4 Single Qubit Gates

As we already learned the Qiskit textbook, the phase of a T-gate is exactly expressed using three bits.

📓Step A. Obtain the least significant phase bit of the $T$-gate by setting up the circuit T_x3 properly and assign the value to the variable x_3.

In the previous section, the first step explains how to construct the circuit to extract the least significant phase bit.

q = QuantumRegister(2)
c = ClassicalRegister(1)

T_x3 = QuantumCircuit(q,c)

########## your code goes here #######

##1 Initialization





##2 Apply control-U operator as many times as needed to get the least significant phase bit





##3 measure the anscillar qubit in x-basis





########## Simulate the circuit and assign the output value to the variable 'x_3' 
job = execute(T_x3, sim, shots=1, memory=True)
x_3 = int(job.result().get_memory()[0])

📓Step B. Extract the middle phase bit of the $T$-gate by creating the circuit T_x2 with phase correction using x_3 value from Step A. Assign the outcome bit to the variable x_2.

Read the the second step in the previous section.

q = QuantumRegister(2)
c = ClassicalRegister(1)

T_x2 = QuantumCircuit(q,c)

########### your code goes here ##########

##1 Initialization





##2 phase correction





##3 Apply control-U operator as many times as needed 





##4 measure the anscillar qubit in x-basis





######## Simulate the circuit and assign the output value to the variable 'x_2' 
job = execute(T_x2, sim, shots=1, memory=True)
x_2 = int(job.result().get_memory()[0])

📓Step C. Find the most significant phase bit of the $T$-gate and assign it to the variable x_1.

q = QuantumRegister(2)
c = ClassicalRegister(1)

T_x1 = QuantumCircuit(q,c)

########### your code goes here #########

##1 Initialization





##2 phase correction





##3 Apply control-U operator as many times as needed to get the least significant phase bit





##4 measure the anscillar qubit in x-basis





########## Simulate the circuit and assign the output value to the variable 'x_1' 
job = execute(T_x1, sim, shots=1, memory=True)
x_1 = int(job.result().get_memory()[0])

Therefore, the $T$-gate phase bit that you found is 0.x_1x_2x_3. Run the following cell to check if your answer is correct by comparing your phase bit x_1x_2x_3 with 001, the answer in the Qiskit textbook, which corresponds to $\frac{1}{8}$ ( = 0.125), the phase of the $T$-gate.

T_phase_bits = '{}{}{}'.format(x_1, x_2, x_3) 
T_phase_bits == '001'

📓Step D. Construct the full IPE circuit and pass it to the variable qc_T ; Put the all steps that you performed into one circuit utilizing conditional operations and reset.

Instead of using three seperate circuits to extract each phase bit value, build one circuit; perform a reset operation on the auxiliary qubit after each bit gets measured into a classical register. Therefore, the circuit requires three classical registers for this example; the least significant bit measured into the classical register, c[0] and the most significant bit measured into c[2]. Implement conditional operator between the auxiliary qubit and the classical register for the phase correction.

##### your code goes here ######

















################
qc_T.draw()

Step E. Excute the following cell to perform the simulation and display the result.

count0 = execute(qc_T, sim).result().get_counts()
n=3
key_new = [str(int(key,2)/2**n) for key in list(count0.keys())]
count1 = dict(zip(key_new, count0.values()))

fig, ax = plt.subplots(1,2)
plot_histogram(count0, ax=ax[0])
plot_histogram(count1, ax=ax[1])
plt.tight_layout()

Part 2: Comparison between QPE and IPE results in the presence of noise

Goal

Understand the importance of implementing shallow circuit algorithms on current noisy quantum computers.

In Part 2 of Lab 4, we executed a Quantum Phase Estimation (QPE) circuit on a real quantum device. Having recognized the limits of the performance due to noise that presents in current quantum system, we utilized several techniques to reduce its influence on the outcome. However, the final result that was obtained, even after all these procedures, is still far from ideal. Here, we implement Iterative Phase Estimation (IPE) algorithm to overcome the effect of noise in phase estimation to a great extent and compare the result with the QPE outcome.

To investigate the impact of the noise from real quantum system on the outcome, we will perform noisy simulations of IPE circuit employing the Qiskit Aer noise module which produces a simplified noise model for an IBM quantum system. To learn more about noisy simulation, read here.

As in Lab 4, we consider to estimate the phase of p gate with $\theta = \frac{1}{3}$. Suppose that the accuracy of the estimation that we desire here is same as when the QPE circuit has four counting qubits, which determines the number of iteration and classical registers required for the IPE circuit.

📓Step A. How many classical registers is needed? Assign the value to the variable n.

## your answer goes here
n = 

📓Step B. Construct the IPE circuit in the following cell.

q = QuantumRegister(2)
c = ClassicalRegister(n)

IPE = QuantumCircuit(q,c)

########## your code goes here ############















#####################
IPE.draw()

Step C. Run the cell below to create the QPE circuit for the comparison.

def qft(n):
    """Creates an n-qubit QFT circuit"""
    circuit = QuantumCircuit(n)
    def swap_registers(circuit, n):
        for qubit in range(n//2):
            circuit.swap(qubit, n-qubit-1)
        return circuit
    def qft_rotations(circuit, n):
        """Performs qft on the first n qubits in circuit (without swaps)"""
        if n == 0:
            return circuit
        n -= 1
        circuit.h(n)
        for qubit in range(n):
            circuit.cp(np.pi/2**(n-qubit), qubit, n)
        qft_rotations(circuit, n)
    
    qft_rotations(circuit, n)
    swap_registers(circuit, n)
    return circuit


# define the parameters
t, psi = 4, 1/3*np.pi*2

# building a circuit
QPE = QuantumCircuit(t+1,t)
QPE.h(range(t))
QPE.x(t)
for idx in range(t):
    QPE.cp(psi*2**idx, idx, t)
    
qft_dag = qft(t).to_gate().inverse()
qft_dag.label = 'QFT+'
QPE.append(qft_dag, range(t))
QPE.measure(range(t), range(t))

QPE.draw()

📓Step D. Transpile the circuits for the backend ibmq_Athens.

Run the following cell to check the properties of the backend, ibmq_Athens. Pick an initial_layout, and transpile the IPE circuit setting optimization_level = 3, and save the transpiled circuit to the variable IPE_trans. Print out the depth of the transpiled circuit.

from qiskit.providers.fake_provider import FakeAthens
import qiskit.tools.jupyter

backend = FakeAthens()
backend

######## your code to transpile IPE circuit goes here ########




#####################
print(IPE_trans.depth())

Execute the cell below to transpile QPE circuit.

num = 500
QPE_trans = transpile([QPE]*num, backend, optimization_level=3)
QPE_trans_depth = np.array([QPE_trans[idx].depth() for idx in range(num)])
print(min(QPE_trans_depth), max(QPE_trans_depth))
best_arg = np.argmin(QPE_trans_depth)
QPE_trans_best = QPE_trans[best_arg]

Step E. Run the following cells to perform the noise simulation of the transipiled circuits.

from qiskit.providers.aer.noise import NoiseModel

noise_model = NoiseModel.from_backend(backend)
shots = 20000

counts = execute([IPE_trans, QPE_trans_best], sim, noise_model=noise_model).result().get_counts()

Step F. Execute the cell below to compute the exact phase estimation results.

from qiskit.quantum_info import Statevector

QPE_exact = QuantumCircuit(t+1)
QPE_exact.h(range(t))
QPE_exact.x(t)
for idx in range(t):
    QPE_exact.cp(psi*2**idx, idx, t)
    
qft_dag = qft(t).to_gate().inverse()
qft_dag.label = 'QFT+'
QPE_exact.append(qft_dag, range(t))

#QPE_exact.draw('mpl')

state = Statevector.from_instruction(QPE_exact)
pmf = state.probabilities_dict(range(4))

Step G. Show the comparison figure by running the following cell.

def count_new(count):
    phi_est = np.array([round(int(key, 2)/2**t, 3) for key in list(count.keys())])
    key_new = list(map(str, phi_est))
    count_new = dict(zip(key_new, count.values()))
    return count_new

pmf_new = count_new(pmf)
count_IPE = count_new(counts[0])
count_QPE = count_new(counts[1])

fig, ax = plt.subplots(1, 2, figsize=(32,10))
fig.suptitle('QPE .vs. IPE for estimating $\\theta=1/3$', fontsize=23)
plot_histogram([pmf_new, count_QPE], ax=ax[0], legend=['No_noise', 'Athens_sim'])
plot_histogram([pmf_new, count_IPE], ax=ax[1], legend=['No_noise', 'Athens_sim'])
ax[0].set_title('QPE circuit result', fontsize=16)
ax[0].set_xlabel('$\\theta$', fontsize=16)
ax[1].set_title('IPE circuit result', fontsize=16)
ax[1].set_xlabel('$\\theta$', fontsize=16)
plt.show()

If you create the IPE circuit successfully to estimate the phase, $\theta = \frac{1}{3}$, you would get the similar plots as shown below.

📓Step G. Discuss about the results.