From 0950594c8d3e6aed637a47255a8a429b64c69874 Mon Sep 17 00:00:00 2001 From: John Lapeyre Date: Thu, 25 Mar 2021 06:44:33 -0700 Subject: [PATCH] Implement quantum phase estimation algorithms (#5642) * Add PhaseEstimator and related classes * Fix linter complaints * Update qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation.py Co-authored-by: Julien Gacon * Fix indentation in doc string * Fix more indentation in doc string * Comment on reasons for removing identity term from Hamiltonian * Replace mistaken description "water" by "H2" * Add default phase shift of 0.0 for identity term * Remove function call from default arg Replace default arg with None, and then set the desired default in the body of the function. * Refactor _run and rename to _result Move all common code out of _run for HamiltonianPhaseEstimation and PhaseEstimation, and rename _run to _result, because it only creates the result class. * Move some inputs and calculation to estimate method * Add cutoff to filter in test_from_bound * An unknown change in the past three weeks caused tests to fail with the default filter threshold of 0 (no filtering). I added a small cutoff so that the tests pass again. * Fix linter complaints * Fix linter complaints * Allow PauliSumOp as input Hamiltonian * Make suggested changes * Move module method to classmethod * Fix indentation * Refactor scale_phase scale_phases * Make HamiltonianPhaseEstimation have a hasa rather than isa use of PE * Remove all computation from PhaseEstimator * Remove ability to set number of computation qubits in PhaseEstimator.estimate() * update to current algorithm paradigm + small fixes * make algorithms stateless (store no problem information) * fix lint and docstrings * fix cyclic imports * Use PauliTrotterEvolution by default in HamiltonianPhaseEstimation * Make HamiltonianPhaseEstmation take a StateFn for state preparation * Add method most_likely_phase to HamiltonianPhaseEstimation * Allow PauliOp as input unitary to HPE * Allow MatrixOp as input Hermitian op to HPE * Fix linter complaints * fix evolution = None case * refactor tests slightly * be explicit about which evolution is used * combine tests using ddt * add some linebreaks * fix import order * Improve docstrings for PE and HPE * attempt to fix sphinx * attempt to fix sphinx #2 * attempt no. 3 to fix sphinx * layout, don't use function calls in defaults * trailing comma Co-authored-by: Julien Gacon --- qiskit/algorithms/__init__.py | 21 ++ .../algorithms/phase_estimators/__init__.py | 29 ++ .../hamiltonian_phase_estimation.py | 211 +++++++++++++++ .../hamiltonian_phase_estimation_result.py | 104 +++++++ .../phase_estimators/phase_estimation.py | 225 ++++++++++++++++ .../phase_estimation_result.py | 154 +++++++++++ .../phase_estimation_scale.py | 141 ++++++++++ .../phase_estimators/phase_estimator.py | 50 ++++ .../python/algorithms/test_phase_estimator.py | 254 ++++++++++++++++++ 9 files changed, 1189 insertions(+) create mode 100644 qiskit/algorithms/phase_estimators/__init__.py create mode 100644 qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation.py create mode 100644 qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation_result.py create mode 100644 qiskit/algorithms/phase_estimators/phase_estimation.py create mode 100644 qiskit/algorithms/phase_estimators/phase_estimation_result.py create mode 100644 qiskit/algorithms/phase_estimators/phase_estimation_scale.py create mode 100644 qiskit/algorithms/phase_estimators/phase_estimator.py create mode 100644 test/python/algorithms/test_phase_estimator.py diff --git a/qiskit/algorithms/__init__.py b/qiskit/algorithms/__init__.py index 68f4d0604c9b..d6333da36758 100644 --- a/qiskit/algorithms/__init__.py +++ b/qiskit/algorithms/__init__.py @@ -116,6 +116,20 @@ QAOA VQE +Phase Estimators +++++++++++++++++ +Algorithms that estimate the phases of eigenstates of a unitary. + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + + HamiltonianPhaseEstimation + HamiltonianPhaseEstimationResult + PhaseEstimationScale + PhaseEstimation + PhaseEstimationResult + Exceptions ========== @@ -141,6 +155,8 @@ from .minimum_eigen_solvers import (VQE, VQEResult, QAOA, NumPyMinimumEigensolver, MinimumEigensolver, MinimumEigensolverResult) +from .phase_estimators import (HamiltonianPhaseEstimation, HamiltonianPhaseEstimationResult, + PhaseEstimationScale, PhaseEstimation, PhaseEstimationResult) from .exceptions import AlgorithmError __all__ = [ @@ -172,5 +188,10 @@ 'NumPyMinimumEigensolver', 'MinimumEigensolver', 'MinimumEigensolverResult', + 'HamiltonianPhaseEstimation', + 'HamiltonianPhaseEstimationResult', + 'PhaseEstimationScale', + 'PhaseEstimation', + 'PhaseEstimationResult', 'AlgorithmError', ] diff --git a/qiskit/algorithms/phase_estimators/__init__.py b/qiskit/algorithms/phase_estimators/__init__.py new file mode 100644 index 000000000000..19567798d9af --- /dev/null +++ b/qiskit/algorithms/phase_estimators/__init__.py @@ -0,0 +1,29 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + +"""Phase Estimators.""" + +from .phase_estimator import PhaseEstimator +from .phase_estimation import PhaseEstimation +from .phase_estimation_result import PhaseEstimationResult +from .phase_estimation_scale import PhaseEstimationScale +from .hamiltonian_phase_estimation import HamiltonianPhaseEstimation +from .hamiltonian_phase_estimation_result import HamiltonianPhaseEstimationResult + +__all__ = [ + 'PhaseEstimator', + 'PhaseEstimation', + 'PhaseEstimationResult', + 'PhaseEstimationScale', + 'HamiltonianPhaseEstimation', + 'HamiltonianPhaseEstimationResult' +] diff --git a/qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation.py b/qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation.py new file mode 100644 index 000000000000..1a3e40273051 --- /dev/null +++ b/qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation.py @@ -0,0 +1,211 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + +"""Phase estimation for the spectrum of a Hamiltonian""" + +from typing import Optional, Union +from qiskit import QuantumCircuit +from qiskit.utils import QuantumInstance +from qiskit.opflow import (EvolutionBase, PauliTrotterEvolution, OperatorBase, + SummedOp, PauliOp, MatrixOp, PauliSumOp, StateFn) +from qiskit.providers import BaseBackend +from .phase_estimation import PhaseEstimation +from .hamiltonian_phase_estimation_result import HamiltonianPhaseEstimationResult +from .phase_estimation_scale import PhaseEstimationScale + + +class HamiltonianPhaseEstimation: + r"""Run the Quantum Phase Estimation algorithm to find the eigenvalues of a Hermitian operator. + + This class is nearly the same as :class:`~qiskit.algorithms.PhaseEstimation`, differing only + in that the input in that class is a unitary operator, whereas here the input is a Hermitian + operator from which a unitary will be obtained by scaling and exponentiating. The scaling is + performed in order to prevent the phases from wrapping around :math:`2\pi`. + The problem of estimating eigenvalues :math:`\lambda_j` of the Hermitian operator + :math:`H` is solved by running a circuit representing + + .. math:: + + \exp(i b H) |\psi\rangle = \sum_j \exp(i b \lambda_j) c_j |\lambda_j\rangle, + + where the input state is + + .. math:: + + |\psi\rangle = \sum_j c_j |\lambda_j\rangle, + + and :math:`\lambda_j` are the eigenvalues of :math:`H`. + + Here, :math:`b` is a scaling factor sufficiently large to map positive :math:`\lambda` to + :math:`[0,\pi)` and negative :math:`\lambda` to :math:`[\pi,2\pi)`. Each time the circuit is + run, one measures a phase corresponding to :math:`lambda_j` with probability :math:`|c_j|^2`. + + If :math:`H` is a Pauli sum, the bound :math:`b` is computed from the sum of the absolute + values of the coefficients of the terms. There is no way to reliably recover eigenvalues + from phases very near the endpoints of these intervals. Because of this you should be aware + that for degenerate cases, such as :math:`H=Z`, the eigenvalues :math:`\pm 1` will be + mapped to the same phase, :math:`\pi`, and so cannot be distinguished. In this case, you need + to specify a larger bound as an argument to the method ``estimate``. + + This class uses and works together with :class:`~qiskit.algorithms.PhaseEstimationScale` to + manage scaling the Hamiltonian and the phases that are obtained by the QPE algorithm. This + includes setting, or computing, a bound on the eigenvalues of the operator, using this + bound to obtain a scale factor, scaling the operator, and shifting and scaling the measured + phases to recover the eigenvalues. + + Note that, although we speak of "evolving" the state according the the Hamiltonian, in the + present algorithm, we are not actually considering time evolution. Rather, the role of time is + played by the scaling factor, which is chosen to best extract the eigenvalues of the + Hamiltonian. + + A few of the ideas in the algorithm may be found in Ref. [1]. + + **Reference:** + + [1]: Quantum phase estimation of multiple eigenvalues for small-scale (noisy) experiments + T.E. O'Brien, B. Tarasinski, B.M. Terhal + `arXiv:1809.09697 `_ + + """ + + def __init__(self, + num_evaluation_qubits: int, + quantum_instance: Optional[Union[QuantumInstance, BaseBackend]] = None) -> None: + """ + Args: + num_evaluation_qubits: The number of qubits used in estimating the phase. The phase will + be estimated as a binary string with this many bits. + quantum_instance: The quantum instance on which the circuit will be run. + """ + self._phase_estimation = PhaseEstimation( + num_evaluation_qubits=num_evaluation_qubits, + quantum_instance=quantum_instance) + + def _get_scale(self, hamiltonian, bound=None) -> None: + if bound is None: + return PhaseEstimationScale.from_pauli_sum(hamiltonian) + + return PhaseEstimationScale(bound) + + def _get_unitary(self, hamiltonian, pe_scale, evolution) -> QuantumCircuit: + """Evolve the Hamiltonian to obtain a unitary. + + Apply the scaling to the Hamiltonian that has been computed from an eigenvalue bound + and compute the unitary by applying the evolution object. + """ + # scale so that phase does not wrap. + scaled_hamiltonian = -pe_scale.scale * hamiltonian + unitary = evolution.convert(scaled_hamiltonian.exp_i()) + if not isinstance(unitary, QuantumCircuit): + unitary_circuit = unitary.to_circuit() + else: + unitary_circuit = unitary + + # Decomposing twice allows some 1Q Hamiltonians to give correct results + # when using MatrixEvolution(), that otherwise would give incorrect results. + # It does not break any others that we tested. + return unitary_circuit.decompose().decompose() + + # pylint: disable=arguments-differ + def estimate(self, hamiltonian: OperatorBase, + state_preparation: Optional[StateFn] = None, + evolution: Optional[EvolutionBase] = None, + bound: Optional[float] = None) -> HamiltonianPhaseEstimationResult: + """Run the Hamiltonian phase estimation algorithm. + + Args: + hamiltonian: A Hermitian operator. + state_preparation: The ``StateFn`` to be prepared, whose eigenphase will be + measured. If this parameter is omitted, no preparation circuit will be run and + input state will be the all-zero state in the computational basis. + evolution: An evolution converter that generates a unitary from ``hamiltonian``. If + ``None``, then the default ``PauliTrotterEvolution`` is used. + bound: An upper bound on the absolute value of the eigenvalues of + ``hamiltonian``. If omitted, then ``hamiltonian`` must be a Pauli sum, or a + ``PauliOp``, in which case a bound will be computed. If ``hamiltonian`` + is a ``MatrixOp``, then ``bound`` may not be ``None``. The tighter the bound, + the higher the resolution of computed phases. + + Returns: + HamiltonianPhaseEstimationResult instance containing the result of the estimation + and diagnostic information. + + Raises: + ValueError: If ``bound`` is ``None`` and ``hamiltonian`` is not a Pauli sum, i.e. a + ``PauliSumOp`` or a ``SummedOp`` whose terms are of type ``PauliOp``. + TypeError: If ``evolution`` is not of type ``EvolutionBase``. + """ + if evolution is None: + evolution = PauliTrotterEvolution() + elif not isinstance(evolution, EvolutionBase): + raise TypeError(f'Expecting type EvolutionBase, got {type(evolution)}') + + if isinstance(hamiltonian, PauliSumOp): + hamiltonian = hamiltonian.to_pauli_op() + elif isinstance(hamiltonian, PauliOp): + hamiltonian = SummedOp([hamiltonian]) + + if isinstance(hamiltonian, SummedOp): + # remove identitiy terms + # The term propto the identity is removed from hamiltonian. + # This is done for three reasons: + # 1. Work around an unknown bug that otherwise causes the energies to be wrong in some + # cases. + # 2. Allow working with a simpler Hamiltonian, one with fewer terms. + # 3. Tighten the bound on the eigenvalues so that the spectrum is better resolved, i.e. + # occupies more of the range of values representable by the qubit register. + # The coefficient of this term will be added to the eigenvalues. + id_coefficient, hamiltonian_no_id = _remove_identity(hamiltonian) + + # get the rescaling object + pe_scale = self._get_scale(hamiltonian_no_id, bound) + + # get the unitary + unitary = self._get_unitary(hamiltonian_no_id, pe_scale, evolution) + + elif isinstance(hamiltonian, MatrixOp): + if bound is None: + raise ValueError('bound must be specified if Hermitian operator is MatrixOp') + + # Do not subtract an identity term from the matrix, so do not compensate. + id_coefficient = 0.0 + pe_scale = self._get_scale(hamiltonian, bound) + unitary = self._get_unitary(hamiltonian, pe_scale, evolution) + else: + raise TypeError(f'Hermitian operator of type {type(hamiltonian)} not supported.') + + if state_preparation is not None: + state_preparation = state_preparation.to_circuit_op().to_circuit() + # run phase estimation + phase_estimation_result = self._phase_estimation.estimate( + unitary=unitary, state_preparation=state_preparation) + + return HamiltonianPhaseEstimationResult( + phase_estimation_result=phase_estimation_result, + id_coefficient=id_coefficient, + phase_estimation_scale=pe_scale) + + +def _remove_identity(pauli_sum): + """Remove any identity operators from `pauli_sum`. Return + the sum of the coefficients of the identities and the new operator. + """ + idcoeff = 0.0 + ops = [] + for op in pauli_sum: + p = op.primitive + if p.x.any() or p.z.any(): + ops.append(op) + else: + idcoeff += op.coeff + + return idcoeff, SummedOp(ops) diff --git a/qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation_result.py b/qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation_result.py new file mode 100644 index 000000000000..f2012780895a --- /dev/null +++ b/qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation_result.py @@ -0,0 +1,104 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + +"""Result from running HamiltonianPhaseEstimation""" + + +from typing import Dict, Union, cast +from qiskit.algorithms.algorithm_result import AlgorithmResult +from .phase_estimation_result import PhaseEstimationResult +from .phase_estimation_scale import PhaseEstimationScale + + +class HamiltonianPhaseEstimationResult(AlgorithmResult): + """Store and manipulate results from running `HamiltonianPhaseEstimation`. + + This API of this class is nearly the same as `PhaseEstimatorResult`, differing only in + the presence of an additional keyword argument in the methods. If `scaled` + is `False`, then the phases are not translated and scaled to recover the + eigenvalues of the Hamiltonian. Instead `phi` in :math:`[0, 1)` is returned, + as is the case when then unitary is not derived from a Hamiltonian. + + This class is meant to be instantiated via `HamiltonianPhaseEstimation.estimate`. + """ + + def __init__(self, + phase_estimation_result: PhaseEstimationResult, + phase_estimation_scale: PhaseEstimationScale, + id_coefficient: float, + ) -> None: + """ + Args: + phase_estimation_result: The result object returned by PhaseEstimation.estimate. + phase_estimation_scale: object used to scale phases to obtain eigenvalues. + id_coefficient: The coefficient of the identity term in the Hamiltonian. + Eigenvalues are computed without this term so that the + coefficient must added to give correct eigenvalues. + This is done automatically when retrieving eigenvalues. + """ + super().__init__() + self._phase_estimation_scale = phase_estimation_scale + self._id_coefficient = id_coefficient + self._phase_estimation_result = phase_estimation_result + + # pylint: disable=arguments-differ + def filter_phases(self, cutoff: float = 0.0, scaled: bool = True, + as_float: bool = True) -> Dict[Union[str, float], float]: + """Filter phases as does `PhaseEstimatorResult.filter_phases`, with + the addition that `phi` is shifted and translated to return eigenvalues + of the Hamiltonian. + + Args: + cutoff: Minimum weight of number of counts required to keep a bit string. + The default value is `0.0`. + scaled: If False, return `phi` in :math:`[0, 1)` rather than the eigenvalues of + the Hamiltonian. + as_float: If `True`, returned keys are floats in :math:`[0.0, 1.0)`. If `False` + returned keys are bit strings. + + Raises: + ValueError: if `as_float` is `False` and `scaled` is `True`. + + Returns: + A dict of filtered phases. + """ + if scaled and not as_float: + raise ValueError('`as_float` must be `True` if `scaled` is `True`.') + + phases = self._phase_estimation_result.filter_phases(cutoff, as_float=as_float) + if scaled: + return cast(Dict, self._phase_estimation_scale.scale_phases(phases, + self._id_coefficient)) + else: + return cast(Dict, phases) + + @property + def most_likely_phase(self) -> float: + """The most likely phase of the unitary corresponding to the Hamiltonian. + + Returns: + The most likely phase. + """ + return self._phase_estimation_result.most_likely_phase + + @property + def most_likely_eigenvalue(self) -> float: + """The most likely eigenvalue of the Hamiltonian. + + This method calls `most_likely_phase` and scales the result to + obtain an eigenvalue. + + Returns: + The most likely eigenvalue of the Hamiltonian. + """ + phase = self._phase_estimation_result.most_likely_phase + return self._phase_estimation_scale.scale_phase(phase, self._id_coefficient) diff --git a/qiskit/algorithms/phase_estimators/phase_estimation.py b/qiskit/algorithms/phase_estimators/phase_estimation.py new file mode 100644 index 000000000000..1d4e542c6b10 --- /dev/null +++ b/qiskit/algorithms/phase_estimators/phase_estimation.py @@ -0,0 +1,225 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + + +"""The Quantum Phase Estimation Algorithm.""" + + +from typing import Optional, Union +import numpy +from qiskit.circuit import QuantumCircuit +import qiskit +import qiskit.circuit as circuit +from qiskit.circuit.classicalregister import ClassicalRegister +from qiskit.providers import BaseBackend, Backend +from qiskit.utils import QuantumInstance +from qiskit.result import Result +from .phase_estimation_result import PhaseEstimationResult, _sort_phases +from .phase_estimator import PhaseEstimator + + +class PhaseEstimation(PhaseEstimator): + r"""Run the Quantum Phase Estimation (QPE) algorithm. + + This runs QPE with a multi-qubit register for reading the phases [1] + of input states. + + The algorithm takes as input a unitary :math:`U` and a state :math:`|\psi\rangle`, + which may be written + + .. math:: + + |\psi\rangle = \sum_j c_j |\phi_j\rangle, + + where :math:`|\phi_j\rangle` are eigenstates of :math:`U`. We prepare the quantum register + in the state :math:`|\psi\rangle` then apply :math:`U` leaving the register in the state + + .. math:: + + U|\psi\rangle = \sum_j \exp(i \phi_j) c_j |\phi_j\rangle. + + In the ideal case, one then measures the phase :math:`\phi_j` with probability + :math:`|c_j|^2`. In practice, many (or all) of the bit strings may be measured due to + noise and the possibility that :math:`\phi_j` may not be representable exactly by the + output register. In the latter case the probability for each eigenphase will be spread + across bitstrings, with amplitudes that decrease with distance from the bitstring most + closely approximating the eigenphase. + + The main inputs are the number of qubits in the phase-reading register, a state preparation + circuit to prepare an input state, and either + + 1) A unitary that will act on the the input state, or + 2) A quantum-phase-estimation circuit in which the unitary is already embedded. + + In case 1), an instance of :class:`qiskit.circuit.PhaseEstimation`, a QPE circuit, containing + the input unitary will be constructed. After construction, the QPE circuit is run on a backend + via the `run` method, and the frequencies or counts of the phases represented by bitstrings + are recorded. The results are returned as an instance of + :class:`~qiskit.algorithms.phase_estimator_result.PhaseEstimationResult`. + + **Reference:** + + [1]: Michael A. Nielsen and Isaac L. Chuang. 2011. + Quantum Computation and Quantum Information: 10th Anniversary Edition (10th ed.). + Cambridge University Press, New York, NY, USA. + + """ + + def __init__(self, + num_evaluation_qubits: int, + quantum_instance: Optional[Union[QuantumInstance, + BaseBackend, Backend]] = None) -> None: + """ + Args: + num_evaluation_qubits: The number of qubits used in estimating the phase. The phase will + be estimated as a binary string with this many bits. + quantum_instance: The quantum instance on which the circuit will be run. + """ + + self._measurements_added = False + if num_evaluation_qubits is not None: + self._num_evaluation_qubits = num_evaluation_qubits + + if isinstance(quantum_instance, (Backend, BaseBackend)): + quantum_instance = QuantumInstance(quantum_instance) + self._quantum_instance = quantum_instance + + def construct_circuit(self, + unitary: QuantumCircuit, + state_preparation: Optional[QuantumCircuit] = None) -> QuantumCircuit: + """Return the circuit to be executed to estimate phases. + + This circuit includes as sub-circuits the core phase estimation circuit, + with the addition of the state-preparation circuit and possibly measurement instructions. + """ + num_evaluation_qubits = self._num_evaluation_qubits + num_unitary_qubits = unitary.num_qubits + + pe_circuit = circuit.library.PhaseEstimation(num_evaluation_qubits, unitary) + + if state_preparation is not None: + pe_circuit.compose( + state_preparation, + qubits=range(num_evaluation_qubits, + num_evaluation_qubits + num_unitary_qubits), + inplace=True, + front=True) + + self._add_measurement_if_required(pe_circuit) + + return pe_circuit + + def _add_measurement_if_required(self, pe_circuit): + if not self._quantum_instance.is_statevector: + # Measure only the evaluation qubits. + regname = 'meas' + creg = ClassicalRegister(self._num_evaluation_qubits, regname) + pe_circuit.add_register(creg) + pe_circuit.barrier() + pe_circuit.measure(range(self._num_evaluation_qubits), + range(self._num_evaluation_qubits)) + + return circuit + + def _compute_phases(self, + num_unitary_qubits: int, + circuit_result: Result) -> Union[numpy.ndarray, qiskit.result.Counts]: + """Compute frequencies/counts of phases from the result of running the QPE circuit. + + How the frequencies are computed depends on whether the backend computes amplitude or + samples outcomes. + + 1) If the backend is a statevector simulator, then the reduced density matrix of the + phase-reading register is computed from the combined phase-reading- and input-state + registers. The elements of the diagonal :math:`(i, i)` give the probability to measure the + each of the states `i`. The index `i` expressed as a binary integer with the LSB rightmost + gives the state of the phase-reading register with the LSB leftmost when interpreted as a + phase. In order to maintain the compact representation, the phases are maintained as decimal + integers. They may be converted to other forms via the results object, + `PhaseEstimationResult` or `HamiltonianPhaseEstimationResult`. + + 2) If the backend samples bitstrings, then the counts are first retrieved as a dict. The + binary strings (the keys) are then reversed so that the LSB is rightmost and the counts are + converted to frequencies. Then the keys are sorted according to increasing phase, so that + they can be easily understood when displaying or plotting a histogram. + + Args: + num_unitary_qubits: The number of qubits in the unitary. + circuit_result: the result object returned by the backend that ran the QPE circuit. + + Returns: + Either a dict or numpy.ndarray representing the frequencies of the phases. + + """ + if self._quantum_instance.is_statevector: + state_vec = circuit_result.get_statevector() + evaluation_density_matrix = qiskit.quantum_info.partial_trace( + state_vec, + range(self._num_evaluation_qubits, + self._num_evaluation_qubits + num_unitary_qubits) + ) + phases = evaluation_density_matrix.probabilities() + else: + # return counts with keys sorted numerically + num_shots = circuit_result.results[0].shots + counts = circuit_result.get_counts() + phases = {k[::-1]: counts[k] / num_shots for k in counts.keys()} + phases = _sort_phases(phases) + phases = qiskit.result.Counts(phases, memory_slots=counts.memory_slots, + creg_sizes=counts.creg_sizes) + + return phases + + def estimate(self, + unitary: Optional[QuantumCircuit] = None, + state_preparation: Optional[QuantumCircuit] = None, + pe_circuit: Optional[QuantumCircuit] = None, + num_unitary_qubits: Optional[int] = None) -> PhaseEstimationResult: + """Run the the phase estimation algorithm. + + Args: + unitary: The circuit representing the unitary operator whose eigenvalues (via phase) + will be measured. Exactly one of `pe_circuit` and `unitary` must be passed. + state_preparation: The circuit that prepares the state whose eigenphase will be + measured. If this parameter is omitted, no preparation circuit + will be run and input state will be the all-zero state in the + computational basis. + pe_circuit: The phase estimation circuit. + num_unitary_qubits: Must agree with the number of qubits in the unitary in `pe_circuit` + if `pe_circuit` is passed. This parameter will be set from `unitary` + if `unitary` is passed. + + Raises: + ValueError: If both `pe_circuit` and `unitary` are passed. + ValueError: If neither `pe_circuit` nor `unitary` are passed. + + Returns: + An instance of qiskit.algorithms.phase_estimator_result.PhaseEstimationResult. + """ + num_evaluation_qubits = self._num_evaluation_qubits + + if unitary is not None: + if pe_circuit is not None: + raise ValueError('Only one of `pe_circuit` and `unitary` may be passed.') + pe_circuit = self.construct_circuit(unitary, state_preparation) + num_unitary_qubits = unitary.num_qubits + + elif pe_circuit is not None: + self._add_measurement_if_required(pe_circuit) + + else: + raise ValueError('One of `pe_circuit` and `unitary` must be passed.') + + circuit_result = self._quantum_instance.execute(pe_circuit) + phases = self._compute_phases(num_unitary_qubits, circuit_result) + return PhaseEstimationResult(num_evaluation_qubits, circuit_result=circuit_result, + phases=phases) diff --git a/qiskit/algorithms/phase_estimators/phase_estimation_result.py b/qiskit/algorithms/phase_estimators/phase_estimation_result.py new file mode 100644 index 000000000000..0322fc24e94f --- /dev/null +++ b/qiskit/algorithms/phase_estimators/phase_estimation_result.py @@ -0,0 +1,154 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + +"""Result of running PhaseEstimation""" + +from typing import Dict, Union +import numpy +from qiskit.result import Result +from .phase_estimator import PhaseEstimatorResult + + +class PhaseEstimationResult(PhaseEstimatorResult): + """Store and manipulate results from running `PhaseEstimation`. + + This class is instantiated by the `PhaseEstimation` class, not via user code. + The `PhaseEstimation` class generates a list of phases and corresponding weights. Upon + completion it returns the results as an instance of this class. The main method for + accessing the results is `filter_phases`. + """ + + def __init__(self, num_evaluation_qubits: int, + circuit_result: Result, + phases: Union[numpy.ndarray, Dict[str, float]]) -> None: + """ + Args: + num_evaluation_qubits: number of qubits in phase-readout register. + circuit_result: result object returned by method running circuit. + phases: ndarray or dict of phases and frequencies determined by QPE. + """ + super().__init__() + + self._phases = phases + # int: number of qubits in phase-readout register + self._num_evaluation_qubits = num_evaluation_qubits + self._circuit_result = circuit_result + + @property + def phases(self) -> Union[numpy.ndarray, dict]: + """Return all phases and their frequencies computed by QPE. + + This is an array or dict whose values correspond to weights on bit strings. + """ + return self._phases + + @property + def circuit_result(self) -> Result: + """Return the result object returned by running the QPE circuit (on hardware or simulator). + + This is useful for inspecting and troubleshooting the QPE algorithm. + """ + return self._circuit_result + + @property + def most_likely_phase(self) -> float: + r"""Return the estimated phase as a number in :math:`[0.0, 1.0)`. + + 1.0 corresponds to a phase of :math:`2\pi`. It is assumed that the input vector is an + eigenvector of the unitary so that the peak of the probability density occurs at the bit + string that most closely approximates the true phase. + """ + if isinstance(self.phases, dict): + binary_phase_string = max(self.phases, key=self.phases.get) + else: + # numpy.argmax ignores complex part of number. But, we take abs anyway + idx = numpy.argmax(abs(self.phases)) + binary_phase_string = numpy.binary_repr(idx, self._num_evaluation_qubits)[::-1] + phase = _bit_string_to_phase(binary_phase_string) + return phase + + def filter_phases(self, cutoff: float = 0.0, + as_float: bool = True) -> Dict: + """Return a filtered dict of phases (keys) and frequencies (values). + + Only phases with frequencies (counts) larger than `cutoff` are included. + It is assumed that the `run` method has been called so that the phases have been computed. + When using a noiseless, shot-based simulator to read a single phase that can + be represented exactly by `num_evaluation_qubits`, all the weight will + be concentrated on a single phase. In all other cases, many, or all, bit + strings will have non-zero weight. This method is useful for filtering + out these uninteresting bit strings. + + Args: + cutoff: Minimum weight of number of counts required to keep a bit string. + The default value is `0.0`. + as_float: If `True`, returned keys are floats in :math:`[0.0, 1.0)`. If `False` + returned keys are bit strings. + + Returns: + A filtered dict of phases (keys) and frequencies (values). + """ + if isinstance(self.phases, dict): + counts = self.phases + if as_float: + phases = {_bit_string_to_phase(k): counts[k] + for k in counts.keys() if counts[k] > cutoff} + else: + phases = {k: counts[k] for k in counts.keys() if counts[k] > cutoff} + + else: + phases = {} + for idx, amplitude in enumerate(self.phases): + if amplitude > cutoff: + # Each index corresponds to a computational basis state with the LSB rightmost. + # But, we chose to apply the unitaries such that the phase is recorded + # in reverse order. So, we reverse the bitstrings here. + binary_phase_string = numpy.binary_repr(idx, self._num_evaluation_qubits)[::-1] + if as_float: + _key = _bit_string_to_phase(binary_phase_string) + else: + _key = binary_phase_string + phases[_key] = amplitude + + phases = _sort_phases(phases) + + return phases + + +def _bit_string_to_phase(binary_string: str) -> float: + """Convert bit string to a normalized phase in :math:`[0,1)`. + + It is assumed that the bit string is correctly padded and that the order of + the bits has been reversed relative to their order when the counts + were recorded. The LSB is the right most when interpreting the bitstring as + a phase. + + Args: + binary_string: A string of characters '0' and '1'. + + Returns: + A phase scaled to :math:`[0,1)`. + """ + n_qubits = len(binary_string) + return int(binary_string, 2) / (2 ** n_qubits) + + +def _sort_phases(phases: Dict) -> Dict: + """Sort a dict of bit strings representing phases (keys) and frequencies (values) by bit string. + + The bit strings are sorted according to increasing phase. This relies on Python + preserving insertion order when building dicts. + """ + pkeys = list(phases.keys()) + pkeys.sort(reverse=False) # Sorts in order of the integer encoded by binary string + phases = {k: phases[k] for k in pkeys} + return phases diff --git a/qiskit/algorithms/phase_estimators/phase_estimation_scale.py b/qiskit/algorithms/phase_estimators/phase_estimation_scale.py new file mode 100644 index 000000000000..d3509e4d56cf --- /dev/null +++ b/qiskit/algorithms/phase_estimators/phase_estimation_scale.py @@ -0,0 +1,141 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + +"""Scaling for Hamiltonian and eigenvalues to avoid phase wrapping""" + +from typing import Union, Dict, List +import numpy +from qiskit.opflow import SummedOp + + +class PhaseEstimationScale(): + """Set and use a bound on eigenvalues of a Hermitian operator in order to ensure phases are in + the desired range and to convert measured phases into eigenvectors. + + The ``bound`` is set when constructing this class. Then the method ``scale`` is used to find the + factor by which to scale the operator. + + If ``bound`` is equal exactly to the largest eigenvalue, and the smallest eigenvalue is minus + the largest, then these two eigenvalues will not be distinguished. For example, if the Hermitian + operator is the Pauli Z operator with eigenvalues :math:`1` and :math:`-1`, and ``bound`` is + :math:`1`, then both eigenvalues will be mapped to :math:`1`. + This can be avoided by making ``bound`` a bit larger. + + Increasing ``bound`` decreases the part of the interval :math:`[0, 1)` that is used to map + eigenvalues to ``phi``. However, sometimes this results in a better determination of the + eigenvalues, because 1) although there are fewer discrete phases in the useful range, it may + shift one of the discrete phases closer to the actual phase. And, 2) If one of the discrete + phases is close to, or exactly equal to the actual phase, then artifacts (probability) in + neighboring phases will be reduced. This is important because the artifacts may be larger than + the probability in a phase representing another eigenvalue of interest whose corresponding + eigenstate has a relatively small weight in the input state. + + """ + + def __init__(self, bound: float) -> None: + """ + Args: + bound: an upper bound on the absolute value of the eigenvalues of a Hermitian operator. + (The operator is not needed here.) + """ + self._bound = bound + + @property + def scale(self) -> float: + r"""Return the Hamiltonian scaling factor. + + Return the scale factor by which a Hermitian operator must be multiplied + so that the phase of the corresponding unitary is restricted to :math:`[-\pi, \pi]`. + This factor is computed from the bound on the absolute values of the eigenvalues + of the operator. The methods ``scale_phase`` and ``scale_phases`` are used recover + the eigenvalues corresponding the original (unscaled) Hermitian operator. + + Returns: + The scale factor. + """ + return numpy.pi / self._bound + + def scale_phase(self, phi: float, id_coefficient: float = 0.0) -> float: + r"""Convert a phase into an eigenvalue. + + The input phase ``phi`` corresponds to the eigenvalue of a unitary obtained by + exponentiating a scaled Hermitian operator. Recall that the phase + is obtained from ``phi`` as :math:`2\pi\phi`. Furthermore, the Hermitian operator + was scaled so that ``phi`` is restricted to :math:`[-1/2, 1/2]`, corresponding to + phases in :math:`[-\pi, \pi]`. But the values of `phi` read from the phase-readout + register are in :math:`[0, 1)`. Any value of ``phi`` greater than :math:`1/2` corresponds + to a raw phase of minus the complement with respect to 1. After this possible + shift, the phase is scaled by the inverse of the factor by which the + Hermitian operator was scaled to recover the eigenvalue of the Hermitian + operator. + + Args: + phi: Normalized phase in :math:`[0, 1)` to be converted to an eigenvalue. + id_coefficient: All eigenvalues are shifted by this value. + + Returns: + An eigenvalue computed from the input phase. + """ + w = 2 * self._bound + if phi <= 0.5: + return phi * w + id_coefficient + else: + return (phi - 1) * w + id_coefficient + + # pylint: disable=unsubscriptable-object + def scale_phases(self, phases: Union[List, Dict], id_coefficient: float = 0.0 + ) -> Union[Dict, List]: + """Convert a list or dict of phases to eigenvalues. + + The values in the list, or keys in the dict, are values of ``phi` and + are converted as described in the description of ``scale_phase``. In case + ``phases`` is a dict, the values of the dict are passed unchanged. + + Args: + phases: a list or dict of values of ``phi``. + id_coefficient: All eigenvalues are shifted by this value. + + Returns: + Eigenvalues computed from phases. + """ + if isinstance(phases, list): + phases = [self.scale_phase(x, id_coefficient) for x in phases] + else: + phases = {self.scale_phase(x, id_coefficient): phases[x] for x in phases.keys()} + + return phases + + @classmethod + def from_pauli_sum(cls, pauli_sum: SummedOp) -> 'PhaseEstimationScale': + """Create a PhaseEstimationScale from a `SummedOp` representing a sum of Pauli Operators. + + It is assumed that the ``pauli_sum`` is the sum of ``PauliOp`` objects. The bound on + the absolute value of the eigenvalues of the sum is obtained as the sum of the + absolute values of the coefficients of the terms. This is the best bound available in + the generic case. A ``PhaseEstimationScale`` object is instantiated using this bound. + + Args: + pauli_sum: A ``SummedOp`` whose terms are ``PauliOp`` objects. + + Raises: + ValueError: if ``pauli_sum`` is not a sum of Pauli operators. + + Returns: + A ``PhaseEstimationScale`` object + """ + if pauli_sum.primitive_strings() != {'Pauli'}: + raise ValueError( + '`pauli_sum` must be a sum of Pauli operators. Got primitives {}.'.format( + pauli_sum.primitive_strings())) + + bound = sum([abs(pauli_sum.coeff * pauli.coeff) for pauli in pauli_sum]) + return PhaseEstimationScale(bound) diff --git a/qiskit/algorithms/phase_estimators/phase_estimator.py b/qiskit/algorithms/phase_estimators/phase_estimator.py new file mode 100644 index 000000000000..31af0b3ad976 --- /dev/null +++ b/qiskit/algorithms/phase_estimators/phase_estimator.py @@ -0,0 +1,50 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + +"""The Phase Estimator interface.""" + +from typing import Optional +from abc import ABC, abstractmethod, abstractproperty +from qiskit.circuit import QuantumCircuit +from qiskit.algorithms.algorithm_result import AlgorithmResult + + +class PhaseEstimator(ABC): + """The Phase Estimator interface. + + Algorithms that can compute a phase for a unitary operator and + initial state may implement this interface to allow different + algorithms to be used interchangeably. + """ + + @abstractmethod + def estimate(self, + unitary: Optional[QuantumCircuit] = None, + state_preparation: Optional[QuantumCircuit] = None, + pe_circuit: Optional[QuantumCircuit] = None, + num_unitary_qubits: Optional[int] = None) -> 'PhaseEstimatorResult': + """Estimate the phase.""" + raise NotImplementedError + + +class PhaseEstimatorResult(AlgorithmResult): + """Phase Estimator Result.""" + + @abstractproperty + def most_likely_phase(self) -> float: + r"""Return the estimated phase as a number in :math:`[0.0, 1.0)`. + + 1.0 corresponds to a phase of :math:`2\pi`. It is assumed that the input vector is an + eigenvector of the unitary so that the peak of the probability density occurs at the bit + string that most closely approximates the true phase. + """ + raise NotImplementedError diff --git a/test/python/algorithms/test_phase_estimator.py b/test/python/algorithms/test_phase_estimator.py new file mode 100644 index 000000000000..4cb48bf34a46 --- /dev/null +++ b/test/python/algorithms/test_phase_estimator.py @@ -0,0 +1,254 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2020. +# +# 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. + +"""Test phase estimation""" + +import unittest +from test.python.algorithms import QiskitAlgorithmsTestCase +from ddt import ddt, data +import numpy as np +from qiskit.algorithms.phase_estimators import PhaseEstimation, HamiltonianPhaseEstimation +from qiskit.opflow.evolutions import PauliTrotterEvolution, MatrixEvolution +import qiskit +from qiskit.opflow import (H, X, Y, Z, I, StateFn) + + +@ddt +class TestHamiltonianPhaseEstimation(QiskitAlgorithmsTestCase): + """Tests for obtaining eigenvalues from phase estimation""" + + def hamiltonian_pe(self, hamiltonian, state_preparation=None, num_evaluation_qubits=6, + backend=None, + evolution=None, + bound=None): + """Run HamiltonianPhaseEstimation and return result with all phases.""" + if backend is None: + backend = qiskit.BasicAer.get_backend('statevector_simulator') + quantum_instance = qiskit.utils.QuantumInstance(backend=backend, shots=10000) + phase_est = HamiltonianPhaseEstimation( + num_evaluation_qubits=num_evaluation_qubits, + quantum_instance=quantum_instance) + result = phase_est.estimate( + hamiltonian=hamiltonian, + state_preparation=state_preparation, evolution=evolution, + bound=bound) + return result + + @data(MatrixEvolution(), PauliTrotterEvolution('suzuki', 4)) + def test_pauli_sum_1(self, evolution): + """Two eigenvalues from Pauli sum with X, Z""" + hamiltonian = 0.5 * X + Z + state_preparation = StateFn(H.to_circuit()) + + result = self.hamiltonian_pe(hamiltonian, state_preparation, evolution=evolution) + phase_dict = result.filter_phases(0.162, as_float=True) + phases = list(phase_dict.keys()) + phases.sort() + + self.assertAlmostEqual(phases[0], -1.125, delta=0.001) + self.assertAlmostEqual(phases[1], 1.125, delta=0.001) + + @data(MatrixEvolution(), PauliTrotterEvolution('suzuki', 3)) + def test_pauli_sum_2(self, evolution): + """Two eigenvalues from Pauli sum with X, Y, Z""" + hamiltonian = 0.5 * X + Y + Z + state_preparation = None + + result = self.hamiltonian_pe(hamiltonian, state_preparation, evolution=evolution) + phase_dict = result.filter_phases(0.1, as_float=True) + phases = list(phase_dict.keys()) + phases.sort() + + self.assertAlmostEqual(phases[0], -1.484, delta=0.001) + self.assertAlmostEqual(phases[1], 1.484, delta=0.001) + + def test_single_pauli_op(self): + """Two eigenvalues from Pauli sum with X, Y, Z""" + hamiltonian = Z + state_preparation = None + + result = self.hamiltonian_pe(hamiltonian, state_preparation, evolution=None) + eigv = result.most_likely_eigenvalue + with self.subTest('First eigenvalue'): + self.assertAlmostEqual(eigv, 1.0, delta=0.001) + + state_preparation = StateFn(X.to_circuit()) + + result = self.hamiltonian_pe(hamiltonian, state_preparation, bound=1.05) + eigv = result.most_likely_eigenvalue + with self.subTest('Second eigenvalue'): + self.assertAlmostEqual(eigv, -0.98, delta=0.01) + + def test_H2_hamiltonian(self): + """Test H2 hamiltonian""" + hamiltonian = (-1.0523732457728587 * (I ^ I)) + (0.3979374248431802 * (I ^ Z)) \ + + (-0.3979374248431802 * (Z ^ I)) + (-0.011280104256235324 * (Z ^ Z)) \ + + (0.18093119978423147 * (X ^ X)) + state_preparation = StateFn((I ^ H).to_circuit()) + evo = PauliTrotterEvolution(trotter_mode='suzuki', reps=4) + result = self.hamiltonian_pe(hamiltonian, state_preparation, evolution=evo) + with self.subTest('Most likely eigenvalues'): + self.assertAlmostEqual(result.most_likely_eigenvalue, -1.855, delta=.001) + with self.subTest('Most likely phase'): + self.assertAlmostEqual(result.most_likely_phase, 0.5937, delta=.001) + with self.subTest('All eigenvalues'): + phase_dict = result.filter_phases(0.1) + phases = list(phase_dict.keys()) + self.assertAlmostEqual(phases[0], -0.8979, delta=0.001) + self.assertAlmostEqual(phases[1], -1.8551, delta=0.001) + self.assertAlmostEqual(phases[2], -1.2376, delta=0.001) + + def test_matrix_evolution(self): + """1Q Hamiltonian with MatrixEvolution""" + hamiltonian = ((0.5 * X) + (0.6 * Y) + (0.7 * I)) + state_preparation = None + result = self.hamiltonian_pe(hamiltonian, state_preparation, evolution=MatrixEvolution()) + phase_dict = result.filter_phases(0.2, as_float=True) + phases = list(phase_dict.keys()) + self.assertAlmostEqual(phases[0], 1.490, delta=0.001) + self.assertAlmostEqual(phases[1], -0.090, delta=0.001) + + def _setup_from_bound(self, evolution, op_class): + hamiltonian = 0.5 * X + Y + Z + state_preparation = None + bound = 1.2 * sum([abs(hamiltonian.coeff * coeff) for coeff in hamiltonian.coeffs]) + if op_class == 'MatrixOp': + hamiltonian = hamiltonian.to_matrix_op() + backend = qiskit.BasicAer.get_backend('statevector_simulator') + qi = qiskit.utils.QuantumInstance(backend=backend, shots=10000) + phase_est = HamiltonianPhaseEstimation(num_evaluation_qubits=6, + quantum_instance=qi) + result = phase_est.estimate(hamiltonian=hamiltonian, + bound=bound, + evolution=evolution, + state_preparation=state_preparation) + return result + + def test_from_bound(self): + """HamiltonianPhaseEstimation with bound""" + for op_class in ('SummedOp', 'MatrixOp'): + result = self._setup_from_bound(MatrixEvolution(), op_class) + cutoff = 0.01 + phases = result.filter_phases(cutoff) + with self.subTest(f'test phases has the correct length: {op_class}'): + self.assertEqual(len(phases), 2) + with self.subTest(f'test scaled phases are correct: {op_class}'): + self.assertEqual(list(phases.keys()), [1.5, -1.5]) + phases = result.filter_phases(cutoff, scaled=False) + with self.subTest(f'test unscaled phases are correct: {op_class}'): + self.assertEqual(list(phases.keys()), [0.25, 0.75]) + + def test_trotter_from_bound(self): + """HamiltonianPhaseEstimation with bound and Trotterization""" + result = self._setup_from_bound(PauliTrotterEvolution(trotter_mode='suzuki', reps=3), + op_class='SummedOp') + phase_dict = result.filter_phases(0.1) + phases = list(phase_dict.keys()) + with self.subTest('test phases has the correct length'): + self.assertEqual(len(phases), 2) + with self.subTest('test phases has correct values'): + self.assertAlmostEqual(phases[0], 1.5, delta=0.001) + self.assertAlmostEqual(phases[1], -1.5, delta=0.001) + + +@ddt +class TestPhaseEstimation(QiskitAlgorithmsTestCase): + """Evolution tests.""" + + # pylint: disable=invalid-name + def one_phase(self, unitary_circuit, state_preparation=None, n_eval_qubits=6, + backend=None): + """Run phase estimation with operator, eigenvalue pair `unitary_circuit`, + `state_preparation`. Return the bit string with the largest amplitude. + """ + if backend is None: + backend = qiskit.BasicAer.get_backend('qasm_simulator') + qi = qiskit.utils.QuantumInstance(backend=backend, shots=10000) + p_est = PhaseEstimation(num_evaluation_qubits=n_eval_qubits, + quantum_instance=qi) + result = p_est.estimate(unitary=unitary_circuit, + state_preparation=state_preparation) + phase = result.most_likely_phase + return phase + + @data('qasm_simulator', 'statevector_simulator') + def test_qpe_Z0(self, backend_type): + """eigenproblem Z, |0>""" + backend = qiskit.BasicAer.get_backend(backend_type) + + unitary_circuit = Z.to_circuit() + state_preparation = None # prepare |0> + phase = self.one_phase(unitary_circuit, state_preparation, backend=backend) + self.assertEqual(phase, 0.0) + + @data('qasm_simulator', 'statevector_simulator') + def test_qpe_Z1(self, backend_type): + """eigenproblem Z, |1>""" + backend = qiskit.BasicAer.get_backend(backend_type) + + unitary_circuit = Z.to_circuit() + state_preparation = X.to_circuit() # prepare |1> + phase = self.one_phase(unitary_circuit, state_preparation, backend=backend) + self.assertEqual(phase, 0.5) + + @data('plus', 'minus') + def test_qpe_Xplus(self, state): + """eigenproblem X, |+>""" + unitary_circuit = X.to_circuit() + if state == 'minus': # prepare |-> + state_preparation = X.to_circuit() + state_preparation.h(0) + else: # prepare |+> + state_preparation = H.to_circuit() + + phase = self.one_phase(unitary_circuit, state_preparation) + if state == 'minus': + self.assertEqual(phase, 0.5) + else: + self.assertEqual(phase, 0.0) + + def phase_estimation(self, unitary_circuit, state_preparation=None, num_evaluation_qubits=6, + backend=None): + """Run phase estimation with operator, eigenvalue pair `unitary_circuit`, + `state_preparation`. Return all results + """ + if backend is None: + backend = qiskit.BasicAer.get_backend('statevector_simulator') + qi = qiskit.utils.QuantumInstance(backend=backend, shots=10000) + phase_est = PhaseEstimation(num_evaluation_qubits=num_evaluation_qubits, + quantum_instance=qi) + result = phase_est.estimate(unitary=unitary_circuit, + state_preparation=state_preparation) + return result + + def test_qpe_Zplus(self): + """superposition eigenproblem Z, |+>""" + unitary_circuit = Z.to_circuit() + state_preparation = H.to_circuit() # prepare |+> + result = self.phase_estimation( + unitary_circuit, state_preparation, + backend=qiskit.BasicAer.get_backend('statevector_simulator')) + + phases = result.filter_phases(1e-15, as_float=True) + with self.subTest('test phases has correct values'): + self.assertEqual(list(phases.keys()), [0.0, 0.5]) + + with self.subTest('test phases has correct probabilities'): + np.testing.assert_allclose(list(phases.values()), [0.5, 0.5]) + + with self.subTest('test bitstring representation'): + phases = result.filter_phases(1e-15, as_float=False) + self.assertEqual(list(phases.keys()), ['000000', '100000']) + + +if __name__ == '__main__': + unittest.main()