Source code for pennylane.templates.subroutines.trotter
# 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.
"""
Contains templates for Suzuki-Trotter approximation based subroutines.
"""
import copy
from collections import defaultdict
from functools import wraps
import pennylane as qml
from pennylane.operation import Operation, Operator
from pennylane.ops import Sum
from pennylane.ops.op_math import SProd
from pennylane.resource import Resources, ResourcesOperation
from pennylane.resource.error import ErrorOperation, SpectralNormError
from pennylane.wires import Wires
def _scalar(order):
"""Compute the scalar used in the recursive expression.
Args:
order (int): order of Trotter product (assume order is an even integer > 2).
Returns:
float: scalar to be used in the recursive expression.
"""
root = 1 / (order - 1)
return (4 - 4**root) ** -1
@qml.QueuingManager.stop_recording()
def _recursive_expression(x, order, ops):
"""Generate a list of operations using the
recursive expression which defines the Trotter product.
Args:
x (complex): the evolution 'time'
order (int): the order of the Trotter expansion
ops (Iterable(~.Operators)): a list of terms in the Hamiltonian
Returns:
list: the approximation as product of exponentials of the Hamiltonian terms
"""
if order == 1:
return [qml.exp(op, x * 1j) for op in ops]
if order == 2:
return [qml.exp(op, x * 0.5j) for op in ops + ops[::-1]]
scalar_1 = _scalar(order)
scalar_2 = 1 - 4 * scalar_1
ops_lst_1 = _recursive_expression(scalar_1 * x, order - 2, ops)
ops_lst_2 = _recursive_expression(scalar_2 * x, order - 2, ops)
return (2 * ops_lst_1) + ops_lst_2 + (2 * ops_lst_1)
[docs]class TrotterProduct(ErrorOperation, ResourcesOperation):
r"""An operation representing the Suzuki-Trotter product approximation for the complex matrix
exponential of a given Hamiltonian.
The Suzuki-Trotter product formula provides a method to approximate the matrix exponential of
Hamiltonian expressed as a linear combination of operands which in general do not commute.
Consider the Hamiltonian :math:`H = \Sigma^{N}_{j=0} O_{j}`, the product formula is constructed
using symmetrized products of the terms in the Hamiltonian. The symmetrized products of order
:math:`m \in [1, 2, 4, ..., 2k]` with :math:`k \in \mathbb{N}` are given by:
.. math::
\begin{align}
S_{1}(t) &= \Pi_{j=0}^{N} \ e^{i t O_{j}} \\
S_{2}(t) &= \Pi_{j=0}^{N} \ e^{i \frac{t}{2} O_{j}} \cdot \Pi_{j=N}^{0} \ e^{i \frac{t}{2} O_{j}} \\
&\vdots \\
S_{m}(t) &= S_{m-2}(p_{m}t)^{2} \cdot S_{m-2}((1-4p_{m})t) \cdot S_{m-2}(p_{m}t)^{2},
\end{align}
where the coefficient is :math:`p_{m} = 1 / (4 - \sqrt[m - 1]{4})`. The :math:`m`th order,
:math:`n`-step Suzuki-Trotter approximation is then defined as:
.. math:: e^{iHt} \approx \left [S_{m}(t / n) \right ]^{n}.
For more details see `J. Math. Phys. 32, 400 (1991) <https://pubs.aip.org/aip/jmp/article-abstract/32/2/400/229229>`_.
Args:
hamiltonian (Union[.Hamiltonian, .Sum, .SProd]): The Hamiltonian written as a linear combination
of operators with known matrix exponentials.
time (float): The time of evolution, namely the parameter :math:`t` in :math:`e^{iHt}`
n (int): An integer representing the number of Trotter steps to perform
order (int): An integer (:math:`m`) representing the order of the approximation (must be 1 or even)
check_hermitian (bool): A flag to enable the validation check to ensure this is a valid unitary operator
Raises:
TypeError: The ``hamiltonian`` is not of type :class:`~.Sum`.
ValueError: The ``hamiltonian`` has only one term or no terms.
ValueError: One or more of the terms in ``hamiltonian`` are not Hermitian
(only for ``check_hermitian=True``)
ValueError: The ``order`` is not one or a positive even integer.
**Example**
.. code-block:: python3
coeffs = [0.25, 0.75]
ops = [qml.X(0), qml.Z(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.TrotterProduct(H, time=2.4, order=2)
# Measure some quantity
return qml.state()
>>> my_circ()
array([-0.13259524+0.59790098j, 0. +0.j , -0.13259524-0.77932754j, 0. +0.j ])
.. warning::
The Trotter-Suzuki decomposition depends on the order of the summed observables. Two
mathematically identical :class:`~.LinearCombination` objects may undergo different time
evolutions due to the order in which those observables are stored. The order of observables
can be queried using the :attr:`~.Sum.operands` attribute. Also see the advanced example
below.
.. warning::
``TrotterProduct`` does not automatically simplify the input Hamiltonian. This allows
for a more fine-grained control over the decomposition but also risks an increased
runtime and/or number of gates. Simplification can be performed manually by
applying :func:`~.simplify` to your Hamiltonian before using it in ``TrotterProduct``.
.. details::
:title: Usage Details
An *upper-bound* for the error in approximating time-evolution using this operator can be
computed by calling :func:`~.TrotterProduct.error()`. It is computed using two different
methods; the "one-norm-bound" scaling method and the "commutator-bound" scaling method.
(see `Childs et al. (2021) <https://arxiv.org/abs/1912.08854>`_)
>>> hamiltonian = qml.dot([1.0, 0.5, -0.25], [qml.X(0), qml.Y(0), qml.Z(0)])
>>> op = qml.TrotterProduct(hamiltonian, time=0.01, order=2)
>>> op.error(method="one-norm-bound")
SpectralNormError(8.039062500000003e-06)
>>> op.error(method="commutator-bound")
SpectralNormError(6.166666666666668e-06)
This operation is similar to the :class:`~.ApproxTimeEvolution`. One can recover the behaviour
of :class:`~.ApproxTimeEvolution` by taking the adjoint:
>>> qml.adjoint(qml.TrotterProduct(hamiltonian, time, order=1, n=n))
The grouping of terms in the ``operands`` attribute of the ``hamiltonian`` impacts
the structure of the gates created by ``TrotterProduct``. To understand this, first
consider this simple two-qubit Hamiltonian with four Pauli word terms:
>>> coeffs = [0.5, 0.2, 0.1, -0.6]
>>> ops = [qml.X(0), qml.Y(1), qml.Y(0) @ qml.Z(1), qml.X(0) @ qml.Y(1)]
>>> H_flat = qml.dot(coeffs, ops)
>>> H_flat
>>> print(*H_flat.operands, sep="\n")
0.5 * X(0)
0.2 * Y(1)
0.1 * (Y(0) @ Z(1))
-0.6 * (X(0) @ Y(1))
As we can see, each Pauli word contributes an individual operand. As a result, the
``TrotterProduct`` (of first order, for simplicity) of this Hamiltonian will contain four
exponentials per Trotter step:
>>> qml.TrotterProduct(H_flat, 1., n=1, order=1).decomposition()
[Exp(1j -0.6 * (X(0) @ Y(1))),
Exp(1j 0.1 * (Y(0) @ Z(1))),
Exp(1j 0.2 * Y(1)),
Exp(1j 0.5 * X(0))]
If we first create two operands with two Pauli words each and then sum those, this is
reflected in the structure of the operator:
>>> H_grouped = qml.sum(qml.dot(coeffs[:2], ops[:2]), qml.dot(coeffs[2:], ops[2:]))
>>> print(*H_grouped.operands, sep="\n")
0.5 * X(0) + 0.2 * Y(1)
0.1 * (Y(0) @ Z(1)) + -0.6 * (X(0) @ Y(1))
The ``TrotterProduct`` accordingly has a different structure as well:
>>> qml.TrotterProduct(H_grouped, 1., n=1, order=1).decomposition()
[Exp(1j 0.1 * (Y(0) @ Z(1)) + -0.6 * (X(0) @ Y(1))),
Exp(1j 0.5 * X(0) + 0.2 * Y(1))]
As we can see, the ``operands`` structure of the Hamiltonian directly impacts the
constructed Trotter circuit, and in general, those circuits will be different
approximations to the true time evolution.
We can also compute the gradient with respect to the coefficients of the Hamiltonian and the
evolution time:
.. code-block:: python3
@qml.qnode(dev)
def my_circ(c1, c2, time):
# Prepare H:
H = qml.dot([c1, c2], [qml.X(0), qml.Z(0)])
# Prepare some state
qml.Hadamard(0)
# Evolve according to H
qml.TrotterProduct(H, time, order=2)
# Measure some quantity
return qml.expval(qml.Z(0) @ qml.Z(1))
>>> args = np.array([1.23, 4.5, 0.1])
>>> qml.grad(my_circ)(*tuple(args))
(tensor(0.00961064, requires_grad=True), tensor(-0.12338274, requires_grad=True), tensor(-5.43401259, requires_grad=True))
"""
@classmethod
def _primitive_bind_call(cls, *args, **kwargs):
# accepts no wires, so bypasses the wire processing.
return cls._primitive.bind(*args, **kwargs)
def __init__( # pylint: disable=too-many-arguments
self, hamiltonian, time, n=1, order=1, check_hermitian=True, id=None
):
r"""Initialize the TrotterProduct class"""
if order <= 0 or order != 1 and order % 2 != 0:
raise ValueError(
f"The order of a TrotterProduct must be 1 or a positive even integer, got {order}."
)
if isinstance(hamiltonian, qml.ops.LinearCombination):
coeffs, ops = hamiltonian.terms()
if len(coeffs) < 2:
raise ValueError(
"There should be at least 2 terms in the Hamiltonian. Otherwise use `qml.exp`"
)
if qml.QueuingManager.recording():
qml.QueuingManager.remove(hamiltonian)
hamiltonian = qml.dot(coeffs, ops)
if isinstance(hamiltonian, SProd):
if qml.QueuingManager.recording():
qml.QueuingManager.remove(hamiltonian)
hamiltonian = hamiltonian.simplify()
if len(hamiltonian.terms()[0]) < 2:
raise ValueError(
"There should be at least 2 terms in the Hamiltonian. Otherwise use `qml.exp`"
)
if not isinstance(hamiltonian, Sum):
raise TypeError(
f"The given operator must be a PennyLane ~.Sum or ~.SProd, got {hamiltonian}"
)
if check_hermitian:
for op in hamiltonian.operands:
if not op.is_hermitian:
raise ValueError(
"One or more of the terms in the Hamiltonian may not be Hermitian"
)
self._hyperparameters = {
"n": n,
"order": order,
"base": hamiltonian,
"check_hermitian": check_hermitian,
}
super().__init__(*hamiltonian.data, time, wires=hamiltonian.wires, id=id)
[docs] def map_wires(self, wire_map: dict):
# pylint: disable=protected-access
new_op = copy.deepcopy(self)
new_op._wires = Wires([wire_map.get(wire, wire) for wire in self.wires])
new_op._hyperparameters["base"] = qml.map_wires(new_op._hyperparameters["base"], wire_map)
return new_op
[docs] def queue(self, context=qml.QueuingManager):
context.remove(self.hyperparameters["base"])
context.append(self)
return self
[docs] def resources(self) -> qml.resource.Resources:
r"""The resource requirements for a given instance of the Suzuki-Trotter product.
Returns:
:class:`~.resource.Resources`: The resources for an instance of ``TrotterProduct``.
"""
with qml.QueuingManager.stop_recording():
decomp = self.compute_decomposition(*self.parameters, **self.hyperparameters)
num_wires = len(self.wires)
num_gates = len(decomp)
depth = qml.tape.QuantumScript(ops=decomp).graph.get_depth()
gate_types = defaultdict(int)
gate_sizes = defaultdict(int)
for op in decomp:
gate_types[op.name] += 1
gate_sizes[len(op.wires)] += 1
return Resources(num_wires, num_gates, gate_types, gate_sizes, depth)
[docs] def error(
self, method: str = "commutator-bound", fast: bool = True
): # pylint: disable=arguments-differ
# pylint: disable=protected-access
r"""Compute an *upper-bound* on the spectral norm error for approximating the
time-evolution of the base Hamiltonian using the Suzuki-Trotter product formula.
The error in the Suzuki-Trotter product formula is defined as
.. math:: || \ e^{iHt} - \left [S_{m}(t / n) \right ]^{n} \ ||,
Where the norm :math:`||\cdot||` is the spectral norm. This function supports two methods
from literature for upper-bounding the error, the "one-norm" error bound and the "commutator"
error bound.
**Example:**
The "one-norm" error bound can be computed by passing the kwarg :code:`method="one-norm-bound"`, and
is defined according to Section 2.3 (lemma 6, equation 22 and 23) of
`Childs et al. (2021) <https://arxiv.org/abs/1912.08854>`_.
>>> hamiltonian = qml.dot([1.0, 0.5, -0.25], [qml.X(0), qml.Y(0), qml.Z(0)])
>>> op = qml.TrotterProduct(hamiltonian, time=0.01, order=2)
>>> op.error(method="one-norm-bound")
SpectralNormError(8.039062500000003e-06)
The "commutator" error bound can be computed by passing the kwarg :code:`method="commutator-bound"`, and
is defined according to Appendix C (equation 189) `Childs et al. (2021) <https://arxiv.org/abs/1912.08854>`_.
>>> hamiltonian = qml.dot([1.0, 0.5, -0.25], [qml.X(0), qml.Y(0), qml.Z(0)])
>>> op = qml.TrotterProduct(hamiltonian, time=0.01, order=2)
>>> op.error(method="commutator-bound")
SpectralNormError(6.166666666666668e-06)
Args:
method (str, optional): Options include "one-norm-bound" and "commutator-bound" and specify the
method with which the error is computed. Defaults to "commutator-bound".
fast (bool, optional): Uses more approximations to speed up computation. Defaults to True.
Raises:
ValueError: The method is not supported.
Returns:
SpectralNormError: The spectral norm error.
"""
base_unitary = self.hyperparameters["base"]
t, p, n = (self.parameters[-1], self.hyperparameters["order"], self.hyperparameters["n"])
parameters = [t] + base_unitary.parameters
if any(
qml.math.get_interface(param) == "tensorflow" for param in parameters
): # TODO: Add TF support
raise TypeError(
"Calculating error bound for Tensorflow objects is currently not supported"
)
terms = base_unitary.operands
if method == "one-norm-bound":
return SpectralNormError(qml.resource.error._one_norm_error(terms, t, p, n, fast=fast))
if method == "commutator-bound":
return SpectralNormError(
qml.resource.error._commutator_error(terms, t, p, n, fast=fast)
)
raise ValueError(
f"The '{method}' method is not supported for computing the error. Please select a valid method for computing the error."
)
def _flatten(self):
"""Serialize the operation into trainable and non-trainable components.
Returns:
data, metadata: The trainable and non-trainable components.
See ``Operator._unflatten``.
The data component can be recursive and include other operations. For example, the trainable component of ``Adjoint(RX(1, wires=0))``
will be the operator ``RX(1, wires=0)``.
The metadata **must** be hashable. If the hyperparameters contain a non-hashable component, then this
method and ``Operator._unflatten`` should be overridden to provide a hashable version of the hyperparameters.
**Example:**
>>> op = qml.Rot(1.2, 2.3, 3.4, wires=0)
>>> qml.Rot._unflatten(*op._flatten())
Rot(1.2, 2.3, 3.4, wires=[0])
>>> op = qml.PauliRot(1.2, "XY", wires=(0,1))
>>> qml.PauliRot._unflatten(*op._flatten())
PauliRot(1.2, XY, wires=[0, 1])
Operators that have trainable components that differ from their ``Operator.data`` must implement their own
``_flatten`` methods.
>>> op = qml.ctrl(qml.U2(3.4, 4.5, wires="a"), ("b", "c") )
>>> op._flatten()
((U2(3.4, 4.5, wires=['a']),),
(Wires(['b', 'c']), (True, True), Wires([])))
"""
hamiltonian = self.hyperparameters["base"]
time = self.data[-1]
hashable_hyperparameters = tuple(
item for item in self.hyperparameters.items() if item[0] != "base"
)
return (hamiltonian, time), hashable_hyperparameters
@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'])
"""
return cls(*data, **dict(metadata))
[docs] @staticmethod
def compute_decomposition(*args, **kwargs):
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
"""
time = args[-1]
n = kwargs["n"]
order = kwargs["order"]
ops = kwargs["base"].operands
decomp = _recursive_expression(time / n, order, ops)[::-1] * n
if qml.QueuingManager.recording():
for op in decomp: # apply operators in reverse order of expression
qml.apply(op)
return decomp
[docs]class TrotterizedQfunc(Operation):
r"""An operation representing the Suzuki-Trotter product approximation applied to a set of
operations defined in a function.
The Suzuki-Trotter product formula provides a method to approximate the matrix exponential of
a Hamiltonian expressed as a linear combination of terms which in general do not commute. Consider
the Hamiltonian :math:`H = \Sigma^{N}_{j=0} O_{j}`, the product formula is constructed using
symmetrized products of the terms in the Hamiltonian. The symmetrized products of order
:math:`m \in [1, 2, 4, ..., 2k]` with :math:`k \in \mathbb{N}` are given by:
.. math::
\begin{align}
S_{1}(t) &= \Pi_{j=0}^{N} \ e^{i t O_{j}} \\
S_{2}(t) &= \Pi_{j=0}^{N} \ e^{i \frac{t}{2} O_{j}} \cdot \Pi_{j=N}^{0} \ e^{i \frac{t}{2} O_{j}} \\
&\vdots \\
S_{m}(t) &= S_{m-2}(p_{m}t)^{2} \cdot S_{m-2}((1-4p_{m})t) \cdot S_{m-2}(p_{m}t)^{2},
\end{align}
where the coefficient is :math:`p_{m} = 1 / (4 - \sqrt[m - 1]{4})`. The :math:`m`th order,
:math:`n`-step Suzuki-Trotter approximation is then defined as:
.. math:: e^{iHt} \approx \left [S_{m}(t / n) \right ]^{n}.
For more details see `J. Math. Phys. 32, 400 (1991) <https://pubs.aip.org/aip/jmp/article-abstract/32/2/400/229229>`_.
Suppose we have direct access to the operators which represent the exponentiated terms of
a hamiltonian:
.. math:: \{ \hat{U}_{j} = e^{i t O_{j}}, j \in [1, N] \}.
Given a quantum circuit which uses these :math:`\hat{U}_{j}` operators to represent the
first order expansion :math:`S_{1}(t)`; this class expands it to any higher order Suzuki-Trotter product.
.. warning::
:code:`TrotterizedQfunc` requires that the input function has a very specific function signature.
The first argument should be a time parameter which will be modified according to the Suzuki-Trotter
product formula. The wires required by the circuit should be either the last explicit argument or the
first keyword argument. :code:`qfunc((time, arg1, ..., arg_n, wires=[...], kwarg_1, ..., kwarg_n))`
.. warning::
:code:`TrotterizedQfunc` currently does not support pickling. Instead please decompose the operation
first before attempting to pickle the quantum circuit.
Args:
time (float): the time of evolution, namely the parameter :math:`t` in :math:`e^{iHt}`
*trainable_args (tuple): the trainable arguments of the first-order expansion function
qfunc (Callable): the first-order expansion given as a callable function which queues operations
wires (Iterable): the set of wires the operation will act upon (should be identical to qfunc wires)
n (int): an integer representing the number of Trotter steps to perform
order (int): an integer (:math:`m`) representing the order of the approximation (must be 1 or even)
reverse (bool): if true, reverse the order of the operations queued by :code:`qfunc`
**non_trainable_kwargs (dict): non-trainable keyword arguments of the first-order expansion function
Raises:
ValueError: A qfunc must be provided to be trotterized.
See also :func:`~.trotterize`.
**Example**
.. code-block:: python3
def first_order_expansion(time, theta, phi, wires=[0, 1, 2], flip=False):
"This is the first order expansion (U_1)."
qml.RX(time*theta, wires[0])
qml.RY(time*phi, wires[1])
if flip:
qml.CNOT(wires=wires[:2])
@qml.qnode(qml.device("default.qubit"))
def my_circuit(time, angles, num_trotter_steps):
qml.TrotterizedQfunc(
time,
*angles,
qfunc=first_order_expansion,
n=num_trotter_steps,
order=2,
wires=['a', 'b'],
flip=True,
)
return qml.state()
We can visualize the circuit to see the Suzuki-Trotter product formula being applied:
>>> time = 0.1
>>> angles = (0.12, -3.45)
>>> print(qml.draw(my_circuit, level="device")(time, angles, num_trotter_steps=1))
a: ──RX(0.01)──╭●─╭●──RX(0.01)──┤ State
b: ──RY(-0.17)─╰X─╰X──RY(-0.17)─┤ State
>>>
>>> print(qml.draw(my_circuit, level="device")(time, angles, num_trotter_steps=3))
a: ──RX(0.00)──╭●─╭●──RX(0.00)───RX(0.00)──╭●─╭●──RX(0.00)───RX(0.00)──╭●─╭●──RX(0.00)──┤ State
b: ──RY(-0.06)─╰X─╰X──RY(-0.06)──RY(-0.06)─╰X─╰X──RY(-0.06)──RY(-0.06)─╰X─╰X──RY(-0.06)─┤ State
"""
# pylint: disable = too-many-arguments
def __init__(
self,
time,
*trainable_args,
qfunc=None,
n=1,
order=2,
reverse=False,
id=None,
**non_trainable_kwargs,
):
# This class requires the input function (qfunc) has a very specific
# signature: qfunc(time, arg1, ..., arg_n, wires=[...], kwarg_1, ..., kwarg_n)
if qfunc is None:
raise ValueError("The qfunc must be provided to be trotterized.")
if order <= 0 or order != 1 and order % 2 != 0:
raise ValueError(f"The order must be 1 or a positive even integer, got {order}.")
try:
wires = non_trainable_kwargs.pop("wires")
except KeyError:
wires = trainable_args[-1]
trainable_args = trainable_args[:-1] # exclude the wires from the args
self._hyperparameters = non_trainable_kwargs
self._hyperparameters["n"] = n
self._hyperparameters["order"] = order
self._hyperparameters["qfunc"] = qfunc
self._hyperparameters["reverse"] = reverse
super().__init__(time, *trainable_args, wires=wires, id=id)
[docs] def decomposition(self) -> list[Operator]:
"""The decomposition"""
n = self.hyperparameters["n"]
order = self.hyperparameters["order"]
qfunc = self.hyperparameters["qfunc"]
reverse = self.hyperparameters["reverse"]
time = self.parameters[0]
qfunc_args = self.parameters[1:]
base_hyper_params = ("n", "order", "qfunc", "reverse")
qfunc_kwargs = {k: v for k, v in self.hyperparameters.items() if not k in base_hyper_params}
decomp = (
_recursive_qfunc(
time / n, order, qfunc, self.wires, reverse, *qfunc_args, **qfunc_kwargs
)
* n
)
if qml.QueuingManager.recording():
for op in decomp:
qml.apply(op)
return decomp
def _flatten(self):
"""Serialize the operation into trainable and non-trainable components.
Returns:
data, metadata: The trainable and non-trainable components.
See ``Operator._unflatten``.
The data component can be recursive and include other operations. For example, the trainable component of ``Adjoint(RX(1, wires=0))``
will be the operator ``RX(1, wires=0)``.
The metadata **must** be hashable. If the hyperparameters contain a non-hashable component, then this
method and ``Operator._unflatten`` should be overridden to provide a hashable version of the hyperparameters.
"""
hashable_hyperparameters = tuple(self.hyperparameters.items()) + (("wires", self.wires),)
return self.data, hashable_hyperparameters
@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``.
"""
return cls(*data, **dict(metadata))
[docs]def trotterize(qfunc, n=1, order=2, reverse=False):
r"""Generates higher order Suzuki-Trotter product formulas from a set of
operations defined in a function.
The Suzuki-Trotter product formula provides a method to approximate the matrix exponential of a
Hamiltonian expressed as a linear combination of terms which in general do not commute. Consider
the Hamiltonian :math:`H = \Sigma^{N}_{j=0} O_{j}`, the product formula is constructed using
symmetrized products of the terms in the Hamiltonian. The symmetrized products of order
:math:`m \in [1, 2, 4, ..., 2k]` with :math:`k \in \mathbb{N}` are given by:
.. math::
\begin{align}
S_{1}(t) &= \Pi_{j=0}^{N} \ e^{i t O_{j}} \\
S_{2}(t) &= \Pi_{j=0}^{N} \ e^{i \frac{t}{2} O_{j}} \cdot \Pi_{j=N}^{0} \ e^{i \frac{t}{2} O_{j}} \\
&\vdots \\
S_{m}(t) &= S_{m-2}(p_{m}t)^{2} \cdot S_{m-2}((1-4p_{m})t) \cdot S_{m-2}(p_{m}t)^{2},
\end{align}
where the coefficient is :math:`p_{m} = 1 / (4 - \sqrt[m - 1]{4})`. The :math:`m`th order,
:math:`n`-step Suzuki-Trotter approximation is then defined as:
.. math:: e^{iHt} \approx \left [S_{m}(t / n) \right ]^{n}.
For more details see `J. Math. Phys. 32, 400 (1991) <https://pubs.aip.org/aip/jmp/article-abstract/32/2/400/229229>`_.
Suppose we have direct access to the operators which represent the exponentiated terms of
a Hamiltonian:
.. math:: \{ \hat{U}_{j} = e^{i t O_{j}}, j \in [1, N] \}.
Given a quantum circuit which uses these :math:`\hat{U}_{j}` operators to represent the
first order expansion :math:`S_{1}(t)`, this function expands it to any higher order Suzuki-Trotter product.
.. warning::
:code:`trotterize()` requires the :code:`qfunc` argument to be a function with a specific call
signature. The first argument of the :code:`qfunc` function should be a time parameter which will be modified according to the
Suzuki-Trotter product formula. The wires required by the circuit should be either the last
positional argument or the first keyword argument:
:code:`qfunc((time, arg1, ..., arg_n, wires=[...], kwarg_1, ..., kwarg_n))`
Args:
qfunc (Callable): the first-order expansion given as a callable function which queues operations
n (int): an integer representing the number of Trotter steps to perform
order (int): an integer (:math:`m`) representing the order of the approximation (must be 1 or even)
reverse (bool): if true, reverse the order of the operations queued by :code:`qfunc`
name (str): an optional name for the instance
Returns:
(Callable): a function with the same signature as :code:`qfunc`, when called it queues an instance of
:class:`~.TrotterizedQfunc`
**Example**
.. code-block:: python3
def first_order_expansion(time, theta, phi, wires, flip=False):
"This is the first order expansion (U_1)."
qml.RX(time*theta, wires[0])
qml.RY(time*phi, wires[1])
if flip:
qml.CNOT(wires=wires[:2])
@qml.qnode(qml.device("default.qubit"))
def my_circuit(time, theta, phi, num_trotter_steps):
qml.trotterize(
first_order_expansion,
n=num_trotter_steps,
order=2,
)(time, theta, phi, wires=['a', 'b'], flip=True)
return qml.state()
We can visualize the circuit to see the Suzuki-Trotter product formula being applied:
>>> time = 0.1
>>> theta, phi = (0.12, -3.45)
>>>
>>> print(qml.draw(my_circuit, level="device")(time, theta, phi, num_trotter_steps=1))
a: ──RX(0.01)──╭●─╭●──RX(0.01)──┤ State
b: ──RY(-0.17)─╰X─╰X──RY(-0.17)─┤ State
>>>
>>> print(qml.draw(my_circuit, level="device")(time, theta, phi, num_trotter_steps=3))
a: ──RX(0.00)──╭●─╭●──RX(0.00)───RX(0.00)──╭●─╭●──RX(0.00)───RX(0.00)──╭●─╭●──RX(0.00)──┤ State
b: ──RY(-0.06)─╰X─╰X──RY(-0.06)──RY(-0.06)─╰X─╰X──RY(-0.06)──RY(-0.06)─╰X─╰X──RY(-0.06)─┤ State
"""
@wraps(qfunc)
def wrapper(time, *args, **kwargs):
special_keys = ["n", "order", "qfunc", "reverse"]
if any(key in kwargs for key in special_keys):
raise ValueError(
f"Cannot use any of the specailized names:\n {special_keys}\nas keyword"
f"arguments to the qfunc: {kwargs}"
)
return TrotterizedQfunc(
time, *args, qfunc=qfunc, n=n, order=order, reverse=reverse, **kwargs
)
return wrapper
@qml.QueuingManager.stop_recording()
def _recursive_qfunc(time, order, qfunc, wires, reverse, *qfunc_args, **qfunc_kwargs):
"""Generate a list of operations using the
recursive expression which defines the Trotter product.
Args:
time (float): the evolution 'time'
order (int): the order of the Trotter expansion
ops (Iterable(~.Operators)): a list of terms in the Hamiltonian
Returns:
list: the approximation as a product of exponentials of the Hamiltonian terms
"""
if order == 1:
tape = qml.tape.make_qscript(qfunc)(time, *qfunc_args, wires=wires, **qfunc_kwargs)
return tape.operations[::-1] if reverse else tape.operations
if order == 2:
tape = qml.tape.make_qscript(qfunc)(time / 2, *qfunc_args, wires=wires, **qfunc_kwargs)
return (
tape.operations[::-1] + tape.operations
if reverse
else tape.operations + tape.operations[::-1]
)
scalar_1 = _scalar(order)
scalar_2 = 1 - 4 * scalar_1
ops_lst_1 = _recursive_qfunc(
scalar_1 * time, order - 2, qfunc, wires, reverse, *qfunc_args, **qfunc_kwargs
)
ops_lst_2 = _recursive_qfunc(
scalar_2 * time, order - 2, qfunc, wires, reverse, *qfunc_args, **qfunc_kwargs
)
return (2 * ops_lst_1) + ops_lst_2 + (2 * ops_lst_1)
_modules/pennylane/templates/subroutines/trotter
Download Python script
Download Notebook
View on GitHub