QUANTUM PHASE ESTIMATION

This tutorial provides a detailed implementation of the Quantum Phase Estimation (QPE) algorithm using the Amazon Braket SDK. The QPE algorithm is designed to estimate the eigenvalues of a unitary operator [1, 2]; it is a very important subroutine to many quantum algorithms, most famously Shor's algorithm for factoring and the HHL algorithm (named after the physicists Harrow, Hassidim and Lloyd) for solving linear systems of equations on a quantum computer [1, 2]. Moreover, eigenvalue problems can be found across many disciplines and application areas, including (for example) principal component analysis (PCA) as used in machine learning or the solution of differential equations as relevant across mathematics, physics, engineering and chemistry. We first review the basics of the QPE algorithm. We then implement the QPE algorithm in code using the Amazon Braket SDK, and we illustrate the application thereof with simple examples. This notebook also showcases the Amazon Braket circuit.subroutine functionality, which allows us to use custom-built gates as if they were any other built-in gates. This tutorial is set up to run either on the local simulator or the on-demand simulators; changing between these devices merely requires changing one line of code as demonstrated as follows in cell [4].

TECHNICAL BACKGROUND OF QPE

Introduction: A unitary matrix is a complex, square matrix whose adjoint (or conjugate transpose) is equal to its inverse. Unitary matrices have many nice properties, including the fact that their eigenvalues are always roots of unity (that is, phases). Given a unitary matrix (satisfying ) and an eigenstate with , the Quantum Phase Estimation (QPE) algorithm provides an estimate for the phase (with since the eigenvalues of a unitary have modulus one). The QPE works with high probability within an additive error using qubits (without counting the qubits used to encode the eigenstate) and controlled- operations [1].

Quantum Phase Estimation Algorithm: The QPE algorithm takes a unitary as input. For the sake of simplicity (we will generalize the discussion below), suppose that the algorithm also takes as input an eigenstate fulfilling

with .

QPE uses two registers of qubits: we refer to the first register as precision qubits (as the number of qubits in the first register sets the achievable precision of our results) and the second register as query qubits (as the second register hosts the eigenstate ). Suppose we have prepared this second register in . We then prepare a uniform superposition of all basis vectors in the first register using a series of Hadamard gates.

Next, we apply a series of controlled-unitaries for different powers of (as illustrated in the circuit diagram that follows). For example, for we get

Note that the second register remains unaffected as it stays in the eigenstate . However, we managed to transfer information about the phase of the eigenvalue of (that is, ) into the first precision register by encoding it as a relative phase in the state of the qubits in the first register.

Similarly, for we obtain

where this time we wrote into the precision register. The process is similar for all .

Introducing the following notation for binary fractions

one can show that the application of a controlled unitary leads to the following transformation

where the first bits of precision in the binary expansion (that is, those bits to the left of the decimal) can be dropped, because for any whole number .

The QPE algorithm implements a series of these transformations for , using qubits in the precision register. In its entirety, this sequence of controlled unitaries leads to the transformation

By inspection, one can see that the state of the register qubits above corresponds to a quantum Fourier transform of the state . Thus, the final step of the QPE algorithm is to run the inverse Quantum Fourier Transform (QFT) algorithm on the precision register to extract the phase information from this state. The resulting state is

Measuring the precision qubits in the computational basis then gives the classical bitstring , from which we can readily infer the phase estimate with the corresponding eigenvalue .

Simple example for illustration: For concreteness, consider a simple example with the unitary given by the Pauli gate, , for which is an eigenstate with eigenvalue , i.e., . This state can be prepared with a Hadamard gate as . We take a precision register consisting of just two qubits ().

Thus, after the first layer of Hadamard gates, the quantum state is

Next, the applications of the controlled- gates (equal to operations, or CNOT gates in this example) leave this state untouched, because is an eigenstate of with eigenvalue . Finally, applying the inverse QFT leads to

from which we deduce and therefore , as expected. Here, in the last step we have used , which makes the effect of the inverse QFT more apparent.

Initial state of query register: So far, we have assumed that the query register is prepared in an eigenstate of . What happens if this is not the case? Let's reconsider the simple example given previously.

Suppose now that the query register is instead prepared in the state . We can always express this state in the eigenbasis of , that is, . By linearity, application of the QPE algorithm then gives (up to normalization)

When we measure the precision qubits in this state, 50% of the time we will observe the eigenphase and 50% of the time we will measure . We illustrate this example numerically as follows.

This example motivates the general case: we can pass a state that is not an eigenstate of to the QPE algorithm, but we may need to repeat our measurements several times in order to obtain an estimate of the desired phase.

