diff --git a/doc/_static/templates/subroutines/qdrift.png b/doc/_static/templates/subroutines/qdrift.png
new file mode 100644
index 00000000000..79bb5a0e35b
Binary files /dev/null and b/doc/_static/templates/subroutines/qdrift.png differ
diff --git a/doc/introduction/templates.rst b/doc/introduction/templates.rst
index 26ee024653c..f5519e9fd22 100644
--- a/doc/introduction/templates.rst
+++ b/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
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 0b139fee1d2..9938017cd61 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -4,6 +4,32 @@
New features since last release
+* 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. ])
+ ```
+
Exponentiate Hamiltonians with flexible Trotter products 🐖
* Higher-order Trotter-Suzuki methods are now easily accessible through a new operation
diff --git a/pennylane/templates/subroutines/__init__.py b/pennylane/templates/subroutines/__init__.py
index 70c64d46fc6..df6f9a4eae2 100644
--- a/pennylane/templates/subroutines/__init__.py
+++ b/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
diff --git a/pennylane/templates/subroutines/qdrift.py b/pennylane/templates/subroutines/qdrift.py
new file mode 100644
index 00000000000..724d4f14604
--- /dev/null
+++ b/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) `_.
+
+ 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), (, ()))
+ >>> qml.Rot._unflatten(*op._flatten())
+ >>> op = qml.PauliRot(1.2, "XY", wires=(0,1))
+ >>> op._flatten()
+ ((1.2,), (, (('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) `_.
+
+ 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)
diff --git a/tests/templates/test_subroutines/test_qdrift.py b/tests/templates/test_subroutines/test_qdrift.py
new file mode 100644
index 00000000000..7a88a9fe99b
--- /dev/null
+++ b/tests/templates/test_subroutines/test_qdrift.py
@@ -0,0 +1,649 @@
+# Copyright 2018-2023 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.
+"""
+Tests for the QDrift template.
+"""
+import copy
+from functools import reduce
+
+import pytest
+
+import pennylane as qml
+from pennylane import numpy as qnp
+from pennylane.math import allclose, isclose
+from pennylane.templates.subroutines.qdrift import _sample_decomposition
+
+
+test_hamiltonians = (
+ (
+ [1, 1, 1],
+ [qml.PauliX(0), qml.PauliY(0), qml.PauliZ(1)],
+ ),
+ (
+ [1.23, -0.45j],
+ [qml.s_prod(0.1, qml.PauliX(0)), qml.prod(qml.PauliX(0), qml.PauliZ(1))],
+ ), # op arith
+ (
+ [1, -0.5, 0.5],
+ [qml.Identity(wires=[0, 1]), qml.PauliZ(0), qml.PauliZ(1)],
+ ),
+)
+
+
+class TestInitialization:
+ """Test that the class is intialized correctly."""
+
+ @pytest.mark.parametrize("n", (1, 2, 3))
+ @pytest.mark.parametrize("time", (0.5, 1, 2))
+ @pytest.mark.parametrize("seed", (None, 1234, 42))
+ @pytest.mark.parametrize("coeffs, ops", test_hamiltonians)
+ def test_init_correctly(self, coeffs, ops, time, n, seed): # pylint: disable=too-many-arguments
+ """Test that all of the attributes are initalized correctly."""
+ h = qml.dot(coeffs, ops)
+ op = qml.QDrift(h, time, n=n, seed=seed)
+
+ assert op.wires == h.wires
+ assert op.parameters == [time]
+ assert op.data == (time,)
+
+ assert op.hyperparameters["n"] == n
+ assert op.hyperparameters["seed"] == seed
+ assert op.hyperparameters["base"] == h
+
+ for term in op.hyperparameters["decomposition"]:
+ # the decomposition is solely made up of exponentials of ops sampled from hamiltonian terms
+ assert term.base in ops
+
+ def test_set_decomp(self):
+ """Test that setting the decomposition works correctly."""
+ h = qml.dot([1.23, -0.45], [qml.PauliX(0), qml.PauliY(0)])
+ decomposition = [
+ qml.exp(qml.PauliX(0), 0.5j * 1.68 / 3),
+ qml.exp(qml.PauliY(0), -0.5j * 1.68 / 3),
+ qml.exp(qml.PauliX(0), 0.5j * 1.68 / 3),
+ ]
+ op = qml.QDrift(h, 0.5, n=3, decomposition=decomposition)
+
+ assert op.hyperparameters["decomposition"] == decomposition
+
+ @pytest.mark.parametrize("n", (1, 2, 3))
+ @pytest.mark.parametrize("time", (0.5, 1, 2))
+ @pytest.mark.parametrize("seed", (None, 1234, 42))
+ @pytest.mark.parametrize("coeffs, ops", test_hamiltonians)
+ def test_copy(self, coeffs, ops, time, n, seed): # pylint: disable=too-many-arguments
+ """Test that we can make copies of QDrift correctly."""
+ h = qml.dot(coeffs, ops)
+ op = qml.QDrift(h, time, n=n, seed=seed)
+ new_op = copy.copy(op)
+
+ assert op.wires == new_op.wires
+ assert op.parameters == new_op.parameters
+ assert op.data == new_op.data
+ assert op.hyperparameters == new_op.hyperparameters
+ assert op is not new_op
+
+ @pytest.mark.parametrize(
+ "hamiltonian, raise_error",
+ (
+ (qml.PauliX(0), True),
+ (qml.prod(qml.PauliX(0), qml.PauliZ(1)), True),
+ (qml.Hamiltonian([1.23, 3.45], [qml.PauliX(0), qml.PauliZ(1)]), False),
+ (qml.dot([1.23, 3.45], [qml.PauliX(0), qml.PauliZ(1)]), False),
+ ),
+ )
+ def test_error_type(self, hamiltonian, raise_error):
+ """Test an error is raised of an incorrect type is passed"""
+ if raise_error:
+ with pytest.raises(
+ TypeError, match="The given operator must be a PennyLane ~.Hamiltonian or ~.Sum"
+ ):
+ qml.QDrift(hamiltonian, time=1.23)
+ else:
+ try:
+ qml.QDrift(hamiltonian, time=1.23)
+ except TypeError:
+ assert False # test should fail if an error was raised when we expect it not to
+
+ def test_error_hamiltonian(self):
+ """Test that a hamiltonian must have atleast 2 terms to be supported."""
+ msg = "There should be atleast 2 terms in the Hamiltonian."
+ with pytest.raises(ValueError, match=msg):
+ h = qml.Hamiltonian([1.0], [qml.PauliX(0)])
+ qml.QDrift(h, 1.23, n=2, seed=None)
+
+ @pytest.mark.parametrize("coeffs, ops", test_hamiltonians)
+ def test_flatten_and_unflatten(self, coeffs, ops):
+ """Test that the flatten and unflatten methods work correctly."""
+ time, n, seed = (0.5, 2, 1234)
+ hamiltonian = qml.dot(coeffs, ops)
+ op = qml.QDrift(hamiltonian, time, n=n, seed=seed)
+ decomp = op.decomposition()
+
+ data, metadata = op._flatten() # pylint: disable=protected-access
+ assert data[0] == time
+ assert metadata[0] == op.wires
+ assert dict(metadata[1]) == {
+ "n": n,
+ "seed": seed,
+ "base": hamiltonian,
+ "decomposition": decomp,
+ }
+
+ new_op = type(op)._unflatten(data, metadata) # pylint: disable=protected-access
+ assert qml.equal(op, new_op)
+ assert new_op is not op
+
+
+class TestDecomposition:
+ """Test decompositions are generated correctly."""
+
+ @pytest.mark.parametrize("n", (1, 2, 3))
+ @pytest.mark.parametrize("time", (0.5, 1, 2))
+ @pytest.mark.parametrize("seed", (None, 1234, 42))
+ @pytest.mark.parametrize("coeffs, ops", test_hamiltonians)
+ def test_private_sample(self, coeffs, ops, time, seed, n): # pylint: disable=too-many-arguments
+ """Test the private function which samples the decomposition"""
+ ops_to_coeffs = dict(zip(ops, coeffs))
+ normalization = qnp.sum(qnp.abs(coeffs))
+
+ with qml.tape.QuantumTape() as tape:
+ decomp = _sample_decomposition(coeffs, ops, time, n, seed)
+
+ assert len(decomp) == n
+ assert len(tape.operations) == 0 # no queuing
+ for term in decomp:
+ coeff = ops_to_coeffs[term.base]
+ s = coeff / qml.math.abs(coeff)
+
+ assert term.base in ops # sample from ops
+ assert term.coeff == (s * normalization * time * 1j / n) # with this exponent
+
+ @pytest.mark.parametrize("coeffs", ([0.99, 0.01], [0.5 + 0.49j, -0.01j]))
+ def test_private_sample_statistics(self, coeffs):
+ """Test the private function samples from the right distribution"""
+ ops = [qml.PauliX(0), qml.PauliZ(1)]
+ decomp = _sample_decomposition(coeffs, ops, 1.23, n=10, seed=1234)
+
+ # High probability we only sample PauliX!
+ assert all(isinstance(op.base, qml.PauliX) for op in decomp)
+
+ @pytest.mark.parametrize("seed", (1234, 42))
+ def test_compute_decomposition(self, seed):
+ """Test that the decomposition is computed and queues correctly."""
+ coeffs = [1, -0.5, 0.5]
+ ops = [qml.Identity(wires=[0, 1]), qml.PauliZ(0), qml.PauliZ(1)]
+
+ h = qml.dot(coeffs, ops)
+ op = qml.QDrift(h, time=1.23, n=10, seed=seed)
+
+ expected_decomp = _sample_decomposition(coeffs, ops, 1.23, 10, seed=seed)
+
+ with qml.tape.QuantumTape() as tape:
+ decomp = op.compute_decomposition(*op.parameters, **op.hyperparameters)
+
+ assert all(decomp == tape.operations) # queue matches decomp with circuit ordering
+ assert all(decomp == expected_decomp) # sample the same ops
+
+
+class TestIntegration:
+ """Test that the QDrift template integrates well with the rest of PennyLane"""
+
+ @pytest.mark.parametrize("n", (1, 2, 3))
+ @pytest.mark.parametrize("time", (0.5, 1, 2))
+ @pytest.mark.parametrize("seed", (1234, 42))
+ @pytest.mark.parametrize("coeffs, ops", test_hamiltonians)
+ def test_execution(self, coeffs, ops, time, n, seed): # pylint: disable=too-many-arguments
+ """Test that the circuit executes as expected"""
+ hamiltonian = qml.dot(coeffs, ops)
+ wires = hamiltonian.wires
+ dev = qml.device("default.qubit", wires=wires)
+
+ @qml.qnode(dev)
+ def circ():
+ qml.QDrift(hamiltonian, time, n=n, seed=seed)
+ return qml.state()
+
+ expected_decomp = _sample_decomposition(coeffs, ops, time, n=n, seed=seed)
+
+ initial_state = qnp.zeros(2 ** (len(wires)))
+ initial_state[0] = 1
+
+ expected_state = (
+ reduce(
+ lambda x, y: x @ y,
+ [qml.matrix(op, wire_order=wires) for op in expected_decomp],
+ )
+ @ initial_state
+ )
+ state = circ()
+
+ assert allclose(expected_state, state)
+
+ @pytest.mark.autograd
+ @pytest.mark.parametrize("seed", (1234, 42))
+ @pytest.mark.parametrize("coeffs, ops", test_hamiltonians)
+ def test_execution_autograd(self, coeffs, ops, seed):
+ """Test that the circuit executes as expected using autograd"""
+
+ time = qnp.array(0.5)
+ coeffs = qnp.array(coeffs, requires_grad=False)
+
+ dev = qml.device("default.qubit", wires=[0, 1])
+
+ @qml.qnode(dev)
+ def circ(time):
+ hamiltonian = qml.dot(coeffs, ops)
+ qml.QDrift(hamiltonian, time, n=2, seed=seed)
+ return qml.state()
+
+ expected_decomp = _sample_decomposition(coeffs, ops, time, n=2, seed=seed)
+
+ initial_state = qnp.array([1.0, 0.0, 0.0, 0.0])
+
+ expected_state = (
+ reduce(
+ lambda x, y: x @ y,
+ [qml.matrix(op, wire_order=[0, 1]) for op in expected_decomp],
+ )
+ @ initial_state
+ )
+ state = circ(time)
+
+ assert allclose(expected_state, state)
+
+ @pytest.mark.torch
+ @pytest.mark.parametrize("seed", (1234, 42))
+ @pytest.mark.parametrize("coeffs, ops", test_hamiltonians)
+ def test_execution_torch(self, coeffs, ops, seed):
+ """Test that the circuit executes as expected using torch"""
+ import torch
+
+ time = torch.tensor(0.5, dtype=torch.complex128, requires_grad=True)
+ dev = qml.device("default.qubit", wires=[0, 1])
+
+ @qml.qnode(dev)
+ def circ(time):
+ hamiltonian = qml.dot(coeffs, ops)
+ qml.QDrift(hamiltonian, time, n=2, seed=seed)
+ return qml.state()
+
+ expected_decomp = _sample_decomposition(coeffs, ops, time, n=2, seed=seed)
+
+ initial_state = torch.tensor([1.0, 0.0, 0.0, 0.0], dtype=torch.complex128)
+
+ expected_state = (
+ reduce(
+ lambda x, y: x @ y,
+ [qml.matrix(op, wire_order=[0, 1]) for op in expected_decomp],
+ )
+ @ initial_state
+ )
+ state = circ(time)
+
+ assert allclose(expected_state, state)
+
+ @pytest.mark.tf
+ @pytest.mark.parametrize("seed", (1234, 42))
+ @pytest.mark.parametrize("coeffs, ops", test_hamiltonians)
+ def test_execution_tf(self, coeffs, ops, seed):
+ """Test that the circuit executes as expected using tensorflow"""
+ import tensorflow as tf
+
+ time = tf.Variable(0.5, dtype=tf.complex128)
+ dev = qml.device("default.qubit", wires=[0, 1])
+
+ @qml.qnode(dev)
+ def circ(time):
+ hamiltonian = qml.dot(coeffs, ops)
+ qml.QDrift(hamiltonian, time, n=2, seed=seed)
+ return qml.state()
+
+ expected_decomp = _sample_decomposition(coeffs, ops, time, n=2, seed=seed)
+
+ initial_state = tf.Variable([1.0, 0.0, 0.0, 0.0], dtype=tf.complex128)
+
+ expected_state = tf.linalg.matvec(
+ reduce(
+ lambda x, y: x @ y,
+ [qml.matrix(op, wire_order=[0, 1]) for op in expected_decomp],
+ ),
+ initial_state,
+ )
+ state = circ(time)
+
+ assert allclose(expected_state, state)
+
+ @pytest.mark.jax
+ @pytest.mark.parametrize("seed", (1234, 42))
+ @pytest.mark.parametrize("coeffs, ops", test_hamiltonians)
+ def test_execution_jax(self, coeffs, ops, seed):
+ """Test that the circuit executes as expected using jax"""
+ from jax import numpy as jnp
+
+ time = jnp.array(0.5)
+ dev = qml.device("default.qubit", wires=[0, 1])
+
+ @qml.qnode(dev)
+ def circ(time):
+ hamiltonian = qml.dot(coeffs, ops)
+ qml.QDrift(hamiltonian, time, n=2, seed=seed)
+ return qml.state()
+
+ expected_decomp = _sample_decomposition(coeffs, ops, time, n=2, seed=seed)
+
+ initial_state = jnp.array([1.0, 0.0, 0.0, 0.0])
+
+ expected_state = (
+ reduce(
+ lambda x, y: x @ y,
+ [qml.matrix(op, wire_order=[0, 1]) for op in expected_decomp],
+ )
+ @ initial_state
+ )
+ state = circ(time)
+
+ assert allclose(expected_state, state)
+
+ @pytest.mark.jax
+ @pytest.mark.parametrize("seed", (1234, 42))
+ @pytest.mark.parametrize("coeffs, ops", test_hamiltonians)
+ def test_execution_jaxjit(self, coeffs, ops, seed):
+ """Test that the circuit executes as expected using jax jit"""
+ import jax
+ from jax import numpy as jnp
+
+ time = jnp.array(0.5)
+ dev = qml.device("default.qubit", wires=[0, 1])
+
+ @jax.jit
+ @qml.qnode(dev, interface="jax")
+ def circ(time):
+ hamiltonian = qml.sum(*(qml.s_prod(coeff, op) for coeff, op in zip(coeffs, ops)))
+ qml.QDrift(hamiltonian, time, n=2, seed=seed)
+ return qml.state()
+
+ expected_decomp = _sample_decomposition(coeffs, ops, time, n=2, seed=seed)
+
+ initial_state = jnp.array([1.0, 0.0, 0.0, 0.0])
+
+ expected_state = (
+ reduce(
+ lambda x, y: x @ y,
+ [qml.matrix(op, wire_order=[0, 1]) for op in expected_decomp],
+ )
+ @ initial_state
+ )
+ state = circ(time)
+
+ assert allclose(expected_state, state)
+
+ @pytest.mark.autograd
+ def test_error_gradient_workflow_autograd(self):
+ """Test that an error is raised if we require a gradient of QDrift with respect to hamiltonian coefficients."""
+ time = qnp.array(1.5)
+ coeffs = qnp.array([1.23, -0.45], requires_grad=True)
+
+ terms = [qml.PauliX(0), qml.PauliZ(0)]
+ dev = qml.device("default.qubit", wires=1)
+
+ @qml.qnode(dev)
+ def circ(time, coeffs):
+ h = qml.dot(coeffs, terms)
+ qml.QDrift(h, time, n=3)
+ return qml.expval(qml.Hadamard(0))
+
+ msg = "The QDrift template currently doesn't support differentiation through the coefficients of the input Hamiltonian."
+ with pytest.raises(qml.QuantumFunctionError, match=msg):
+ qml.grad(circ)(time, coeffs)
+
+ @pytest.mark.torch
+ def test_error_gradient_workflow_torch(self):
+ """Test that an error is raised if we require a gradient of QDrift with respect to hamiltonian coefficients."""
+ import torch
+
+ time = torch.tensor(1.5, dtype=torch.complex128, requires_grad=True)
+ coeffs = torch.tensor([1.23, -0.45], dtype=torch.complex128, requires_grad=True)
+
+ terms = [qml.PauliX(0), qml.PauliZ(0)]
+ dev = qml.device("default.qubit", wires=1)
+
+ @qml.qnode(dev)
+ def circ(time, coeffs):
+ h = qml.dot(coeffs, terms)
+ qml.QDrift(h, time, n=3)
+ return qml.expval(qml.Hadamard(0))
+
+ msg = "The QDrift template currently doesn't support differentiation through the coefficients of the input Hamiltonian."
+ with pytest.raises(qml.QuantumFunctionError, match=msg):
+ res_circ = circ(time, coeffs)
+ res_circ.backward()
+
+ @pytest.mark.tf
+ def test_error_gradient_workflow_tf(self):
+ """Test that an error is raised if we require a gradient of QDrift with respect to hamiltonian coefficients."""
+ import tensorflow as tf
+
+ time = tf.Variable(1.5, dtype=tf.complex128)
+ coeffs = tf.Variable([1.23, -0.45], dtype=tf.complex128)
+
+ terms = [qml.PauliX(0), qml.PauliZ(0)]
+ dev = qml.device("default.qubit", wires=1)
+
+ @qml.qnode(dev)
+ def circ(time, coeffs):
+ h = qml.sum(
+ qml.s_prod(coeffs[0], terms[0]),
+ qml.s_prod(coeffs[1], terms[1]),
+ )
+ qml.QDrift(h, time, n=3)
+ return qml.expval(qml.Hadamard(0))
+
+ msg = "The QDrift template currently doesn't support differentiation through the coefficients of the input Hamiltonian."
+ with pytest.raises(qml.QuantumFunctionError, match=msg):
+ with tf.GradientTape() as tape:
+ result = circ(time, coeffs)
+ tape.gradient(result, coeffs)
+
+ @pytest.mark.jax
+ def test_error_gradient_workflow_jax(self):
+ """Test that an error is raised if we require a gradient of QDrift with respect to hamiltonian coefficients."""
+ import jax
+ from jax import numpy as jnp
+
+ time = jnp.array(1.5)
+ coeffs = jnp.array([1.23, -0.45])
+
+ terms = [qml.PauliX(0), qml.PauliZ(0)]
+ dev = qml.device("default.qubit", wires=1)
+
+ @qml.qnode(dev)
+ def circ(time, coeffs):
+ h = qml.dot(coeffs, terms)
+ qml.QDrift(h, time, n=3)
+ return qml.expval(qml.Hadamard(0))
+
+ msg = "The QDrift template currently doesn't support differentiation through the coefficients of the input Hamiltonian."
+ with pytest.raises(qml.QuantumFunctionError, match=msg):
+ jax.grad(circ, argnums=[1])(time, coeffs)
+
+ @pytest.mark.autograd
+ @pytest.mark.parametrize("n", (1, 5, 10))
+ @pytest.mark.parametrize("seed", (1234, 42))
+ def test_autograd_gradient(self, n, seed):
+ """Test that the gradient is computed correctly"""
+ time = qnp.array(1.5)
+ coeffs = qnp.array([1.23, -0.45], requires_grad=False)
+ terms = [qml.PauliX(0), qml.PauliZ(0)]
+
+ dev = qml.device("default.qubit", wires=1)
+
+ @qml.qnode(dev)
+ def circ(time, coeffs):
+ h = qml.dot(coeffs, terms)
+ qml.QDrift(h, time, n=n, seed=seed)
+ return qml.expval(qml.Hadamard(0))
+
+ @qml.qnode(dev)
+ def reference_circ(time, coeffs):
+ with qml.QueuingManager.stop_recording():
+ decomp = _sample_decomposition(coeffs, terms, time, n, seed)
+
+ for op in decomp:
+ qml.apply(op)
+
+ return qml.expval(qml.Hadamard(0))
+
+ measured_grad = qml.grad(circ)(time, coeffs)
+ reference_grad = qml.grad(reference_circ)(time, coeffs)
+ assert allclose(measured_grad, reference_grad)
+
+ @pytest.mark.torch
+ @pytest.mark.parametrize("n", (1, 5, 10))
+ @pytest.mark.parametrize("seed", (1234, 42))
+ def test_torch_gradient(self, n, seed):
+ """Test that the gradient is computed correctly using torch"""
+ import torch
+
+ time = torch.tensor(1.5, dtype=torch.complex128, requires_grad=True)
+ coeffs = torch.tensor([1.23, -0.45], dtype=torch.complex128, requires_grad=False)
+ ref_time = torch.tensor(1.5, dtype=torch.complex128, requires_grad=True)
+ ref_coeffs = torch.tensor([1.23, -0.45], dtype=torch.complex128, requires_grad=False)
+ terms = [qml.PauliX(0), qml.PauliZ(0)]
+
+ dev = qml.device("default.qubit", wires=1)
+
+ @qml.qnode(dev)
+ def circ(time, coeffs):
+ h = qml.dot(coeffs, terms)
+ qml.QDrift(h, time, n=n, seed=seed)
+ return qml.expval(qml.Hadamard(0))
+
+ @qml.qnode(dev)
+ def reference_circ(time, coeffs):
+ with qml.QueuingManager.stop_recording():
+ decomp = _sample_decomposition(coeffs, terms, time, n, seed)
+
+ for op in decomp:
+ qml.apply(op)
+
+ return qml.expval(qml.Hadamard(0))
+
+ res_circ = circ(time, coeffs)
+ res_circ.backward()
+ measured_grad = time.grad
+
+ ref_circ = reference_circ(ref_time, ref_coeffs)
+ ref_circ.backward()
+ reference_grad = ref_time.grad
+
+ assert allclose(measured_grad, reference_grad)
+
+ @pytest.mark.tf
+ @pytest.mark.parametrize("n", (1, 5, 10))
+ @pytest.mark.parametrize("seed", (1234, 42))
+ def test_tf_gradient(self, n, seed):
+ """Test that the gradient is computed correctly using tensorflow"""
+ import tensorflow as tf
+
+ time = tf.Variable(1.5, dtype=tf.complex128)
+ coeffs = [1.23, -0.45]
+ terms = [qml.PauliX(0), qml.PauliZ(0)]
+
+ dev = qml.device("default.qubit", wires=1)
+
+ @qml.qnode(dev)
+ def circ(time, coeffs):
+ h = qml.dot(coeffs, terms)
+ qml.QDrift(h, time, n=n, seed=seed)
+ return qml.expval(qml.Hadamard(0))
+
+ @qml.qnode(dev)
+ def reference_circ(time, coeffs):
+ with qml.QueuingManager.stop_recording():
+ decomp = _sample_decomposition(coeffs, terms, time, n, seed)
+
+ for op in decomp:
+ qml.apply(op)
+
+ return qml.expval(qml.Hadamard(0))
+
+ with tf.GradientTape() as tape:
+ result = circ(time, coeffs)
+ measured_grad = tape.gradient(result, time)
+
+ with tf.GradientTape() as tape:
+ result = reference_circ(time, coeffs)
+ reference_grad = tape.gradient(result, time)
+
+ assert allclose(measured_grad, reference_grad)
+
+ @pytest.mark.jax
+ @pytest.mark.parametrize("n", (1, 5, 10))
+ @pytest.mark.parametrize("seed", (1234, 42))
+ def test_jax_gradient(self, n, seed):
+ """Test that the gradient is computed correctly using jax"""
+ import jax
+ from jax import numpy as jnp
+
+ time = jnp.array(1.5)
+ coeffs = jnp.array([1.23, -0.45])
+ terms = [qml.PauliX(0), qml.PauliZ(0)]
+
+ dev = qml.device("default.qubit", wires=1)
+
+ @qml.qnode(dev)
+ def circ(time, coeffs):
+ h = qml.dot(coeffs, terms)
+ qml.QDrift(h, time, n=n, seed=seed)
+ return qml.expval(qml.Hadamard(0))
+
+ @qml.qnode(dev)
+ def reference_circ(time, coeffs):
+ with qml.QueuingManager.stop_recording():
+ decomp = _sample_decomposition(coeffs, terms, time, n, seed)
+
+ for op in decomp:
+ qml.apply(op)
+
+ return qml.expval(qml.Hadamard(0))
+
+ measured_grad = jax.grad(circ, argnums=[0])(time, coeffs)
+ reference_grad = jax.grad(reference_circ, argnums=[0])(time, coeffs)
+ assert allclose(measured_grad, reference_grad)
+
+
+test_error_data = ( # Computed by hand
+ (qml.dot([1.23, -0.45j], [qml.PauliX(0), qml.PauliZ(1)]), 0.5, 5, 0.3949494464),
+ (qml.Hamiltonian([1.23, -0.45], [qml.PauliX(0), qml.PauliZ(1)]), 0.5, 5, 0.3949494464),
+ (
+ qml.dot([1, -0.5, 0.5j], [qml.Identity(wires=[0, 1]), qml.PauliZ(0), qml.Hadamard(1)]),
+ 3,
+ 100,
+ 0.81179773314,
+ ),
+)
+
+
+@pytest.mark.parametrize("h, time, n, expected_error", test_error_data)
+def test_error_func(h, time, n, expected_error):
+ """Test that the error function computes the expected precision correctly"""
+ computed_error = qml.QDrift.error(h, time, n)
+ assert isclose(computed_error, expected_error)
+
+
+def test_error_func_type_error():
+ """Test that an error is raised if the wrong type is passed for hamiltonian"""
+ msg = "The given operator must be a PennyLane ~.Hamiltonian or ~.Sum"
+ with pytest.raises(TypeError, match=msg):
+ qml.QDrift.error(qml.PauliX(0), time=1.23, n=10)