qml.QDrift template (#4671)

Note: I have left in comments some doubts 

## Context:

Third of 3 PRs adding the new QDrift template to allow for advanced
Trotter methods in Pennylane

## Description of the Change:

Implement the template

## Quick test to show that the solution is correct

It has been checked if at infinity (n = 10000) it converges to the exact
value.

```python
coeffs = [-1.1, .5]
ops = [qml.PauliX(0)@qml.PauliY(1), qml.PauliZ(1)]

H = qml.dot(coeffs, ops)
H1 = qml.Hamiltonian(coeffs, ops)

n = 10000
time = 0.1

qdrift_matrix1 = qml.matrix(qml.QDrift(H, time, n))
qdrift_matrix2 = qml.matrix(qml.QDrift(H1, time, n))
trotter_matrix = qml.matrix(qml.exp(H1, 1j * time))

print(np.round(qdrift_matrix1, 2))
print(np.round(qdrift_matrix2, 2))
print(np.round(trotter_matrix, 2))
```

Running it shows that the three matrices are equal

---------

Co-authored-by: soranjh <soran.jahangiri@gmail.com>
Co-authored-by: Jay Soni <jbsoni@uwaterloo.ca>

doc/introduction/templates.rst

@@ -227,6 +227,10 @@ Other useful templates which do not belong to the previous categories can be fou
   :description: :doc:`ApproxTimeEvolution <../code/api/pennylane.ApproxTimeEvolution>`
   :figure: _static/templates/subroutines/approx_time_evolution.png
 
+.. gallery-item::
+  :description: :doc:`QDrift <../code/api/pennylane.QDrift>`
+  :figure: _static/templates/subroutines/qdrift.png
+
 .. gallery-item::
   :description: :doc:`TrotterProduct <../code/api/pennylane.TrotterProduct>`
   :figure: _static/templates/subroutines/trotter_product.png

doc/releases/changelog-dev.md

@@ -4,6 +4,32 @@
 
 <h3>New features since last release</h3>
 
+* Approximate metrix exponentiation with random product formulas is available with the new `QDrift`
+  operation.
+  [(#4671)](https://github.com/PennyLaneAI/pennylane/pull/4671)
+
+  The product terms in `QDrift` are generated by randomly sampling from the Hamiltonian terms with
+  a probability that depends on the Hamiltonian coefficients.
+
+  ```python
+  coeffs = [0.25, 0.75]
+  ops = [qml.PauliX(0), qml.PauliZ(0)]
+  H = qml.dot(coeffs, ops)
+
+  dev = qml.device("default.qubit", wires=2)
+
+  @qml.qnode(dev)
+  def circuit():
+      qml.Hadamard(0)
+      qml.QDrift(H, time=1.2, n = 10)
+      return qml.probs()
+  ```
+
+  ```pycon
+  >>> circuit()
+  array([0.61814334, 0.        , 0.38185666, 0.        ])
+  ```
+
 <h4>Exponentiate Hamiltonians with flexible Trotter products 🐖</h4>
 
 * Higher-order Trotter-Suzuki methods are now easily accessible through a new operation

pennylane/templates/subroutines/__init__.py

@@ -35,4 +35,5 @@
 from .basis_rotation import BasisRotation
 from .qsvt import QSVT, qsvt
 from .select import Select
+from .qdrift import QDrift
 from .trotter import TrotterProduct

pennylane/templates/subroutines/qdrift.py

@@ -0,0 +1,289 @@
+# 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 template for QDrift subroutine."""
+
+import pennylane as qml
+from pennylane.operation import Operation
+from pennylane.math import requires_grad, unwrap
+from pennylane.ops import Sum, SProd, Hamiltonian
+
+
+@qml.QueuingManager.stop_recording()
+def _sample_decomposition(coeffs, ops, time, n=1, seed=None):
+    """Generate the randomly sampled decomposition
+
+    Args:
+        coeffs (array): the coefficients of the operations from each term in the Hamiltonian
+        ops (list[~.Operator]): the normalized operations from each term in the Hamiltonian
+        time (float): time to evolve under the target Hamiltonian
+        n (int): number of samples in the product, defaults to 1
+        seed (int): random seed. defaults to None
+
+    Returns:
+        list[~.Operator]: the decomposition of operations as per the approximation
+    """
+    normalization_factor = qml.math.sum(qml.math.abs(coeffs))
+    probs = qml.math.abs(coeffs) / normalization_factor
+    exps = [
+        qml.exp(base, (coeff / qml.math.abs(coeff)) * normalization_factor * time * 1j / n)
+        for base, coeff in zip(ops, coeffs)
+    ]
+
+    choice_rng = qml.math.random.default_rng(seed)
+    return choice_rng.choice(exps, p=probs, size=n, replace=True)
+
+
+class QDrift(Operation):
+    r"""An operation representing the QDrift approximation for the complex matrix exponential
+    of a given Hamiltonian.
+
+    The QDrift subroutine provides a method to approximate the matrix exponential of a Hamiltonian
+    expressed as a linear combination of terms which in general do not commute. For the Hamiltonian
+    :math:`H = \Sigma_j h_j H_{j}`, the product formula is constructed by random sampling from the
+    terms of the Hamiltonian with the probability :math:`p_j = h_j / \sum_{j} hj` as:
+
+    .. math::
+
+        \prod_{j}^{n} e^{i \lambda H_j \tau / n},
+
+    where :math:`\tau` is time, :math:`\lambda = \sum_j |h_j|` and :math:`n` is the total number of
+    terms to be sampled and added to the product. Note, the terms :math:`H_{j}` are assumed to be
+    normalized such that the "impact" of each term is fully encoded in the magnitude of :math:`h_{j}`.
+
+    The number of samples :math:`n` required for a given error threshold can be approximated by:
+
+    .. math::
+
+        n \ \approx \ \frac{2\lambda^{2}t^{2}}{\epsilon}
+
+    For more details see `Phys. Rev. Lett. 123, 070503 (2019) <https://arxiv.org/abs/1811.08017>`_.
+
+    Args:
+        hamiltonian (Union[.Hamiltonian, .Sum]): The Hamiltonian written as a sum of operations
+        time (float): The time of evolution, namely the parameter :math:`t` in :math:`e^{iHt}`
+        n (int): An integer representing the number of exponentiated terms
+        seed (int): The seed for the random number generator
+
+    Raises:
+        TypeError: The ``hamiltonian`` is not of type :class:`~.Hamiltonian`, or :class:`~.Sum`
+        QuantumFunctionError: If the coefficients of ``hamiltonian`` are trainable and are used
+            in a differentiable workflow.
+
+    **Example**
+
+    .. code-block:: python3
+
+        coeffs = [0.25, 0.75]
+        ops = [qml.PauliX(0), qml.PauliZ(0)]
+        H = qml.dot(coeffs, ops)
+
+        dev = qml.device("default.qubit", wires=2)
+        @qml.qnode(dev)
+        def my_circ():
+            # Prepare some state
+            qml.Hadamard(0)
+
+            # Evolve according to H
+            qml.QDrift(H, time=1.2, n=10, seed=10)
+
+            # Measure some quantity
+            return qml.probs()
+
+    >>> my_circ()
+    array([0.65379493, 0.        , 0.34620507, 0.        ])
+
+
+    .. details::
+        :title: Usage Details
+
+        We currently **Do NOT** support computing gradients with respect to the
+        coefficients of the input Hamiltonian. We can however compute the gradient
+        with respect to the evolution time:
+
+        .. code-block:: python3
+
+            dev = qml.device("default.qubit", wires=2)
+
+            @qml.qnode(dev)
+            def my_circ(time):
+                # Prepare H:
+                H = qml.dot([0.2, -0.1], [qml.PauliY(0), qml.PauliZ(1)])
+
+                # Prepare some state
+                qml.Hadamard(0)
+
+                # Evolve according to H
+                qml.QDrift(H, time, n=10, seed=10)
+
+                # Measure some quantity
+                return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))
+
+
+        >>> time = np.array(1.23)
+        >>> print(qml.grad(my_circ)(time))
+        0.27980654844422853
+
+        The error in the approximation of time evolution with the QDrift protocol is
+        directly related to the number of samples used in the product. We provide a
+        method to upper-bound the error:
+
+        >>> H = qml.dot([0.25, 0.75], [qml.PauliX(0), qml.PauliZ(0)])
+        >>> print(qml.QDrift.error(H, time=1.2, n=10))
+        0.3661197552925645
+
+    """
+
+    def __init__(  # pylint: disable=too-many-arguments
+        self, hamiltonian, time, n=1, seed=None, decomposition=None, id=None
+    ):
+        r"""Initialize the QDrift class"""
+
+        if isinstance(hamiltonian, Hamiltonian):
+            coeffs, ops = hamiltonian.terms()
+
+        elif isinstance(hamiltonian, Sum):
+            coeffs, ops = [], []
+            for op in hamiltonian:
+                try:
+                    coeffs.append(op.scalar)
+                    ops.append(op.base)
+                except AttributeError:  # coefficient is 1.0
+                    coeffs.append(1.0)
+                    ops.append(op)
+
+        else:
+            raise TypeError(
+                f"The given operator must be a PennyLane ~.Hamiltonian or ~.Sum got {hamiltonian}"
+            )
+
+        if len(ops) < 2:
+            raise ValueError(
+                "There should be atleast 2 terms in the Hamiltonian. Otherwise use `qml.exp`"
+            )
+
+        if any(requires_grad(coeff) for coeff in coeffs):
+            raise qml.QuantumFunctionError(
+                "The QDrift template currently doesn't support differentiation through the "
+                "coefficients of the input Hamiltonian."
+            )
+
+        if decomposition is None:  # need to do this to allow flatten and _unflatten
+            unwrapped_coeffs = unwrap(coeffs)
+            decomposition = _sample_decomposition(unwrapped_coeffs, ops, time, n=n, seed=seed)
+
+        self._hyperparameters = {
+            "n": n,
+            "seed": seed,
+            "base": hamiltonian,
+            "decomposition": decomposition,
+        }
+        super().__init__(time, wires=hamiltonian.wires, id=id)
+
+    @classmethod
+    def _unflatten(cls, data, metadata):
+        """Recreate an operation from its serialized format.
+
+        Args:
+            data: the trainable component of the operation
+            metadata: the non-trainable component of the operation
+
+        The output of ``Operator._flatten`` and the class type must be sufficient to reconstruct the original
+        operation with ``Operator._unflatten``.
+
+        **Example:**
+
+        >>> op = qml.Rot(1.2, 2.3, 3.4, wires=0)
+        >>> op._flatten()
+        ((1.2, 2.3, 3.4), (<Wires = [0]>, ()))
+        >>> qml.Rot._unflatten(*op._flatten())
+        >>> op = qml.PauliRot(1.2, "XY", wires=(0,1))
+        >>> op._flatten()
+        ((1.2,), (<Wires = [0, 1]>, (('pauli_word', 'XY'),)))
+        >>> op = qml.ctrl(qml.U2(3.4, 4.5, wires="a"), ("b", "c") )
+        >>> type(op)._unflatten(*op._flatten())
+        Controlled(U2(3.4, 4.5, wires=['a']), control_wires=['b', 'c'])
+
+        """
+        hyperparameters_dict = dict(metadata[1])
+        hamiltonian = hyperparameters_dict.pop("base")
+        return cls(hamiltonian, *data, **hyperparameters_dict)
+
+    @staticmethod
+    def compute_decomposition(*args, **kwargs):  # pylint: disable=unused-argument
+        r"""Representation of the operator as a product of other operators (static method).
+
+        .. math:: O = O_1 O_2 \dots O_n.
+
+        .. note::
+
+            Operations making up the decomposition should be queued within the
+            ``compute_decomposition`` method.
+
+        .. seealso:: :meth:`~.Operator.decomposition`.
+
+        Args:
+            *params (list): trainable parameters of the operator, as stored in the ``parameters`` attribute
+            wires (Iterable[Any], Wires): wires that the operator acts on
+            **hyperparams (dict): non-trainable hyperparameters of the operator, as stored in the ``hyperparameters`` attribute
+
+        Returns:
+            list[Operator]: decomposition of the operator
+        """
+        decomp = kwargs["decomposition"]
+
+        if qml.QueuingManager.recording():
+            for op in decomp:
+                qml.apply(op)
+
+        return decomp
+
+    @staticmethod
+    def error(hamiltonian, time, n=1):
+        r"""A method for determining the upper-bound for the error in the approximation of
+        the true matrix exponential.
+
+        The error is bounded according to the following expression:
+
+        .. math::
+
+            \epsilon \ \leq \ \frac{2\lambda^{2}t^{2}}{n}  e^{\frac{2 \lambda t}{n}},
+
+        where :math:`t` is time, :math:`\lambda = \sum_j |h_j|` and :math:`n` is the total number of
+        terms to be added to the product. For more details see `Phys. Rev. Lett. 123, 070503 (2019) <https://arxiv.org/abs/1811.08017>`_.
+
+        Args:
+        hamiltonian (Union[.Hamiltonian, .Sum]): The Hamiltonian written as a sum of operations
+        time (float): The time of evolution, namely the parameter :math:`t` in :math:`e^{-iHt}`
+        n (int): An integer representing the number of exponentiated terms. default is 1
+
+        Raises:
+            TypeError: The given operator must be a PennyLane .Hamiltonian or .Sum
+
+        Returns:
+            float: upper bound on the precision achievable using the QDrift protocol
+        """
+        if isinstance(hamiltonian, Hamiltonian):
+            lmbda = qml.math.sum(qml.math.abs(hamiltonian.coeffs))
+
+        elif isinstance(hamiltonian, Sum):
+            lmbda = qml.math.sum(
+                qml.math.abs(op.scalar) if isinstance(op, SProd) else 1.0 for op in hamiltonian
+            )
+
+        else:
+            raise TypeError(
+                f"The given operator must be a PennyLane ~.Hamiltonian or ~.Sum got {hamiltonian}"
+            )
+
+        return (2 * lmbda**2 * time**2 / n) * qml.math.exp(2 * lmbda * time / n)