CIRCUIT IMPLEMENTATION OF QPE

The QPE circuit can be implemented using Hadamard gates, controlled- unitaries, and the inverse QFT (denoted as ). The details of the calculation can be found in a number of resources (such as, [1]); we omit them here. Following the previous discussion, the circuit that implements the QPE algorithm reads as below, where m is the size of lower query register and n is the size of upper precision register.

IMPORTS and SETUP

Pauli Matrices:

In some of our examples, we choose the unitary to be given by the Pauli Matrices, which we thus define as follows:

IMPLEMENTATION OF THE QPE CIRCUIT

In utils_qpe.py we provide simple helper functions to implement the quantum circuit for the QPE algorithm. Specifically, we demonstrate that such modular building blocks can be registered as subroutines, using @circuit.subroutine(register=True). Moreover, we provide a helper function (called get_qpe_phases) to perform postprocessing based on the measurement results to extract the phase. The details of utils_qpe.py are shown in the Appendix.

To implement the unitary , one can use the fact that , so that can be constructed by repeatedly applying the core building block . However, the circuit generated using this approach will have a significantly larger depth. In our implementation, we instead define the matrix and create the controlled gate from that.

VISUALIZATION OF THE QFT CIRCUIT

To check our implementation of the QPE circuit, we visualize this circuit for a small number of qubits.

As shown in the following code, the two registers can be distributed anywhere across the circuit, with arbitrary indices for the precision and the query registers.

As follows, we set up the same circuit, this time implementing the unitary , by repeatedly applying the core building block . This operation can be done by setting the parameter control_unitary=False (default is True).

In the circuit diagram, we can visually infer the exponents for , at the expense of a larger circuit depth.

NUMERICAL TEST EXPERIMENTS

In the following section, we verify that our QFT implementation works as expected with a few test examples:

1. We run QPE with and prepare the eigenstate with phase and eigenvalue .
2. We run QPE with and prepare the eigenstate with phase and eigenvalue .
3. We run QPE with and prepare which is not an eigenstate of . Because , we expect to measure both and associated with the two eigenstates .
4. We run QPE with unitary , and prepare the query register in the eigenstate . Here, we expect to measure the phase (giving the corresponding eigenvalue ).
5. We run QPE with a random two qubit unitary, diagonal in the computational basis, and prepare the query register in the eigenstate . In this case, we should be able to read off the eigenvalue and phase from and verify QPE gives the right answer (with high probability) up to a small error (that depends on the number of qubits in the precision register).

HELPER FUNCTIONS FOR NUMERICAL TESTS

Because we will run the same code repeatedly, let's first create a helper function we can use to keep the notebook clean.

NUMERICAL TEST EXAMPLE 1

First, apply the QPE algorithm to the simple single-qubit unitary , with eigenstate . Here, we expect to measure the phase (giving the corresponding eigenvalue ). We show that this result stays the same as we increase the number of qubits for the top register.

Next, check that we get the same result for a larger precision (top) register.

NUMERICAL TEST EXAMPLE 2

Next, apply the QPE algorithm to the simple single-qubit unitary , with eigenstate . Here, we expect to measure the phase (giving the corresponding eigenvalue ).

NUMERICAL TEST EXAMPLE 3

Next, apply the QPE algorithm again to the simple single-qubit unitary , but we initialize the query register in the state which is not an eigenstate of . Here, following the previous discussion, we expect to measure the phases (giving the corresponding eigenvalue ). Accordingly, here we set items_to_keep=2.

NUMERICAL TEST EXAMPLE 4

Next, apply the QPE algorithm to the two-qubit unitary , and prepare the query register in the eigenstate . Here, we expect to measure the phase (giving the corresponding eigenvalue ).

NUMERICAL TEST EXAMPLE 5

In this example, we choose the unitary to be a random two-qubit unitary, diagonal in the computational basis. We initialize the query register to be in the eigenstate of , which we can prepare using that . In this case we should be able to read off the eigenvalue and phase from and verify that QPE gives the right answer.

When we execute the QPE circuit, we expect the following (approximate) result for the eigenvalue estimate:

We can easily improve the precision of our parameter estimate by increasing the number of qubits in the precision register, as shown in the following example.

---

APPENDIX

Details of the utiles_qpe.py module

Imports, including inverse QFT

# general imports
import numpy as np
import math
from collections import Counter
from datetime import datetime
import pickle

# AWS imports: Import Braket SDK modules
from braket.circuits import Circuit, circuit

# local imports
from utils_qft import inverse_qft

QPE Subroutine

