Add quantum phase estimation template

.github/CHANGELOG.md

@@ -2,6 +2,53 @@
 
 <h3>New features since last release</h3>
 
+- Added the `QuantumPhaseEstimation` template for performing quantum phase estimation for an input
+  unitary matrix.
+  [(#1095)](https://github.com/PennyLaneAI/pennylane/pull/1095)
+
+  Consider the matrix corresponding to a rotation from an `RX` gate:
+
+  ```pycon
+  >>> phase = 5
+  >>> target_wires = [0]
+  >>> u = qml.RX(phase, wires=0).matrix
+  ```
+
+  The ``phase`` parameter can be estimated using ``QuantumPhaseEstimation``. For example, using five
+  phase-estimation qubits:
+
+  ```python
+  n_estimation_wires = 5
+  estimation_wires = range(1, n_estimation_wires + 1)
+
+  dev = qml.device("default.qubit", wires=n_estimation_wires + 1)
+
+  @qml.qnode(dev)
+  def circuit():
+      # Start in the |+> eigenstate of the unitary
+      qml.Hadamard(wires=target_wires)
+
+      QuantumPhaseEstimation(
+          u,
+          target_wires=target_wires,
+          estimation_wires=estimation_wires,
+      )
+
+      return qml.probs(estimation_wires)
+
+  phase_estimated = np.argmax(circuit()) / 2 ** n_estimation_wires
+
+  # Need to rescale phase due to convention of RX gate
+  phase_estimated = 4 * np.pi * (1 - phase)
+  ```
+
+  The resulting phase is a close approximation to the true value:
+
+  ```pycon
+  >>> phase_estimated
+  5.105088062083414
+  ```
+
 * Batches of shots can now be specified as a list, allowing measurement statistics
   to be course-grained with a single QNode evaluation.
   [(#1103)](https://github.com/PennyLaneAI/pennylane/pull/1103)
@@ -172,9 +219,7 @@
 
 This release contains contributions from (in alphabetical order):
 
-Thomas Bromley, Kyle Godbey, Josh Izaac, Daniel Polatajko, Chase Roberts, Maria Schuld.
-
-
+Thomas Bromley, Diego Guala, Kyle Godbey, Josh Izaac, Daniel Polatajko, Chase Roberts, Maria Schuld.
 
 # Release 0.14.1 (current release)
 

doc/_static/templates/subroutines/qpe.tex

@@ -0,0 +1,25 @@
+\documentclass{standalone}
+\usepackage{qcircuit}
+
+\begin{document}
+\Qcircuit @C=1em @R=0.5em {
+& & & & & & & & &\\
+& & & & \hspace{-0.5cm}\mbox{target wires} & &\\
+\\
+& \qw & \multigate{3}{U^{2^{n - 1}}} & \multigate{3}{U^{2^{n - 2}}} & \qw & \cdots & & \multigate{3}{U} & \qw & \qw  \\
+& \qw & \ghost{U^{2^{n - 1}}} & \ghost{U^{2^{n - 2}}} & \qw & \cdots & & \ghost{U} & \qw & \qw  \\
+& \qw & \ghost{U^{2^{n - 1}}} & \ghost{U^{2^{n - 2}}} & \qw & \cdots & & \ghost{U} & \qw & \qw  \\
+& \qw & \ghost{U^{2^{n - 1}}} & \ghost{U^{2^{n - 2}}} & \qw & \cdots & & \ghost{U} & \qw & \qw \\
+& & & & & & & & &\\
+& & & & & & & & &\\
+& & & & & & & & &\\
+& \gate{H} & \ctrl{-4} & \qw & \qw & \cdots & & \qw & \multigate{4}{QFT^{-1}} & \qw \\
+& \gate{H} & \qw & \ctrl{-5} & \qw & \cdots & & \qw & \ghost{QFT^{-1}} & \qw\\
+& \vdots & & & & \ddots & & & \\
+& & & & & & & &\\
+& \gate{H} \gategroup{10}{1}{17}{10}{.8em}{--}\gategroup{1}{1}{8}{10}{.8em}{--} & \qw & \qw & \qw & \cdots & & \ctrl{-8} & \ghost{QFT^{-1}} & \qw
+ \\
+ & & & \hspace{1.6cm}\mbox{estimation wires} & & & & & & \\
+  & & & & & & & & & \\
+}
+\end{document}

doc/introduction/templates.rst

@@ -193,6 +193,11 @@ of other templates.
   :description: Permute
   :figure: ../_static/templates/subroutines/permute.png
 
+.. customgalleryitem::
+  :link: ../code/api/pennylane.templates.subroutines.QuantumPhaseEstimation.html
+  :description: QuantumPhaseEstimation
+  :figure: ../_static/templates/subroutines/qpe.svg
+
 .. raw:: html
 
         <div style='clear:both'></div>

pennylane/templates/subroutines/__init__.py

@@ -23,3 +23,4 @@
 from .uccsd import UCCSD
 from .approx_time_evolution import ApproxTimeEvolution
 from .permute import Permute
+from .qpe import QuantumPhaseEstimation

pennylane/templates/subroutines/qpe.py

@@ -0,0 +1,123 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Contains the ``QuantumPhaseEstimation`` template.
+"""
+import pennylane as qml
+from pennylane.templates.decorator import template
+from pennylane.wires import Wires
+
+
+@template
+def QuantumPhaseEstimation(unitary, target_wires, estimation_wires):
+    r"""Performs the
+    `quantum phase estimation <https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm>`__
+    circuit.
+
+    Given a unitary matrix :math:`U`, this template applies the circuit for quantum phase
+    estimation. The unitary is applied to the qubits specified by ``target_wires`` and :math:`n`
+    qubits are used for phase estimation as specified by ``estimation_wires``.
+
+    .. figure:: ../../_static/templates/subroutines/qpe.svg
+        :align: center
+        :width: 60%
+        :target: javascript:void(0);
+
+    This circuit can be used to perform the standard quantum phase estimation algorithm, consisting
+    of the following steps:
+
+    #. Prepare ``target_wires`` in an eigenstate of :math:`U`. If that eigenstate has a
+       corresponding eigenvalue :math:`e^{2 \pi i \theta}` with phase :math:`\theta \in [0, 1)`,
+       this algorithm will measure :math:`\theta`.
+    #. Apply the ``QuantumPhaseEstimation`` circuit.
+    #. Measure ``estimation_wires`` using :func:`~.probs`, giving a probability distribution over
+       measurement outcomes in the computational basis.
+    #. Find the index of the largest value in the probability distribution and divide that number by
+       :math:`2^{n}`. This number will be an estimate of :math:`\theta` with an error that decreases
+       exponentially with the number of qubits :math:`n`.
+
+    Note that if :math:`\theta \in (-1, 0]`, we can estimate the phase by again finding the index
+    :math:`i` found in step 4 and calculating :math:`\theta \approx \frac{1 - i}{2^{n}}`. The
+    usage details below give an example of this case.
+
+    Args:
+        unitary (array): the phase estimation unitary, specified as a matrix
+        target_wires (Union[Wires, Sequence[int], or int]): the target wires to apply the unitary
+        estimation_wires (Union[Wires, Sequence[int], or int]): the wires to be used for phase
+            estimation
+
+    Raises:
+        QuantumFunctionError: if the ``target_wires`` and ``estimation_wires`` share a common
+            element
+
+    .. UsageDetails::
+
+        Consider the matrix corresponding to a rotation from an :class:`~.RX` gate:
+
+        .. code-block:: python
+
+            import pennylane as qml
+            from pennylane.templates import QuantumPhaseEstimation
+            from pennylane import numpy as np
+
+            phase = 5
+            target_wires = [0]
+            u = qml.RX(phase, wires=0).matrix
+
+        The ``phase`` parameter can be estimated using ``QuantumPhaseEstimation``. An example is
+        shown below using a register of five phase-estimation qubits:
+
+        .. code-block:: python
+
+            n_estimation_wires = 5
+            estimation_wires = range(1, n_estimation_wires + 1)
+
+            dev = qml.device("default.qubit", wires=n_estimation_wires + 1)
+
+            @qml.qnode(dev)
+            def circuit():
+                # Start in the |+> eigenstate of the unitary
+                qml.Hadamard(wires=target_wires)
+
+                QuantumPhaseEstimation(
+                    u,
+                    target_wires=target_wires,
+                    estimation_wires=estimation_wires,
+                )
+
+                return qml.probs(estimation_wires)
+
+            phase_estimated = np.argmax(circuit()) / 2 ** n_estimation_wires
+
+            # Need to rescale phase due to convention of RX gate
+            phase_estimated = 4 * np.pi * (1 - phase)
+    """
+
+    target_wires = Wires(target_wires)
+    estimation_wires = Wires(estimation_wires)
+
+    if len(Wires.shared_wires([target_wires, estimation_wires])) != 0:
+        raise qml.QuantumFunctionError("The target wires and estimation wires must be different")
+
+    unitary_powers = [unitary]
+
+    for _ in range(len(estimation_wires) - 1):
+        new_power = unitary_powers[-1] @ unitary_powers[-1]
+        unitary_powers.append(new_power)
+
+    for wire in estimation_wires:
+        qml.Hadamard(wire)
+        qml.ControlledQubitUnitary(unitary_powers.pop(), control_wires=wire, wires=target_wires)
+
+    qml.QFT(wires=estimation_wires).inv()

tests/templates/test_integration.py

@@ -155,6 +155,14 @@
                                {'weights': [[[ 0.17586701, -0.20382066]]]},
                                {'wires': [0, 1], 'init_state': np.array([1, 0], requires_grad=False)},
                                2),
+                              (qml.templates.QuantumPhaseEstimation,
+                               {},
+                               {
+                                   "unitary": np.array([[0, 1], [1, 0]]),
+                                   "target_wires": [0],
+                                   "estimation_wires": [1, 2],
+                               },
+                               3),
                               ]
 
 CV_DIFFABLE_NONDIFFABLE = [(qml.templates.DisplacementEmbedding,
@@ -706,6 +714,9 @@ def circuit_consec():
         if template.__name__ == 'UCCSD':
              kwargs2['s_wires'] = [nonconsecutive_wires[:3], nonconsecutive_wires[1:]]
              kwargs2['d_wires'] = [[nonconsecutive_wires[:2], nonconsecutive_wires[2:]]]
+        if template.__name__ == "QuantumPhaseEstimation":
+            kwargs2["target_wires"] = [nonconsecutive_wires[0]]
+            kwargs2["estimation_wires"] = nonconsecutive_wires[1:3]
 
         dev_nonconsec = qml.device('default.qubit', wires=nonconsecutive_wires)
 

tests/templates/test_subroutines.py

@@ -30,6 +30,7 @@
     UCCSD,
     ApproxTimeEvolution,
     Permute,
+    QuantumPhaseEstimation,
 )
 
 from pennylane.templates.subroutines.arbitrary_unitary import (
@@ -1477,3 +1478,75 @@ def test_subset_permutations_tape(
         # Make sure to start comparison after the set of RZs have been applied
         assert all(op.name == "SWAP" for op in tape.operations[len(wire_labels) :])
         assert [op.wires.labels for op in tape.operations[len(wire_labels) :]] == expected_wires
+
+
+class TestQuantumPhaseEstimation:
+    """Tests for the QuantumPhaseEstimation template from the pennylane.templates.subroutine
+    module."""
+
+    def test_same_wires(self):
+        """Tests if a QuantumFunctionError is raised if target_wires and estimation_wires contain a
+        common element"""
+
+        with pytest.raises(qml.QuantumFunctionError, match="The target wires and estimation wires"):
+            QuantumPhaseEstimation(np.eye(2), target_wires=[0, 1], estimation_wires=[1, 2])
+
+    def test_expected_tape(self):
+        """Tests if QuantumPhaseEstimation populates the tape as expected for a fixed example"""
+
+        m = qml.RX(0.3, wires=0).matrix
+
+        with qml.tape.QuantumTape() as tape:
+            QuantumPhaseEstimation(m, target_wires=[0], estimation_wires=[1, 2])
+
+        with qml.tape.QuantumTape() as tape2:
+            qml.Hadamard(1),
+            qml.ControlledQubitUnitary(m @ m, control_wires=[1], wires=[0]),
+            qml.Hadamard(2),
+            qml.ControlledQubitUnitary(m, control_wires=[2], wires=[0]),
+            qml.QFT(wires=[1, 2]).inv()
+
+        assert len(tape2.queue) == len(tape.queue)
+        assert all([op1.name == op2.name for op1, op2 in zip(tape.queue, tape2.queue)])
+        assert all([op1.wires == op2.wires for op1, op2 in zip(tape.queue, tape2.queue)])
+        assert np.allclose(tape.queue[1].matrix, tape2.queue[1].matrix)
+        assert np.allclose(tape.queue[3].matrix, tape2.queue[3].matrix)
+
+    def test_phase_estimated(self):
+        """Tests that the QPE circuit can correctly estimate the phase of a simple RX rotation."""
+        phase = 6
+        estimates = []
+        wire_range = range(2, 10)
+
+        for wires in wire_range:
+            dev = qml.device("default.qubit", wires=wires)
+            m = qml.RX(phase, wires=0).matrix
+            target_wires = [0]
+            estimation_wires = range(1, wires)
+
+            with qml.tape.QuantumTape() as tape:
+                # We want to prepare ourselves in an eigenstate of RX, in this case |+>
+                qml.Hadamard(wires=target_wires)
+
+                qml.templates.QuantumPhaseEstimation(
+                    m, target_wires=target_wires, estimation_wires=estimation_wires
+                )
+                qml.probs(estimation_wires)
+
+            res = tape.execute(dev).flatten()
+            initial_estimate = np.argmax(res) / 2 ** (wires - 1)
+
+            # We need to rescale because RX is exp(- i theta X / 2) and we expect a unitary of the
+            # form exp(2 pi i theta X)
+            rescaled_estimate = (1 - initial_estimate) * np.pi * 4
+            estimates.append(rescaled_estimate)
+
+        # Check that the error is monotonically decreasing
+        for i in range(len(estimates) - 1):
+            err1 = np.abs(estimates[i] - phase)
+            err2 = np.abs(estimates[i + 1] - phase)
+            assert err1 >= err2
+
+        # This is quite a large error, but we'd need to push the qubit number up more to get it
+        # lower
+        assert np.allclose(estimates[-1], phase, rtol=1e-2)