Implement quantum phase estimation algorithms

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__ = [
@@ -171,5 +187,10 @@
     'NumPyMinimumEigensolver',
     'MinimumEigensolver',
     'MinimumEigensolverResult',
+    'HamiltonianPhaseEstimation',
+    'HamiltonianPhaseEstimationResult',
+    'PhaseEstimationScale',
+    'PhaseEstimation',
+    'PhaseEstimationResult',
     'AlgorithmError',
 ]

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'
+]

qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation.py

@@ -0,0 +1,150 @@
+# 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, OperatorBase, SummedOp
+from qiskit.providers import BaseBackend, Backend
+from .phase_estimation import PhaseEstimation
+from . import phase_estimation_scale
+from .hamiltonian_phase_estimation_result import HamiltonianPhaseEstimationResult
+from .phase_estimation_scale import PhaseEstimationScale
+
+
+class HamiltonianPhaseEstimation(PhaseEstimation):
+    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.PhaseEstimator`, 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`. 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
+         https://arxiv.org/abs/1809.09697
+    """
+    def __init__(self,
+                 num_evaluation_qubits: int,
+                 hamiltonian: OperatorBase,
+                 evolution: EvolutionBase,
+                 state_preparation: Optional[QuantumCircuit] = None,
+                 bound: Optional[float] = None,
+                 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.
+            hamiltonian: a Hamiltonian or Hermitian operator
+            evolution: An evolution object that generates a unitary from `hamiltonian`.
+            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.
+            bound: An upper bound on the absolute value of the eigenvalues of
+                `hamiltonian`. If omitted, then `hamiltonian` must be a Pauli sum, in which case
+                then a bound will be computed.
+            quantum_instance: The quantum instance on which the circuit will be run.
+
+        Raises:
+            ValueError: if `bound` is `None` and `hamiltonian` is not a Pauli sum (i.e. a
+            `SummedOp` whose terms are `PauliOp`s.)
+        """
+
+        self._evolution = evolution
+        self._bound = bound
+
+        # The term propto the identity is removed from hamiltonian.
+        # The coefficient of this term will be added to the eigenvalues.
+        id_coefficient, hamiltonian_no_id = _remove_identity(hamiltonian)
+        self._hamiltonian = hamiltonian_no_id
+        self._id_coefficient = id_coefficient
+
+        self._set_scale()
+        unitary = self._get_unitary()
+
+        super().__init__(num_evaluation_qubits,
+                         unitary=unitary,
+                         pe_circuit=None,
+                         num_unitary_qubits=None,
+                         state_preparation=state_preparation,
+                         quantum_instance=quantum_instance)
+
+    def _set_scale(self) -> None:
+        if self._bound is None:
+            pe_scale = phase_estimation_scale.from_pauli_sum(self._hamiltonian)
+            self._pe_scale = pe_scale
+        else:
+            self._pe_scale = PhaseEstimationScale(self._bound)
+
+    def _get_unitary(self) -> 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 = -self._pe_scale.scale * self._hamiltonian
+        unitary = self._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()
+
+    def _run(self) -> HamiltonianPhaseEstimationResult:
+        """Run the circuit and return and return `HamiltonianPhaseEstimationResult`.
+        """
+
+        circuit_result = self._quantum_instance.execute(self.construct_circuit())
+        phases = self._compute_phases(circuit_result)
+        return HamiltonianPhaseEstimationResult(
+            self._num_evaluation_qubits, phases=phases, id_coefficient=self._id_coefficient,
+            circuit_result=circuit_result, phase_estimation_scale=self._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)

qiskit/algorithms/phase_estimators/hamiltonian_phase_estimation_result.py

@@ -0,0 +1,94 @@
+# 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
+import numpy
+from qiskit.result import Result
+from .phase_estimation_result import PhaseEstimationResult
+from .phase_estimation_scale import PhaseEstimationScale
+
+
+class HamiltonianPhaseEstimationResult(PhaseEstimationResult):
+    """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.
+    """
+    def __init__(self, num_evaluation_qubits: int, circuit_result: Result,
+                 phase_estimation_scale: PhaseEstimationScale,
+                 id_coefficient: float,
+                 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.
+            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.
+            phases: ndarray or dict of phases and frequencies determined by QPE.
+        """
+        self._phase_estimation_scale = phase_estimation_scale
+        self._id_coefficient = id_coefficient
+
+        super().__init__(num_evaluation_qubits, circuit_result, phases)
+
+    # pylint: disable=arguments-differ
+    def filter_phases(self, cutoff: float = 0.0, scaled: bool = True,  # type: ignore
+                      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 = super().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_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 = super().most_likely_phase
+        return self._phase_estimation_scale.scale_phase(phase, self._id_coefficient)