English
Languages
English
Japanese
Spanish
Table of Contents

Note

This page was generated from docs/tutorials/qpe_with_sampler.ipynb.

Quantum phase estimation using the Sampler primitive

The quantum phase estimation (QPE) algorithm is an important subroutine in quantum computation. It serves as a central building block of many quantum algorithms including the Shor’s factoring algorithm and the quantum algorithm for linear systems of equations (HHL algorithm). In this tutorial, you will build a parameterized version of QPE circuit and run it using the Sampler primitive, which makes running parameterized circuits very easy.

Before you begin

This tutorial requires a Qiskit Runtime service instance. If you haven’t done so already, follow the instructions in the Qiskit Getting started guide to set one up.

Background information

Quantum phase estimation algorithm

The QPE algorithm [1][2] estimates the value of theta, where a unitary operator \(U\) has an eigenvector \(|\psi \rangle\) with eigenvalue \(e^{2\pi i \theta}\). To find the estimation, we assume we have black boxes (sometimes called oracles), that can prepare the state \(|\psi \rangle\) and run a controlled-\(U^{2^j}\) operation.

Because it relies on black boxes, the QPE algorithm is not actually a complete algorithm, but can be thought of as a subroutine. It can work with other subroutines to perform other computations, such as Shor’s factoring algorithm and the HHL algorithm.

We are not going to explain the details of the QPE algorithm here. If you want to learn more, you can read the chapter about the QPE algorithm in the Qiskit textbook [3].

Overview

As explained earlier, there are black boxes in the QPE algorithm to prepare the state \(|\psi \rangle\) and perform the controlled-\(U^{2^j}\) operation. In this tutorial, you will prepare a series of QPE circuits containing different black boxes corresponding to different phases. You will then run these circuits using the Sampler primitive and analyze the results. As you will see, the Sampler primitive makes running parameterized circuits very easy. Rather than create a series of QPE circuits, you only need to create one QPE circuit with a parameter specifying the phase the black boxes generate and a series of phase values for the parameter.

Note

In a typical usage of the QPE algorithm, such as the Shor’s factoring algorithm, the black boxes are not configured by the user but instead are specified by other subroutines. However, the purpose of this tutorial is to use the QPE algorithm to illustrate the ease of running parameterized circuits by using the Sampler primitive.

Step 1: Create QPE circuits

Create a parameterized QPE circuit

The procedure of the QPE algorithm is as follows:

  1. Create a circuit with two quantum registers (the first for estimating the phase and the second for storing the eigenvector \(|\psi\rangle\)) and one classical register for readout.

  2. Initialize the qubits in the first register as \(|0\rangle\) and the second register as \(|\psi\rangle\).

  3. Create a superposition in the first register.

  4. Apply the controlled-\(U^{2^j}\) black box.

  5. Apply an inverse quantum Fourier transform (QFT) to the first register.

  6. Measure the first register.

The following code creates a function create_qpe_circuit for creating a QPE circuit given the keyword arguments theta and num_qubits. Note that theta doesn’t contain the \(2\pi\) factor. The num_qubits argument specifies the number of qubits in the first register. More qubits will provide more precision for the phase estimation.

[1]:

import numpy as np

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from qiskit.circuit.library import QFT


def create_qpe_circuit(theta, num_qubits):
    """Creates a QPE circuit given theta and num_qubits."""

    # Step 1: Create a circuit with two quantum registers and one classical register.
    first = QuantumRegister(
        size=num_qubits, name="first"
    )  # the first register for phase estimation
    second = QuantumRegister(
        size=1, name="second"
    )  # the second register for storing eigenvector |psi>
    classical = ClassicalRegister(
        size=num_qubits, name="readout"
    )  # classical register for readout
    qpe_circuit = QuantumCircuit(first, second, classical)

    # Step 2: Initialize the qubits.
    # All qubits are initialized in |0> by default, no extra code is needed to initialize the first register.
    qpe_circuit.x(
        second
    )  # Initialize the second register with state |psi>, which is |1> in this example.

    # Step 3: Create superposition in the first register.
    qpe_circuit.barrier()  # Add barriers to separate each step of the algorithm for better visualization.
    qpe_circuit.h(first)

    # Step 4: Apply a controlled-U^(2^j) black box.
    qpe_circuit.barrier()
    for j in range(num_qubits):
        qpe_circuit.cp(
            theta * 2 * np.pi * (2**j), j, num_qubits
        )  # Theta doesn't contain the 2 pi factor.

    # Step 5: Apply an inverse QFT to the first register.
    qpe_circuit.barrier()
    qpe_circuit.compose(QFT(num_qubits, inverse=True), inplace=True)

    # Step 6: Measure the first register.
    qpe_circuit.barrier()
    qpe_circuit.measure(first, classical)

    return qpe_circuit

If you pass a real number to the theta argument in the create_qpe_circuit function, it will return a circuit with a fixed phase encoded.

[2]:

