Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

qml.QDrift template #4671

Merged
merged 39 commits into from
Oct 20, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
39 commits
Select commit Hold shift + click to select a range
89c5e12
QDrift
KetpuntoG Oct 12, 2023
e533cef
Update qdrift.py
KetpuntoG Oct 12, 2023
8c115e3
Update qdrift.py
KetpuntoG Oct 12, 2023
fd46345
Update qdrift.py
KetpuntoG Oct 12, 2023
5ac1d16
sign typo
KetpuntoG Oct 13, 2023
49c747d
black
KetpuntoG Oct 13, 2023
96d3c7a
Update qdrift.py
KetpuntoG Oct 13, 2023
83d1a3c
Update qdrift.py
KetpuntoG Oct 13, 2023
1447038
Update qdrift.py
KetpuntoG Oct 13, 2023
2701987
add figure
soranjh Oct 15, 2023
214fdff
add figure
soranjh Oct 15, 2023
fcf1d9e
add figure
soranjh Oct 15, 2023
9880bd4
use random from math
soranjh Oct 15, 2023
ac8fda6
modify docstring
soranjh Oct 15, 2023
38965f7
[skip ci] correct seed
soranjh Oct 15, 2023
945e22d
[skip ci] modify docstring
soranjh Oct 15, 2023
96edcdd
correct math
soranjh Oct 15, 2023
cd735d1
[skip ci] modify docs
soranjh Oct 15, 2023
e8f433c
Added error computation
Jaybsoni Oct 17, 2023
c52774b
adding base tests
Jaybsoni Oct 18, 2023
bcbe8bb
Added more integration tests
Jaybsoni Oct 18, 2023
5e2f998
lint
Jaybsoni Oct 18, 2023
092c20d
Merge branch 'master' into qdrift
Jaybsoni Oct 18, 2023
623f3b7
update changelog
soranjh Oct 19, 2023
ef6437d
update image
soranjh Oct 19, 2023
ceea397
add grad tests
Jaybsoni Oct 19, 2023
110e822
add interface tests for execution
Jaybsoni Oct 19, 2023
d84b04c
Merge branch 'master' into qdrift
Jaybsoni Oct 19, 2023
2cfa58f
Added tests
Jaybsoni Oct 19, 2023
48f8066
final tests
Jaybsoni Oct 19, 2023
63a0c1e
Merge branch 'master' into qdrift
Jaybsoni Oct 19, 2023
464057c
fix error func + lint
Jaybsoni Oct 19, 2023
7a89ce7
lint
Jaybsoni Oct 19, 2023
6152f4d
Merge branch 'master' into qdrift
Jaybsoni Oct 19, 2023
a7b3864
lint + code review commentsa
Jaybsoni Oct 20, 2023
bb71916
lint docstring
Jaybsoni Oct 20, 2023
c676250
lint
Jaybsoni Oct 20, 2023
34c6cd2
Merge branch 'master' into qdrift
Jaybsoni Oct 20, 2023
559ee72
remove flaky import
Jaybsoni Oct 20, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file added doc/_static/templates/subroutines/qdrift.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
4 changes: 4 additions & 0 deletions doc/introduction/templates.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
26 changes: 26 additions & 0 deletions doc/releases/changelog-dev.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions pennylane/templates/subroutines/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
289 changes: 289 additions & 0 deletions pennylane/templates/subroutines/qdrift.py
Original file line number Diff line number Diff line change
@@ -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)
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
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
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved

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
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
):
r"""Initialize the QDrift class"""

if isinstance(hamiltonian, Hamiltonian):
coeffs, ops = hamiltonian.terms()

elif isinstance(hamiltonian, Sum):
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
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)
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved

self._hyperparameters = {
"n": n,
"seed": seed,
"base": hamiltonian,
"decomposition": decomposition,
}
super().__init__(time, wires=hamiltonian.wires, id=id)

@classmethod
def _unflatten(cls, data, metadata):
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
"""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():
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
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)
Loading
Loading