@circuit.subroutine(register=True)
def controlled_unitary(control, target_qubits, unitary):
    """
    Construct a circuit object corresponding to the controlled unitary

    Args:
        control: The qubit on which to control the gate

        target_qubits: List of qubits on which the unitary U acts

        unitary: matrix representation of the unitary we wish to implement in a controlled way
    """

    # Define projectors onto the computational basis
    p0 = np.array([[1., 0.],
                   [0., 0.]])

    p1 = np.array([[0., 0.],
                   [0., 1.]])

    # Instantiate circuit object
    circ = Circuit()

    # Construct numpy matrix
    id_matrix = np.eye(len(unitary))
    controlled_matrix = np.kron(p0, id_matrix) + np.kron(p1, unitary)

    # Set all target qubits
    targets = [control] + target_qubits

    # Add controlled unitary
    circ.unitary(matrix=controlled_matrix, targets=targets)

    return circ


@circuit.subroutine(register=True)
def qpe(precision_qubits, query_qubits, unitary, control_unitary=True):
    """
    Function to implement the QPE algorithm using two registers for precision (read-out) and query.
    Register qubits need not be contiguous.

    Args:
        precision_qubits: list of qubits defining the precision register

        query_qubits: list of qubits defining the query register

        unitary: Matrix representation of the unitary whose eigenvalues we wish to estimate

        control_unitary: Optional boolean flag for controlled unitaries,
                         with C-(U^{2^k}) by default (default is True),
                         or C-U controlled-unitary (2**power) times
    """
    qpe_circ = Circuit()

    # Get number of qubits
    num_precision_qubits = len(precision_qubits)
    num_query_qubits = len(query_qubits)

    # Apply Hadamard across precision register
    qpe_circ.h(precision_qubits)

    # Apply controlled unitaries. Start with the last precision_qubit, and end with the first
    for ii, qubit in enumerate(reversed(precision_qubits)):
        # Set power exponent for unitary
        power = ii

        # Alterantive 1: Implement C-(U^{2^k})
        if control_unitary:
            # Define the matrix U^{2^k}
            Uexp = np.linalg.matrix_power(unitary,2**power)

            # Apply the controlled unitary C-(U^{2^k})
            qpe_circ.controlled_unitary(qubit, query_qubits, Uexp)
        # Alterantive 2: One can instead apply controlled-unitary (2**power) times to get C-U^{2^power}
        else:
            for _ in range(2**power):
                qpe_circ.controlled_unitary(qubit, query_qubits, unitary)

    # Apply inverse qft to the precision_qubits
    qpe_circ.inverse_qft(precision_qubits)

    return qpe_circ

QPE postprocessing helper functions

# helper function to remove query bits from bitstrings
def substring(key, precision_qubits):
    """
    Helper function to get substring from keys for dedicated string positions as given by precision_qubits.
    This function is necessary to allow for arbitrary qubit mappings in the precision and query registers
    (that is, so that the register qubits need not be contiguous.)

    Args:
        key: string from which we want to extract the substring supported only on the precision qubits

        precision_qubits: List of qubits corresponding to precision_qubits.
                          Currently assumed to be a list of integers corresponding to the indices of the qubits
    """
    short_key = ''
    for idx in precision_qubits:
        short_key = short_key + key[idx]

    return short_key


# helper function to convert binary fractional to decimal
# reference: https://www.geeksforgeeks.org/convert-binary-fraction-decimal/
def binaryToDecimal(binary):
    """
    Helper function to convert binary string (example: '01001') to decimal

    Args:
        binary: string which to convert to decimal fraction
    """

    length = len(binary)
    fracDecimal = 0

    # Convert fractional part of binary to decimal equivalent
    twos = 2

    for ii in range(length):
        fracDecimal += ((ord(binary[ii]) - ord('0')) / twos);
        twos *= 2.0

    # return fractional part
    return fracDecimal