num_qubits = 4
qpe_circuit_fixed_phase = create_qpe_circuit(
    1 / 2, num_qubits
)  # Create a QPE circuit with fixed theta=1/2.
qpe_circuit_fixed_phase.draw("mpl")

[2]:
../_images/tutorials_qpe_with_sampler_9_0.png

However, if you pass a `Parameter <https://qiskit.org/documentation/stubs/qiskit.circuit.Parameter.html>`__ object instead, it will return a parameterized circuit whose parameter values can be assigned after the circuit has been created. You will use the parameterized version of the QPE circuit for the rest of the tutorial.

[3]:

from qiskit.circuit import Parameter

theta = Parameter(
    "theta"
)  # Create a parameter `theta` whose values can be assigned later.
qpe_circuit_parameterized = create_qpe_circuit(theta, num_qubits)
qpe_circuit_parameterized.draw("mpl")

[3]:
../_images/tutorials_qpe_with_sampler_11_0.png

Create a list of phase values to be assigned later

After creating the parameterized QPE circuit, you will create a list of phase values to be assigned to the circuit in the next step. You can use the following code to create a list of 21 phase values that range from 0 to 2 with equal spacing, that is, 0, 0.1, 0.2, …, 1.9, 2.0.

[4]:

number_of_phases = 21
phases = np.linspace(0, 2, number_of_phases)
individual_phases = [
    [ph] for ph in phases
]  # Phases need to be expressed as a list of lists.

Step 2: Submit the circuits to a quantum computer on the cloud

Connect to the Qiskit Runtime service

First, you will connect to the Qiskit Runtime service instance that you created in the first step. We will use ibmq_qasm_simulator to run this tutorial.

[5]:

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = "ibmq_qasm_simulator"  # use the simulator

Run the parameterized circuits by using the Sampler primitive

With the Qiskit Runtime service connected, a parameterized QPE circuit, and a list of parameter values, you are now ready to use the Sampler primitive to run the QPE circuits. The with ... syntax opens a Qiskit Runtime session. Within the session, you will run the parameterized QPE circuit with 21 parameter values with just one line of code. The Sampler primitive will take the parameter values, assign them to the circuit and run it as 21 circuits with fixed parameter values and return a single job result containing the (quasi-)probabilities of bit strings. This saves you from writing dozens of lines of code.

The Estimator primitive has a similar API interface for parameterized circuits. See the API reference for more details.

[6]:

from qiskit_ibm_runtime import Sampler, Session

with Session(service=service, backend=backend):
    results = (
        Sampler()
        .run(
            [qpe_circuit_parameterized] * len(individual_phases),
            parameter_values=individual_phases,
        )
        .result()
    )

Step 3: Analyze the results

Analyze the results of one circuit

After the job has been completed, you can start analyzing the results by looking at the histogram of one of the circuits. The (quasi-)probability distribution of the measurement outcome is stored in results.quasi_dists and you can access results from an individual circuit by passing the index of the circuit (\(idx=6\) in the following example) you are interested in.

[7]:

from qiskit.tools.visualization import plot_histogram

idx = 6
plot_histogram(
    results.quasi_dists[idx].binary_probabilities(),
    legend=[f"$\\theta$={phases[idx]:.3f}"],
)

[7]:
../_images/tutorials_qpe_with_sampler_21_0.png

To estimate \(\theta\), you need to divide the most likely outcome \(N_1\) (1010 in binary or 10 in decimal) by \(2^n\), where \(n\) is the number of qubits in the first register (\(n=4\) in our example):

\[\theta = \frac{N_1}{2^n} = \frac{10}{2^4} = 0.625.\]

This is close to the expected value of \(0.6\), but you can get a better estimate by taking the weighted average of the most likely outcome (1010 in binary or 10 in decimal) and the second most likely outcome (1001 in binary or 9 in decimal):

\[\theta = \frac{P_1 \times \frac{N_1}{2^n} + P_2 \times \frac{N_2}{2^n}}{P_1 + P_2} = \frac{0.554 \times \frac{10}{2^4} + 0.269 \times \frac{9}{2^4}}{0.554 + 0.269} = 0.605,\]

where \(N_1\) and \(P_1\) are the most likely outcome and its probability and \(N_2\) and \(P_2\) are the second most likely outcome and its probability. The result, \(0.605\), is much closer to the expected value of \(0.6\). We will use this method to analyze the results from all circuits.

Analyze the results of all circuits

You can use the following helper functions to find the values for \(N_1\), \(P_1\), \(N_2\), and \(P_2\) for a phase and calculate the estimated \(\theta\). Ideally \(N_2\) should be a “neighbor” of \(N_1\) (for example, the neighbors of 1010 are 1001 and 1011). However, due to the presence of noise, the second most likely outcome may not be a neighbor of \(N_1\) if the results were obtained from a real quantum system. The helper functions take the \(N_2\) value only from \(N_1\)’s neighbors.

[8]:

def most_likely_bitstring(results_dict):
    """Finds the most likely outcome bit string from a result dictionary."""
    return max(results_dict, key=results_dict.get)