# helper function for postprocessing based on measurement shots
def get_qpe_phases(measurement_counts, precision_qubits, items_to_keep=1):
    """
    Get QPE phase estimate from measurement_counts for given number of precision qubits

    Args:
        measurement_counts: measurement results from a device run

        precision_qubits: List of qubits corresponding to precision_qubits.
                          Currently assumed to be a list of integers corresponding to the indices of the qubits

        items_to_keep: number of items to return (topmost measurement counts for precision register)
    """

    # Aggregate the results (that is, ignore the query register qubits):

    # First get bitstrings with corresponding counts for precision qubits only
    bitstrings_precision_register = [substring(key, precision_qubits)  for key in measurement_counts.keys()]
    # Then keep only the unique strings
    bitstrings_precision_register_set = set(bitstrings_precision_register)
    # Cast as a list for later use
    bitstrings_precision_register_list = list(bitstrings_precision_register_set)

    # Now create a new dict to collect measurement results on the precision_qubits.
    # Keys are given by the measurement count substrings on the register qubits. Initialize the counts to zero.
    precision_results_dic = {key: 0 for key in bitstrings_precision_register_list}

    # Loop over all measurement outcomes
    for key in measurement_counts.keys():
        # Save the measurement count for this outcome
        counts = measurement_counts[key]
        # Generate the corresponding shortened key (supported only on the precision_qubits register)
        count_key = substring(key, precision_qubits)
        # Add these measurement counts to the corresponding key in our new dict
        precision_results_dic[count_key] += counts

    # Get topmost values only
    c = Counter(precision_results_dic)
    topmost= c.most_common(items_to_keep)
    # get decimal phases from bitstrings for topmost bitstrings
    phases_decimal = [binaryToDecimal(item[0]) for item in topmost]

    # Get decimal phases from bitstrings for all bitstrings
    # number_precision_qubits = len(precision_qubits)
    # Generate binary decimal expansion
    # phases_decimal = [int(key, 2)/(2**number_precision_qubits) for key in precision_results_dic]
    # phases_decimal = [binaryToDecimal(key) for key in precision_results_dic]

    return phases_decimal, precision_results_dic

Run QPE experiments:

def run_qpe(unitary, precision_qubits, query_qubits, query_circuit,
            device, items_to_keep=1, shots=1000, save_to_pck=False):
    """
    Function to run QPE algorithm end-to-end and return measurement counts.

    Args:
        precision_qubits: list of qubits defining the precision register

        query_qubits: list of qubits defining the query register

        unitary: Matrix representation of the unitary whose eigenvalues we wish to estimate

        query_circuit: query circuit for state preparation of query register

        items_to_keep: (optional) number of items to return (topmost measurement counts for precision register)

        device: Braket device backend

        shots: (optional) number of measurement shots (default is 1000)

        save_to_pck: (optional) save results to pickle file if True (default is False)
    """

    # get size of precision register and total number of qubits
    number_precision_qubits = len(precision_qubits)
    num_qubits = len(precision_qubits) + len(query_qubits)

    # Define the circuit. Start by copying the query_circuit, then add the QPE:
    circ = query_circuit
    circ.qpe(precision_qubits, query_qubits, unitary)

    # Add desired results_types
    circ.probability()

    # Run the circuit with all zeros input.
    # The query_circuit subcircuit generates the desired input from all zeros.
    # The following code executes the correct device.run call, depending on whether the backend is local or on-demand
    task = device.run(circ, shots=shots)

    # get result for this task
    result = task.result()

    # get metadata
    metadata = result.task_metadata

    # get output probabilities (see result_types above)
    probs_values = result.values[0]

    # get measurement results
    measurements = result.measurements
    measured_qubits = result.measured_qubits
    measurement_counts = result.measurement_counts
    measurement_probabilities = result.measurement_probabilities

    # bitstrings
    format_bitstring = '{0:0' + str(num_qubits) + 'b}'
    bitstring_keys = [format_bitstring.format(ii) for ii in range(2**num_qubits)]

    # QPE postprocessing
    phases_decimal, precision_results_dic = get_qpe_phases(measurement_counts, precision_qubits, items_to_keep)
    eigenvalues = [np.exp(2*np.pi*1j*phase) for phase in phases_decimal]

    # aggregate results
    out = {'circuit': circ,
           'task_metadata': metadata,
           'measurements': measurements,
           'measured_qubits': measured_qubits,
           'measurement_counts': measurement_counts,
           'measurement_probabilities': measurement_probabilities,
           'probs_values': probs_values,
           'bitstring_keys': bitstring_keys,
           'precision_results_dic': precision_results_dic,
           'phases_decimal': phases_decimal,
           'eigenvalues': eigenvalues}

    if save_to_pck:
        # store results: dump output to pickle with timestamp in filename
        time_now = datetime.strftime(datetime.now(), '%Y%m%d%H%M%S')
        results_file = 'results-'+time_now+'.pck'
        pickle.dump(out, open(results_file, "wb"))
        # you can load results as follows
        # out = pickle.load(open(results_file, "rb"))

    return out

---

REFERENCES

[1] Wikipedia: https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm

[2] Nielsen, Michael A., Chuang, Isaac L. (2010). Quantum Computation and Quantum Information (2nd ed.). Cambridge: Cambridge University Press.