def find_neighbors(bitstring):
    """Finds the neighbors of a bit string.

    Example:
        For bit string '1010', this function returns ('1001', '1011')
    """
    if bitstring == len(bitstring) * "0":
        neighbor_left = len(bitstring) * "1"
    else:
        neighbor_left = format((int(bitstring, 2) - 1), "0%sb" % len(bitstring))

    if bitstring == len(bitstring) * "1":
        neighbor_right = len(bitstring) * "0"
    else:
        neighbor_right = format((int(bitstring, 2) + 1), "0%sb" % len(bitstring))

    return (neighbor_left, neighbor_right)


def estimate_phase(results_dict):
    """Estimates the phase from a result dictionary of a QPE circuit."""

    # Find the most likely outcome bit string N1 and its neighbors.
    num_1_key = most_likely_bitstring(results_dict)
    neighbor_left, neighbor_right = find_neighbors(num_1_key)

    # Get probabilities of N1 and its neighbors.
    num_1_prob = results_dict.get(num_1_key)
    neighbor_left_prob = results_dict.get(neighbor_left)
    neighbor_right_prob = results_dict.get(neighbor_right)

    # Find the second most likely outcome N2 and its probability P2 among the neighbors.
    if neighbor_left_prob is None:
        # neighbor_left doesn't exist
        if neighbor_right_prob is None:
            # both neighbors don't exist, N2 is N1
            num_2_key = num_1_key
            num_2_prob = num_1_prob
        else:
            # If only neighbor_left doesn't exist, N2 is neighbor_right.
            num_2_key = neighbor_right
            num_2_prob = neighbor_right_prob
    elif neighbor_right_prob is None:
        # If only neighbor_right doesn't exist, N2 is neighbor_left.
        num_2_key = neighbor_left
        num_2_prob = neighbor_left_prob
    elif neighbor_left_prob > neighbor_right_prob:
        # Both neighbors exist and neighbor_left has higher probability, so N2 is neighbor_left.
        num_2_key = neighbor_left
        num_2_prob = neighbor_left_prob
    else:
        # Both neighbors exist and neighbor_right has higher probability, so N2 is neighor_right.
        num_2_key = neighbor_right
        num_2_prob = neighbor_right_prob

    # Calculate the estimated phases for N1 and N2.
    num_qubits = len(num_1_key)
    num_1_phase = int(num_1_key, 2) / 2**num_qubits
    num_2_phase = int(num_2_key, 2) / 2**num_qubits

    # Calculate the weighted average phase from N1 and N2.
    phase_estimated = (num_1_phase * num_1_prob + num_2_phase * num_2_prob) / (
        num_1_prob + num_2_prob
    )

    return phase_estimated

You can use the estimate_phase helper function to estimate phases from the results of all circuits.

[9]:

qpe_solutions = []
for idx, result_dict in enumerate(results.quasi_dists):
    qpe_solutions.append(estimate_phase(result_dict.binary_probabilities()))

The ideal solutions for the phases \(\theta\) are periodic with a period of 1 because the eigenvalue \(e^{2 \pi i \theta}\) is \(2 \pi\) periodic. You can run the following cell to generate the ideal solutions for comparison with the solutions obtained from the QPE algorithm.

[10]:

ideal_solutions = np.append(
    phases[: (number_of_phases - 1) // 2],  # first period
    np.subtract(phases[(number_of_phases - 1) // 2 : -1], 1),  # second period
)
ideal_solutions = np.append(
    ideal_solutions, np.subtract(phases[-1], 2)
)  # starting point of the third period

You can run the following code to visualize the solutions.

[11]:

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(10, 6))
plt.plot(phases, ideal_solutions, "--", label="Ideal solutions")
plt.plot(phases, qpe_solutions, "o", label="QPE solutions")

plt.title("Quantum Phase Estimation Algorithm")
plt.xlabel("Input Phase")
plt.ylabel("Output Phase")
plt.legend()

[11]:

<matplotlib.legend.Legend at 0x130347bb0>
../_images/tutorials_qpe_with_sampler_30_1.png

As you can see, the solutions obtained from the QPE algorithm follow closely to the ideal solutions. Congratulations! You have successfully run the QPE algorithm and obtained good solutions!

Summary

In this tutorial, you have created a parameterized QPE circuit, run it using the Sampler primitive, analyzed the results, and obtained solutions that are closed to the ideal solutions. You can see that the Sampler primitive makes running parameterized circuits very easy.

References

  1. Kitaev, A. Y. (1995). Quantum measurements and the Abelian stabilizer problem. arXiv preprint quant-ph/9511026.

  2. Nielsen, M., & Chuang, I. (2010). Quantum computation and quantum information (2nd ed., pp. 221-226). Cambridge University Press.

  3. Abbas, A., Andersson, S., Asfaw, A., Córcoles, A., Bello, L., Ben-Haim, Y., … & Wootton, J. (2020). Learn quantum computation using qiskit. URL: https://qiskit.org/textbook/preface.html (accessed 07/14/2022).