diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index 2b995645da5..b7bd88e9fcb 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -21,8 +21,29 @@ The `SingleExcitation` operation supports analytic gradients on hardware using only four expectation value calculations, following results from - [Kottmann et al.](https://arxiv.org/abs/2011.05938) + [Kottmann et al.](https://arxiv.org/abs/2011.05938) + * Added the `DoubleExcitation` four-qubit operation, which is useful for quantum + chemistry applications. [(#1123)](https://github.com/PennyLaneAI/pennylane/pull/1123) + + It can be used to perform an SO(2) rotation in the subspace + spanned by the states :math:`|1100\rangle` and :math:`|0011\rangle`. + For example, the following circuit performs the transformation + :math:`|1100\rangle\rightarrow \cos(\phi/2)|1100\rangle - \sin(\phi/2)|0011\rangle`: + + ```python + dev = qml.device('default.qubit', wires=2) + + @qml.qnode(dev) + def circuit(phi): + qml.PauliX(wires=0) + qml.PauliX(wires=1) + qml.DoubleExcitation(phi, wires=[0, 1, 2, 3]) + ``` + + The `DoubleExcitation` operation supports analytic gradients on hardware using only + four expectation value calculations, following results from + [Kottmann et al.](https://arxiv.org/abs/2011.05938). * Adds a new function ``qml.math.conj``. [(#1143)](https://github.com/PennyLaneAI/pennylane/pull/1143) diff --git a/doc/introduction/operations.rst b/doc/introduction/operations.rst index 9de7a57b87a..4a4e3e523b9 100644 --- a/doc/introduction/operations.rst +++ b/doc/introduction/operations.rst @@ -88,6 +88,10 @@ Qubit gates ~pennylane.SingleExcitation ~pennylane.SingleExcitationPlus ~pennylane.SingleExcitationMinus + ~pennylane.DoubleExcitation + ~pennylane.DoubleExcitationPlus + ~pennylane.DoubleExcitationMinus + :html:`` diff --git a/pennylane/devices/autograd_ops.py b/pennylane/devices/autograd_ops.py index 905df14b769..fa69464fb9a 100644 --- a/pennylane/devices/autograd_ops.py +++ b/pennylane/devices/autograd_ops.py @@ -1,226 +1,330 @@ -# Copyright 2018-2020 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. -r""" -Utility functions and numerical implementations of quantum operations for Autograd-based devices. -""" -from autograd import numpy as np -from numpy import kron -from pennylane.utils import pauli_eigs - -C_DTYPE = np.complex128 -R_DTYPE = np.float64 - -I = np.array([[1, 0], [0, 1]], dtype=C_DTYPE) -X = np.array([[0, 1], [1, 0]], dtype=C_DTYPE) -Y = np.array([[0j, -1j], [1j, 0j]], dtype=C_DTYPE) -Z = np.array([[1, 0], [0, -1]], dtype=C_DTYPE) - -II = np.eye(4, dtype=C_DTYPE) -ZZ = np.array(kron(Z, Z), dtype=C_DTYPE) - -IX = np.array(kron(I, X), dtype=C_DTYPE) -IY = np.array(kron(I, Y), dtype=C_DTYPE) -IZ = np.array(kron(I, Z), dtype=C_DTYPE) - -ZI = np.array(kron(Z, I), dtype=C_DTYPE) -ZX = np.array(kron(Z, X), dtype=C_DTYPE) -ZY = np.array(kron(Z, Y), dtype=C_DTYPE) - - -def PhaseShift(phi): - r"""One-qubit phase shift. - - Args: - phi (float): phase shift angle - - Returns: - array[complex]: diagonal part of the phase shift matrix - """ - return np.array([1.0, np.exp(1j * phi)]) - - -def ControlledPhaseShift(phi): - r"""Two-qubit controlled phase shift. - - Args: - phi (float): phase shift angle - - Returns: - array[complex]: diagonal part of the controlled phase shift matrix - """ - return np.array([1.0, 1.0, 1.0, np.exp(1j * phi)]) - - -def RX(theta): - r"""One-qubit rotation about the x axis. - - Args: - theta (float): rotation angle - - Returns: - array[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_x \theta/2}` - """ - return np.cos(theta / 2) * I + 1j * np.sin(-theta / 2) * X - - -def RY(theta): - r"""One-qubit rotation about the y axis. - - Args: - theta (float): rotation angle - - Returns: - array[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_y \theta/2}` - """ - return np.cos(theta / 2) * I + 1j * np.sin(-theta / 2) * Y - - -def RZ(theta): - r"""One-qubit rotation about the z axis. - - Args: - theta (float): rotation angle - - Returns: - array[complex]: the diagonal part of the rotation matrix :math:`e^{-i \sigma_z \theta/2}` - """ - p = np.exp(-0.5j * theta) - return np.array([p, np.conj(p)]) - - -def Rot(a, b, c): - r"""Arbitrary one-qubit rotation using three Euler angles. - - Args: - a,b,c (float): rotation angles - - Returns: - array[complex]: unitary 2x2 rotation matrix ``rz(c) @ ry(b) @ rz(a)`` - """ - return np.diag(RZ(c)) @ RY(b) @ np.diag(RZ(a)) - - -def CRX(theta): - r"""Two-qubit controlled rotation about the x axis. - - Args: - theta (float): rotation angle - - Returns: - array[complex]: unitary 4x4 rotation matrix - :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_x(\theta)` - """ - return ( - np.cos(theta / 4) ** 2 * II - - 1j * np.sin(theta / 2) / 2 * IX - + np.sin(theta / 4) ** 2 * ZI - + 1j * np.sin(theta / 2) / 2 * ZX - ) - - -def CRY(theta): - r"""Two-qubit controlled rotation about the y axis. - - Args: - theta (float): rotation angle - Returns: - array[complex]: unitary 4x4 rotation matrix :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_y(\theta)` - """ - return ( - np.cos(theta / 4) ** 2 * II - - 1j * np.sin(theta / 2) / 2 * IY - + np.sin(theta / 4) ** 2 * ZI - + 1j * np.sin(theta / 2) / 2 * ZY - ) - - -def CRZ(theta): - r"""Two-qubit controlled rotation about the z axis. - - Args: - theta (float): rotation angle - Returns: - array[complex]: diagonal part of the 4x4 rotation matrix - :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_z(\theta)` - """ - p = np.exp(-0.5j * theta) - return np.array([1.0, 1.0, p, np.conj(p)]) - - -def CRot(a, b, c): - r"""Arbitrary two-qubit controlled rotation using three Euler angles. - - Args: - a,b,c (float): rotation angles - Returns: - array[complex]: unitary 4x4 rotation matrix - :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R(a,b,c)` - """ - return np.diag(CRZ(c)) @ (CRY(b) @ np.diag(CRZ(a))) - - -def MultiRZ(theta, n): - r"""Arbitrary multi Z rotation. - - Args: - theta (float): rotation angle :math:`\theta` - wires (Sequence[int] or int): the wires the operation acts on - Returns: - array[complex]: diagonal part of the multi-qubit rotation matrix - """ - return np.exp(-1j * theta / 2 * pauli_eigs(n)) - - -def SingleExcitation(phi): - r"""Single excitation rotation. - - Args: - phi (float): rotation angle - - Returns: - array[float]: Single excitation rotation matrix - """ - c = np.cos(phi / 2) - s = np.sin(phi / 2) - return np.array([[1, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, 1]]) - - -def SingleExcitationPlus(phi): - r"""Single excitation rotation with positive phase-shift outside the rotation subspace. - - Args: - phi (float): rotation angle - - Returns: - array[complex]: Single excitation rotation matrix with positive phase-shift - """ - c = np.cos(phi / 2) - s = np.sin(phi / 2) - e = np.exp(1j * phi / 2) - return np.array([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) - - -def SingleExcitationMinus(phi): - r"""Single excitation rotation with negative phase-shift outside the rotation subspace. - - Args: - phi (float): rotation angle - - Returns: - array[complex]: Single excitation rotation matrix with negative phase-shift - """ - c = np.cos(phi / 2) - s = np.sin(phi / 2) - e = np.exp(-1j * phi / 2) - return np.array([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) +# Copyright 2018-2020 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. +r""" +Utility functions and numerical implementations of quantum operations for Autograd-based devices. +""" +from autograd import numpy as np +from numpy import kron +from pennylane.utils import pauli_eigs + +C_DTYPE = np.complex128 +R_DTYPE = np.float64 + +I = np.array([[1, 0], [0, 1]], dtype=C_DTYPE) +X = np.array([[0, 1], [1, 0]], dtype=C_DTYPE) +Y = np.array([[0j, -1j], [1j, 0j]], dtype=C_DTYPE) +Z = np.array([[1, 0], [0, -1]], dtype=C_DTYPE) + +II = np.eye(4, dtype=C_DTYPE) +ZZ = np.array(kron(Z, Z), dtype=C_DTYPE) + +IX = np.array(kron(I, X), dtype=C_DTYPE) +IY = np.array(kron(I, Y), dtype=C_DTYPE) +IZ = np.array(kron(I, Z), dtype=C_DTYPE) + +ZI = np.array(kron(Z, I), dtype=C_DTYPE) +ZX = np.array(kron(Z, X), dtype=C_DTYPE) +ZY = np.array(kron(Z, Y), dtype=C_DTYPE) + + +def PhaseShift(phi): + r"""One-qubit phase shift. + + Args: + phi (float): phase shift angle + + Returns: + array[complex]: diagonal part of the phase shift matrix + """ + return np.array([1.0, np.exp(1j * phi)]) + + +def ControlledPhaseShift(phi): + r"""Two-qubit controlled phase shift. + + Args: + phi (float): phase shift angle + + Returns: + array[complex]: diagonal part of the controlled phase shift matrix + """ + return np.array([1.0, 1.0, 1.0, np.exp(1j * phi)]) + + +def RX(theta): + r"""One-qubit rotation about the x axis. + + Args: + theta (float): rotation angle + + Returns: + array[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_x \theta/2}` + """ + return np.cos(theta / 2) * I + 1j * np.sin(-theta / 2) * X + + +def RY(theta): + r"""One-qubit rotation about the y axis. + + Args: + theta (float): rotation angle + + Returns: + array[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_y \theta/2}` + """ + return np.cos(theta / 2) * I + 1j * np.sin(-theta / 2) * Y + + +def RZ(theta): + r"""One-qubit rotation about the z axis. + + Args: + theta (float): rotation angle + + Returns: + array[complex]: the diagonal part of the rotation matrix :math:`e^{-i \sigma_z \theta/2}` + """ + p = np.exp(-0.5j * theta) + return np.array([p, np.conj(p)]) + + +def Rot(a, b, c): + r"""Arbitrary one-qubit rotation using three Euler angles. + + Args: + a,b,c (float): rotation angles + + Returns: + array[complex]: unitary 2x2 rotation matrix ``rz(c) @ ry(b) @ rz(a)`` + """ + return np.diag(RZ(c)) @ RY(b) @ np.diag(RZ(a)) + + +def CRX(theta): + r"""Two-qubit controlled rotation about the x axis. + + Args: + theta (float): rotation angle + + Returns: + array[complex]: unitary 4x4 rotation matrix + :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_x(\theta)` + """ + return ( + np.cos(theta / 4) ** 2 * II + - 1j * np.sin(theta / 2) / 2 * IX + + np.sin(theta / 4) ** 2 * ZI + + 1j * np.sin(theta / 2) / 2 * ZX + ) + + +def CRY(theta): + r"""Two-qubit controlled rotation about the y axis. + + Args: + theta (float): rotation angle + Returns: + array[complex]: unitary 4x4 rotation matrix :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_y(\theta)` + """ + return ( + np.cos(theta / 4) ** 2 * II + - 1j * np.sin(theta / 2) / 2 * IY + + np.sin(theta / 4) ** 2 * ZI + + 1j * np.sin(theta / 2) / 2 * ZY + ) + + +def CRZ(theta): + r"""Two-qubit controlled rotation about the z axis. + + Args: + theta (float): rotation angle + Returns: + array[complex]: diagonal part of the 4x4 rotation matrix + :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_z(\theta)` + """ + p = np.exp(-0.5j * theta) + return np.array([1.0, 1.0, p, np.conj(p)]) + + +def CRot(a, b, c): + r"""Arbitrary two-qubit controlled rotation using three Euler angles. + + Args: + a,b,c (float): rotation angles + Returns: + array[complex]: unitary 4x4 rotation matrix + :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R(a,b,c)` + """ + return np.diag(CRZ(c)) @ (CRY(b) @ np.diag(CRZ(a))) + + +def MultiRZ(theta, n): + r"""Arbitrary multi Z rotation. + + Args: + theta (float): rotation angle :math:`\theta` + wires (Sequence[int] or int): the wires the operation acts on + Returns: + array[complex]: diagonal part of the multi-qubit rotation matrix + """ + return np.exp(-1j * theta / 2 * pauli_eigs(n)) + + +def SingleExcitation(phi): + r"""Single excitation rotation. + + Args: + phi (float): rotation angle + + Returns: + array[float]: Single excitation rotation matrix + """ + c = np.cos(phi / 2) + s = np.sin(phi / 2) + return np.array([[1, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, 1]]) + + +def SingleExcitationPlus(phi): + r"""Single excitation rotation with positive phase-shift outside the rotation subspace. + + Args: + phi (float): rotation angle + + Returns: + array[complex]: Single excitation rotation matrix with positive phase-shift + """ + c = np.cos(phi / 2) + s = np.sin(phi / 2) + e = np.exp(1j * phi / 2) + return np.array([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) + + +def SingleExcitationMinus(phi): + r"""Single excitation rotation with negative phase-shift outside the rotation subspace. + + Args: + phi (float): rotation angle + + Returns: + array[complex]: Single excitation rotation matrix with negative phase-shift + """ + c = np.cos(phi / 2) + s = np.sin(phi / 2) + e = np.exp(-1j * phi / 2) + return np.array([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) + + +def DoubleExcitation(phi): + r"""Double excitation rotation. + + Args: + phi (float): rotation angle + Returns: + array[float]: Double excitation rotation matrix + + """ + c = np.cos(phi / 2) + s = np.sin(phi / 2) + + U = [ + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, c, 0, 0, 0, 0, 0, 0, 0, 0, -s, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, s, 0, 0, 0, 0, 0, 0, 0, 0, c, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], + ] + + return np.array(U) + + +def DoubleExcitationPlus(phi): + r"""Double excitation rotation with positive phase-shift. + + Args: + phi (float): rotation angle + Returns: + array[complex]: rotation matrix + """ + + c = np.cos(phi / 2) + s = np.sin(phi / 2) + e = np.exp(1j * phi / 2) + + U = [ + [e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, c, 0, 0, 0, 0, 0, 0, 0, 0, -s, 0, 0, 0], + [0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0], + [0, 0, 0, s, 0, 0, 0, 0, 0, 0, 0, 0, c, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e], + ] + + return np.array(U) + + +def DoubleExcitationMinus(phi): + r"""Double excitation rotation with negative phase-shift. + + Args: + phi (float): rotation angle + Returns: + array[complex]: rotation matrix + """ + + c = np.cos(phi / 2) + s = np.sin(phi / 2) + e = np.exp(-1j * phi / 2) + + U = [ + [e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, c, 0, 0, 0, 0, 0, 0, 0, 0, -s, 0, 0, 0], + [0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0], + [0, 0, 0, s, 0, 0, 0, 0, 0, 0, 0, 0, c, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e], + ] + + return np.array(U) diff --git a/pennylane/devices/default_mixed.py b/pennylane/devices/default_mixed.py index 3a5a126b3bf..920218066a7 100644 --- a/pennylane/devices/default_mixed.py +++ b/pennylane/devices/default_mixed.py @@ -95,6 +95,9 @@ class DefaultMixed(QubitDevice): "SingleExcitation", "SingleExcitationPlus", "SingleExcitationMinus", + "DoubleExcitation", + "DoubleExcitationPlus", + "DoubleExcitationMinus", } def __init__(self, wires, *, shots=None, cache=0): diff --git a/pennylane/devices/default_qubit.py b/pennylane/devices/default_qubit.py index 1c681a13a57..7d6d86f91f7 100644 --- a/pennylane/devices/default_qubit.py +++ b/pennylane/devices/default_qubit.py @@ -124,6 +124,9 @@ class DefaultQubit(QubitDevice): "SingleExcitation", "SingleExcitationPlus", "SingleExcitationMinus", + "DoubleExcitation", + "DoubleExcitationPlus", + "DoubleExcitationMinus", } observables = {"PauliX", "PauliY", "PauliZ", "Hadamard", "Hermitian", "Identity"} diff --git a/pennylane/devices/default_qubit_autograd.py b/pennylane/devices/default_qubit_autograd.py index ce3caee2c0e..89291730fd2 100644 --- a/pennylane/devices/default_qubit_autograd.py +++ b/pennylane/devices/default_qubit_autograd.py @@ -1,178 +1,181 @@ -# Copyright 2018-2020 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. -"""This module contains an autograd implementation of the :class:`~.DefaultQubit` -reference plugin. -""" -from pennylane.operation import DiagonalOperation -from pennylane import numpy as np - -from pennylane.devices import DefaultQubit -from pennylane.devices import autograd_ops - - -class DefaultQubitAutograd(DefaultQubit): - """Simulator plugin based on ``"default.qubit"``, written using Autograd. - - **Short name:** ``default.qubit.autograd`` - - This device provides a pure-state qubit simulator written using Autograd. As a result, it - supports classical backpropagation as a means to compute the gradient. This can be faster than - the parameter-shift rule for analytic quantum gradients when the number of parameters to be - optimized is large. - - To use this device, you will need to install Autograd: - - .. code-block:: console - - pip install autograd - - **Example** - - The ``default.qubit.autograd`` is designed to be used with end-to-end classical backpropagation - (``diff_method="backprop"``) with the Autograd interface. This is the default method of - differentiation when creating a QNode with this device. - - Using this method, the created QNode is a 'white-box', and is - tightly integrated with your Autograd computation: - - >>> dev = qml.device("default.qubit.autograd", wires=1) - >>> @qml.qnode(dev, interface="autograd", diff_method="backprop") - ... def circuit(x): - ... qml.RX(x[1], wires=0) - ... qml.Rot(x[0], x[1], x[2], wires=0) - ... return qml.expval(qml.PauliZ(0)) - >>> weights = np.array([0.2, 0.5, 0.1]) - >>> grad_fn = qml.grad(circuit) - >>> print(grad_fn(weights)) - array([-2.2526717e-01 -1.0086454e+00 1.3877788e-17]) - - There are a couple of things to keep in mind when using the ``"backprop"`` - differentiation method for QNodes: - - * You must use the ``"autograd"`` interface for classical backpropagation, as Autograd is - used as the device backend. - - * Only exact expectation values, variances, and probabilities are differentiable. - When instantiating the device with ``analytic=False``, differentiating QNode - outputs will result in an error. - - Args: - wires (int): the number of wires to initialize the device with - shots (None, int): How many times the circuit should be evaluated (or sampled) to estimate - the expectation values. Defaults to ``None`` if not specified, which means that the device - returns analytical results. - analytic (bool): Indicates if the device should calculate expectations - and variances analytically. In non-analytic mode, the ``diff_method="backprop"`` - QNode differentiation method is not supported and it is recommended to consider - switching device to ``default.qubit`` and using ``diff_method="parameter-shift"``. - """ - - name = "Default qubit (Autograd) PennyLane plugin" - short_name = "default.qubit.autograd" - - parametric_ops = { - "PhaseShift": autograd_ops.PhaseShift, - "ControlledPhaseShift": autograd_ops.ControlledPhaseShift, - "RX": autograd_ops.RX, - "RY": autograd_ops.RY, - "RZ": autograd_ops.RZ, - "Rot": autograd_ops.Rot, - "CRX": autograd_ops.CRX, - "CRY": autograd_ops.CRY, - "CRZ": autograd_ops.CRZ, - "CRot": autograd_ops.CRot, - "MultiRZ": autograd_ops.MultiRZ, - "SingleExcitation": autograd_ops.SingleExcitation, - "SingleExcitationPlus": autograd_ops.SingleExcitationPlus, - "SingleExcitationMinus": autograd_ops.SingleExcitationMinus, - } - - C_DTYPE = np.complex128 - R_DTYPE = np.float64 - _dot = staticmethod(np.dot) - _abs = staticmethod(np.abs) - _reduce_sum = staticmethod(lambda array, axes: np.sum(array, axis=tuple(axes))) - _reshape = staticmethod(np.reshape) - _flatten = staticmethod(lambda array: array.flatten()) - _gather = staticmethod(lambda array, indices: array[indices]) - _einsum = staticmethod(np.einsum) - _cast = staticmethod(np.asarray) - _transpose = staticmethod(np.transpose) - _tensordot = staticmethod(np.tensordot) - _conj = staticmethod(np.conj) - _imag = staticmethod(np.imag) - _roll = staticmethod(np.roll) - _stack = staticmethod(np.stack) - - @staticmethod - def _asarray(array, dtype=None): - res = np.asarray(array, dtype=dtype) - - if res.dtype is np.dtype("O"): - return np.hstack(array).flatten().astype(dtype) - - return res - - def __init__(self, wires, *, shots=None): - super().__init__(wires, shots=shots, cache=0) - - # prevent using special apply methods for these gates due to slowdown in Autograd - # implementation - del self._apply_ops["PauliY"] - del self._apply_ops["Hadamard"] - del self._apply_ops["CZ"] - - @classmethod - def capabilities(cls): - capabilities = super().capabilities().copy() - capabilities.update( - passthru_interface="autograd", - supports_reversible_diff=False, - ) - return capabilities - - @staticmethod - def _scatter(indices, array, new_dimensions): - new_array = np.zeros(new_dimensions, dtype=array.dtype.type) - new_array[indices] = array - return new_array - - def _get_unitary_matrix(self, unitary): - """Return the matrix representing a unitary operation. - - Args: - unitary (~.Operation): a PennyLane unitary operation - - Returns: - array[complex]: Returns a 2D matrix representation of - the unitary in the computational basis, or, in the case of a diagonal unitary, - a 1D array representing the matrix diagonal. - """ - op_name = unitary.name.split(".inv")[0] - - if op_name in self.parametric_ops: - if op_name == "MultiRZ": - mat = self.parametric_ops[op_name](*unitary.parameters, len(unitary.wires)) - else: - mat = self.parametric_ops[op_name](*unitary.parameters) - - if unitary.inverse: - mat = self._transpose(self._conj(mat)) - - return mat - - if isinstance(unitary, DiagonalOperation): - return unitary.eigvals - - return unitary.matrix +# Copyright 2018-2020 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. +"""This module contains an autograd implementation of the :class:`~.DefaultQubit` +reference plugin. +""" +from pennylane.operation import DiagonalOperation +from pennylane import numpy as np + +from pennylane.devices import DefaultQubit +from pennylane.devices import autograd_ops + + +class DefaultQubitAutograd(DefaultQubit): + """Simulator plugin based on ``"default.qubit"``, written using Autograd. + + **Short name:** ``default.qubit.autograd`` + + This device provides a pure-state qubit simulator written using Autograd. As a result, it + supports classical backpropagation as a means to compute the gradient. This can be faster than + the parameter-shift rule for analytic quantum gradients when the number of parameters to be + optimized is large. + + To use this device, you will need to install Autograd: + + .. code-block:: console + + pip install autograd + + **Example** + + The ``default.qubit.autograd`` is designed to be used with end-to-end classical backpropagation + (``diff_method="backprop"``) with the Autograd interface. This is the default method of + differentiation when creating a QNode with this device. + + Using this method, the created QNode is a 'white-box', and is + tightly integrated with your Autograd computation: + + >>> dev = qml.device("default.qubit.autograd", wires=1) + >>> @qml.qnode(dev, interface="autograd", diff_method="backprop") + ... def circuit(x): + ... qml.RX(x[1], wires=0) + ... qml.Rot(x[0], x[1], x[2], wires=0) + ... return qml.expval(qml.PauliZ(0)) + >>> weights = np.array([0.2, 0.5, 0.1]) + >>> grad_fn = qml.grad(circuit) + >>> print(grad_fn(weights)) + array([-2.2526717e-01 -1.0086454e+00 1.3877788e-17]) + + There are a couple of things to keep in mind when using the ``"backprop"`` + differentiation method for QNodes: + + * You must use the ``"autograd"`` interface for classical backpropagation, as Autograd is + used as the device backend. + + * Only exact expectation values, variances, and probabilities are differentiable. + When instantiating the device with ``analytic=False``, differentiating QNode + outputs will result in an error. + + Args: + wires (int): the number of wires to initialize the device with + shots (None, int): How many times the circuit should be evaluated (or sampled) to estimate + the expectation values. Defaults to ``None`` if not specified, which means that the device + returns analytical results. + analytic (bool): Indicates if the device should calculate expectations + and variances analytically. In non-analytic mode, the ``diff_method="backprop"`` + QNode differentiation method is not supported and it is recommended to consider + switching device to ``default.qubit`` and using ``diff_method="parameter-shift"``. + """ + + name = "Default qubit (Autograd) PennyLane plugin" + short_name = "default.qubit.autograd" + + parametric_ops = { + "PhaseShift": autograd_ops.PhaseShift, + "ControlledPhaseShift": autograd_ops.ControlledPhaseShift, + "RX": autograd_ops.RX, + "RY": autograd_ops.RY, + "RZ": autograd_ops.RZ, + "Rot": autograd_ops.Rot, + "CRX": autograd_ops.CRX, + "CRY": autograd_ops.CRY, + "CRZ": autograd_ops.CRZ, + "CRot": autograd_ops.CRot, + "MultiRZ": autograd_ops.MultiRZ, + "SingleExcitation": autograd_ops.SingleExcitation, + "SingleExcitationPlus": autograd_ops.SingleExcitationPlus, + "SingleExcitationMinus": autograd_ops.SingleExcitationMinus, + "DoubleExcitation": autograd_ops.DoubleExcitation, + "DoubleExcitationPlus": autograd_ops.DoubleExcitationPlus, + "DoubleExcitationMinus": autograd_ops.DoubleExcitationMinus, + } + + C_DTYPE = np.complex128 + R_DTYPE = np.float64 + _dot = staticmethod(np.dot) + _abs = staticmethod(np.abs) + _reduce_sum = staticmethod(lambda array, axes: np.sum(array, axis=tuple(axes))) + _reshape = staticmethod(np.reshape) + _flatten = staticmethod(lambda array: array.flatten()) + _gather = staticmethod(lambda array, indices: array[indices]) + _einsum = staticmethod(np.einsum) + _cast = staticmethod(np.asarray) + _transpose = staticmethod(np.transpose) + _tensordot = staticmethod(np.tensordot) + _conj = staticmethod(np.conj) + _imag = staticmethod(np.imag) + _roll = staticmethod(np.roll) + _stack = staticmethod(np.stack) + + @staticmethod + def _asarray(array, dtype=None): + res = np.asarray(array, dtype=dtype) + + if res.dtype is np.dtype("O"): + return np.hstack(array).flatten().astype(dtype) + + return res + + def __init__(self, wires, *, shots=None): + super().__init__(wires, shots=shots, cache=0) + + # prevent using special apply methods for these gates due to slowdown in Autograd + # implementation + del self._apply_ops["PauliY"] + del self._apply_ops["Hadamard"] + del self._apply_ops["CZ"] + + @classmethod + def capabilities(cls): + capabilities = super().capabilities().copy() + capabilities.update( + passthru_interface="autograd", + supports_reversible_diff=False, + ) + return capabilities + + @staticmethod + def _scatter(indices, array, new_dimensions): + new_array = np.zeros(new_dimensions, dtype=array.dtype.type) + new_array[indices] = array + return new_array + + def _get_unitary_matrix(self, unitary): + """Return the matrix representing a unitary operation. + + Args: + unitary (~.Operation): a PennyLane unitary operation + + Returns: + array[complex]: Returns a 2D matrix representation of + the unitary in the computational basis, or, in the case of a diagonal unitary, + a 1D array representing the matrix diagonal. + """ + op_name = unitary.name.split(".inv")[0] + + if op_name in self.parametric_ops: + if op_name == "MultiRZ": + mat = self.parametric_ops[op_name](*unitary.parameters, len(unitary.wires)) + else: + mat = self.parametric_ops[op_name](*unitary.parameters) + + if unitary.inverse: + mat = self._transpose(self._conj(mat)) + + return mat + + if isinstance(unitary, DiagonalOperation): + return unitary.eigvals + + return unitary.matrix diff --git a/pennylane/devices/default_qubit_jax.py b/pennylane/devices/default_qubit_jax.py index 8953be06cf9..9b0a4a8d54e 100644 --- a/pennylane/devices/default_qubit_jax.py +++ b/pennylane/devices/default_qubit_jax.py @@ -1,264 +1,267 @@ -# Copyright 2018-2020 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. -"""This module contains an jax implementation of the :class:`~.DefaultQubit` -reference plugin. -""" - - -from pennylane.operation import DiagonalOperation -from pennylane.devices import DefaultQubit -from pennylane.devices import jax_ops -import numpy as np - -try: - import jax.numpy as jnp - import jax - -except ImportError as e: # pragma: no cover - raise ImportError("default.qubit.jax device requires installing jax>0.2.0") from e - - -class DefaultQubitJax(DefaultQubit): - """Simulator plugin based on ``"default.qubit"``, written using jax. - - **Short name:** ``default.qubit.jax`` - - This device provides a pure-state qubit simulator written using jax. As a result, it - supports classical backpropagation as a means to compute the gradient. This can be faster than - the parameter-shift rule for analytic quantum gradients when the number of parameters to be - optimized is large. - - To use this device, you will need to install jax: - - .. code-block:: console - - pip install jax jaxlib - - **Example** - - The ``default.qubit.jax`` device is designed to be used with end-to-end classical backpropagation - (``diff_method="backprop"``) with the JAX interface. This is the default method of - differentiation when creating a QNode with this device. - - Using this method, the created QNode is a 'white-box', and is - tightly integrated with your JAX computation: - - >>> dev = qml.device("default.qubit.jax", wires=1) - >>> @qml.qnode(dev, interface="jax", diff_method="backprop") - ... def circuit(x): - ... qml.RX(x[1], wires=0) - ... qml.Rot(x[0], x[1], x[2], wires=0) - ... return qml.expval(qml.PauliZ(0)) - >>> weights = jnp.array([0.2, 0.5, 0.1]) - >>> grad_fn = jax.grad(circuit) - >>> print(grad_fn(weights)) - array([-2.2526717e-01 -1.0086454e+00 1.3877788e-17]) - - There are a couple of things to keep in mind when using the ``"backprop"`` - differentiation method for QNodes: - - * You must use the ``"jax"`` interface for classical backpropagation, as JAX is - used as the device backend. - - .. UsageDetails:: - - JAX does randomness in a special way when compared to NumPy, in that all randomness needs to - be seeded. While we handle this for you automatically in op-by-op mode, when using ``jax.jit``, - the automatically generated seed gets constantant compiled. - - Example: - - .. code-block:: python - - dev = qml.device("default.qubit.jax", wires=1) - - @jax.jit - @qml.qnode(dev, interface="jax", diff_method="backprop") - def circuit(): - qml.Hadamard(0) - return qml.sample(qml.PauliZ(wires=0)) - - a = circuit() - b = circuit() # Bad! b will be the exact same samples as a. - - - To fix this, you should wrap your qnode in another function that takes a PRNGKey, and pass - that in during your device construction. - - .. code-block:: python - - @jax.jit - def keyed_circuit(key): - dev = qml.device("default.qubit.jax", interface="jax", prng_key=key) - @qml.qnode(dev, interface="jax", diff_method="backprop") - def circuit(): - qml.Hadamard(0) - return qml.sample(qml.PauliZ(wires=0)) - return circuit() - - key1 = jax.random.PRNGKey(0) - key2 = jax.random.PRNGKey(1) - a = keyed_circuit(key1) - b = keyed_circuit(key2) # b will be different samples now. - - Check out out the `JAX random documentation `__ - for more information. - - Args: - wires (int): The number of wires to initialize the device with. - shots (None, int): How many times the circuit should be evaluated (or sampled) to estimate - the expectation values. Defaults to ``None`` if not specified, which means that the device - returns analytical results. - analytic (bool): Indicates if the device should calculate expectations - and variances analytically. In non-analytic mode, the ``diff_method="backprop"`` - QNode differentiation method is not supported and it is recommended to consider - switching device to ``default.qubit`` and using ``diff_method="parameter-shift"``. - prng_key (Optional[jax.random.PRNGKey]): An optional ``jax.random.PRNGKey``. This is the key to the - pseudo random number generator. If None, a random key will be generated. - - """ - - name = "Default qubit (jax) PennyLane plugin" - short_name = "default.qubit.jax" - - parametric_ops = { - "PhaseShift": jax_ops.PhaseShift, - "ControlledPhaseShift": jax_ops.ControlledPhaseShift, - "RX": jax_ops.RX, - "RY": jax_ops.RY, - "RZ": jax_ops.RZ, - "Rot": jax_ops.Rot, - "CRX": jax_ops.CRX, - "CRY": jax_ops.CRY, - "CRZ": jax_ops.CRZ, - "MultiRZ": jax_ops.MultiRZ, - "SingleExcitation": jax_ops.SingleExcitation, - "SingleExcitationPlus": jax_ops.SingleExcitationPlus, - "SingleExcitationMinus": jax_ops.SingleExcitationMinus, - } - - C_DTYPE = jnp.complex64 - R_DTYPE = jnp.float32 - _asarray = staticmethod(jnp.array) - _dot = staticmethod(jnp.dot) - _abs = staticmethod(jnp.abs) - _reduce_sum = staticmethod(lambda array, axes: jnp.sum(array, axis=tuple(axes))) - _reshape = staticmethod(jnp.reshape) - _flatten = staticmethod(lambda array: array.ravel()) - _gather = staticmethod(lambda array, indices: array[indices]) - _einsum = staticmethod(jnp.einsum) - _cast = staticmethod(jnp.array) - _transpose = staticmethod(jnp.transpose) - _tensordot = staticmethod( - lambda a, b, axes: jnp.tensordot( - a, b, axes if isinstance(axes, int) else list(map(tuple, axes)) - ) - ) - _conj = staticmethod(jnp.conj) - _imag = staticmethod(jnp.imag) - _roll = staticmethod(jnp.roll) - _stack = staticmethod(jnp.stack) - - def __init__(self, wires, *, shots=None, prng_key=None): - super().__init__(wires, shots=shots, cache=0) - - # prevent using special apply methods for these gates due to slowdown in jax - # implementation - del self._apply_ops["PauliY"] - del self._apply_ops["Hadamard"] - del self._apply_ops["CZ"] - self._prng_key = prng_key - - @classmethod - def capabilities(cls): - capabilities = super().capabilities().copy() - capabilities.update( - passthru_interface="jax", - supports_reversible_diff=False, - ) - return capabilities - - @staticmethod - def _scatter(indices, array, new_dimensions): - new_array = jnp.zeros(new_dimensions, dtype=array.dtype.type) - new_array = new_array.at[indices].set(array) - return new_array - - def _get_unitary_matrix(self, unitary): - """Return the matrix representing a unitary operation. - - Args: - unitary (~.Operation): a PennyLane unitary operation - - Returns: - array[complex]: Returns a 2D matrix representation of - the unitary in the computational basis, or, in the case of a diagonal unitary, - a 1D array representing the matrix diagonal. - """ - op_name = unitary.name.split(".inv")[0] - - if op_name in self.parametric_ops: - if op_name == "MultiRZ": - mat = self.parametric_ops[op_name](*unitary.parameters, len(unitary.wires)) - else: - mat = self.parametric_ops[op_name](*unitary.parameters) - - if unitary.inverse: - mat = self._transpose(self._conj(mat)) - - return mat - - if isinstance(unitary, DiagonalOperation): - return unitary.eigvals - - return unitary.matrix - - def sample_basis_states(self, number_of_states, state_probability): - """Sample from the computational basis states based on the state - probability. - - This is an auxiliary method to the generate_samples method. - - Args: - number_of_states (int): the number of basis states to sample from - - Returns: - List[int]: the sampled basis states - """ - if self._prng_key is None: - # Assuming op-by-op, so we'll just make one. - key = jax.random.PRNGKey(np.random.randint(0, 2 ** 31)) - else: - key = self._prng_key - return jax.random.choice(key, number_of_states, shape=(self.shots,), p=state_probability) - - @staticmethod - def states_to_binary(samples, num_wires, dtype=jnp.int32): - """Convert basis states from base 10 to binary representation. - - This is an auxiliary method to the generate_samples method. - - Args: - samples (List[int]): samples of basis states in base 10 representation - num_wires (int): the number of qubits - dtype (type): Type of the internal integer array to be used. Can be - important to specify for large systems for memory allocation - purposes. - - Returns: - List[int]: basis states in binary representation - """ - powers_of_two = 1 << jnp.arange(num_wires, dtype=dtype) - states_sampled_base_ten = samples[:, None] & powers_of_two - return (states_sampled_base_ten > 0).astype(dtype)[:, ::-1] +# Copyright 2018-2020 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. +"""This module contains an jax implementation of the :class:`~.DefaultQubit` +reference plugin. +""" + + +from pennylane.operation import DiagonalOperation +from pennylane.devices import DefaultQubit +from pennylane.devices import jax_ops +import numpy as np + +try: + import jax.numpy as jnp + import jax + +except ImportError as e: # pragma: no cover + raise ImportError("default.qubit.jax device requires installing jax>0.2.0") from e + + +class DefaultQubitJax(DefaultQubit): + """Simulator plugin based on ``"default.qubit"``, written using jax. + + **Short name:** ``default.qubit.jax`` + + This device provides a pure-state qubit simulator written using jax. As a result, it + supports classical backpropagation as a means to compute the gradient. This can be faster than + the parameter-shift rule for analytic quantum gradients when the number of parameters to be + optimized is large. + + To use this device, you will need to install jax: + + .. code-block:: console + + pip install jax jaxlib + + **Example** + + The ``default.qubit.jax`` device is designed to be used with end-to-end classical backpropagation + (``diff_method="backprop"``) with the JAX interface. This is the default method of + differentiation when creating a QNode with this device. + + Using this method, the created QNode is a 'white-box', and is + tightly integrated with your JAX computation: + + >>> dev = qml.device("default.qubit.jax", wires=1) + >>> @qml.qnode(dev, interface="jax", diff_method="backprop") + ... def circuit(x): + ... qml.RX(x[1], wires=0) + ... qml.Rot(x[0], x[1], x[2], wires=0) + ... return qml.expval(qml.PauliZ(0)) + >>> weights = jnp.array([0.2, 0.5, 0.1]) + >>> grad_fn = jax.grad(circuit) + >>> print(grad_fn(weights)) + array([-2.2526717e-01 -1.0086454e+00 1.3877788e-17]) + + There are a couple of things to keep in mind when using the ``"backprop"`` + differentiation method for QNodes: + + * You must use the ``"jax"`` interface for classical backpropagation, as JAX is + used as the device backend. + + .. UsageDetails:: + + JAX does randomness in a special way when compared to NumPy, in that all randomness needs to + be seeded. While we handle this for you automatically in op-by-op mode, when using ``jax.jit``, + the automatically generated seed gets constantant compiled. + + Example: + + .. code-block:: python + + dev = qml.device("default.qubit.jax", wires=1) + + @jax.jit + @qml.qnode(dev, interface="jax", diff_method="backprop") + def circuit(): + qml.Hadamard(0) + return qml.sample(qml.PauliZ(wires=0)) + + a = circuit() + b = circuit() # Bad! b will be the exact same samples as a. + + + To fix this, you should wrap your qnode in another function that takes a PRNGKey, and pass + that in during your device construction. + + .. code-block:: python + + @jax.jit + def keyed_circuit(key): + dev = qml.device("default.qubit.jax", interface="jax", prng_key=key) + @qml.qnode(dev, interface="jax", diff_method="backprop") + def circuit(): + qml.Hadamard(0) + return qml.sample(qml.PauliZ(wires=0)) + return circuit() + + key1 = jax.random.PRNGKey(0) + key2 = jax.random.PRNGKey(1) + a = keyed_circuit(key1) + b = keyed_circuit(key2) # b will be different samples now. + + Check out out the `JAX random documentation `__ + for more information. + + Args: + wires (int): The number of wires to initialize the device with. + shots (None, int): How many times the circuit should be evaluated (or sampled) to estimate + the expectation values. Defaults to ``None`` if not specified, which means that the device + returns analytical results. + analytic (bool): Indicates if the device should calculate expectations + and variances analytically. In non-analytic mode, the ``diff_method="backprop"`` + QNode differentiation method is not supported and it is recommended to consider + switching device to ``default.qubit`` and using ``diff_method="parameter-shift"``. + prng_key (Optional[jax.random.PRNGKey]): An optional ``jax.random.PRNGKey``. This is the key to the + pseudo random number generator. If None, a random key will be generated. + + """ + + name = "Default qubit (jax) PennyLane plugin" + short_name = "default.qubit.jax" + + parametric_ops = { + "PhaseShift": jax_ops.PhaseShift, + "ControlledPhaseShift": jax_ops.ControlledPhaseShift, + "RX": jax_ops.RX, + "RY": jax_ops.RY, + "RZ": jax_ops.RZ, + "Rot": jax_ops.Rot, + "CRX": jax_ops.CRX, + "CRY": jax_ops.CRY, + "CRZ": jax_ops.CRZ, + "MultiRZ": jax_ops.MultiRZ, + "SingleExcitation": jax_ops.SingleExcitation, + "SingleExcitationPlus": jax_ops.SingleExcitationPlus, + "SingleExcitationMinus": jax_ops.SingleExcitationMinus, + "DoubleExcitation": jax_ops.DoubleExcitation, + "DoubleExcitationPlus": jax_ops.DoubleExcitationPlus, + "DoubleExcitationMinus": jax_ops.DoubleExcitationMinus, + } + + C_DTYPE = jnp.complex64 + R_DTYPE = jnp.float32 + _asarray = staticmethod(jnp.array) + _dot = staticmethod(jnp.dot) + _abs = staticmethod(jnp.abs) + _reduce_sum = staticmethod(lambda array, axes: jnp.sum(array, axis=tuple(axes))) + _reshape = staticmethod(jnp.reshape) + _flatten = staticmethod(lambda array: array.ravel()) + _gather = staticmethod(lambda array, indices: array[indices]) + _einsum = staticmethod(jnp.einsum) + _cast = staticmethod(jnp.array) + _transpose = staticmethod(jnp.transpose) + _tensordot = staticmethod( + lambda a, b, axes: jnp.tensordot( + a, b, axes if isinstance(axes, int) else list(map(tuple, axes)) + ) + ) + _conj = staticmethod(jnp.conj) + _imag = staticmethod(jnp.imag) + _roll = staticmethod(jnp.roll) + _stack = staticmethod(jnp.stack) + + def __init__(self, wires, *, shots=None, prng_key=None): + super().__init__(wires, shots=shots, cache=0) + + # prevent using special apply methods for these gates due to slowdown in jax + # implementation + del self._apply_ops["PauliY"] + del self._apply_ops["Hadamard"] + del self._apply_ops["CZ"] + self._prng_key = prng_key + + @classmethod + def capabilities(cls): + capabilities = super().capabilities().copy() + capabilities.update( + passthru_interface="jax", + supports_reversible_diff=False, + ) + return capabilities + + @staticmethod + def _scatter(indices, array, new_dimensions): + new_array = jnp.zeros(new_dimensions, dtype=array.dtype.type) + new_array = new_array.at[indices].set(array) + return new_array + + def _get_unitary_matrix(self, unitary): + """Return the matrix representing a unitary operation. + + Args: + unitary (~.Operation): a PennyLane unitary operation + + Returns: + array[complex]: Returns a 2D matrix representation of + the unitary in the computational basis, or, in the case of a diagonal unitary, + a 1D array representing the matrix diagonal. + """ + op_name = unitary.name.split(".inv")[0] + + if op_name in self.parametric_ops: + if op_name == "MultiRZ": + mat = self.parametric_ops[op_name](*unitary.parameters, len(unitary.wires)) + else: + mat = self.parametric_ops[op_name](*unitary.parameters) + + if unitary.inverse: + mat = self._transpose(self._conj(mat)) + + return mat + + if isinstance(unitary, DiagonalOperation): + return unitary.eigvals + + return unitary.matrix + + def sample_basis_states(self, number_of_states, state_probability): + """Sample from the computational basis states based on the state + probability. + + This is an auxiliary method to the generate_samples method. + + Args: + number_of_states (int): the number of basis states to sample from + + Returns: + List[int]: the sampled basis states + """ + if self._prng_key is None: + # Assuming op-by-op, so we'll just make one. + key = jax.random.PRNGKey(np.random.randint(0, 2 ** 31)) + else: + key = self._prng_key + return jax.random.choice(key, number_of_states, shape=(self.shots,), p=state_probability) + + @staticmethod + def states_to_binary(samples, num_wires, dtype=jnp.int32): + """Convert basis states from base 10 to binary representation. + + This is an auxiliary method to the generate_samples method. + + Args: + samples (List[int]): samples of basis states in base 10 representation + num_wires (int): the number of qubits + dtype (type): Type of the internal integer array to be used. Can be + important to specify for large systems for memory allocation + purposes. + + Returns: + List[int]: basis states in binary representation + """ + powers_of_two = 1 << jnp.arange(num_wires, dtype=dtype) + states_sampled_base_ten = samples[:, None] & powers_of_two + return (states_sampled_base_ten > 0).astype(dtype)[:, ::-1] diff --git a/pennylane/devices/default_qubit_tf.py b/pennylane/devices/default_qubit_tf.py index 879e947cd6d..5551f457c31 100644 --- a/pennylane/devices/default_qubit_tf.py +++ b/pennylane/devices/default_qubit_tf.py @@ -1,234 +1,237 @@ -# Copyright 2018-2020 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. -"""This module contains a TensorFlow implementation of the :class:`~.DefaultQubit` -reference plugin. -""" -import numpy as np -import semantic_version - -from pennylane.operation import DiagonalOperation - -try: - import tensorflow as tf - - if tf.__version__[0] == "1": - raise ImportError("default.qubit.tf device requires TensorFlow>=2.0") - - from tensorflow.python.framework.errors_impl import InvalidArgumentError - - SUPPORTS_APPLY_OPS = semantic_version.match(">=2.3.0", tf.__version__) - -except ImportError as e: - raise ImportError("default.qubit.tf device requires TensorFlow>=2.0") from e - - -# With TF 2.1+, the legacy tf.einsum was renamed to _einsum_v1, while -# the replacement tf.einsum introduced the bug. This try-except block -# will dynamically patch TensorFlow versions where _einsum_v1 exists, to make it the -# default einsum implementation. -# -# For more details, see https://github.com/tensorflow/tensorflow/issues/37307 -try: - from tensorflow.python.ops.special_math_ops import _einsum_v1 - - tf.einsum = _einsum_v1 -except ImportError: - pass - -from . import DefaultQubit -from . import tf_ops - - -class DefaultQubitTF(DefaultQubit): - """Simulator plugin based on ``"default.qubit"``, written using TensorFlow. - - **Short name:** ``default.qubit.tf`` - - This device provides a pure-state qubit simulator written using TensorFlow. - As a result, it supports classical backpropagation as a means to compute the Jacobian. This can - be faster than the parameter-shift rule for analytic quantum gradients - when the number of parameters to be optimized is large. - - To use this device, you will need to install TensorFlow: - - .. code-block:: console - - pip install tensorflow>=2.0 - - **Example** - - The ``default.qubit.tf`` is designed to be used with end-to-end classical backpropagation - (``diff_method="backprop"``) with the TensorFlow interface. This is the default method - of differentiation when creating a QNode with this device. - - Using this method, the created QNode is a 'white-box', and is - tightly integrated with your TensorFlow computation: - - >>> dev = qml.device("default.qubit.tf", wires=1) - >>> @qml.qnode(dev, interface="tf", diff_method="backprop") - ... def circuit(x): - ... qml.RX(x[1], wires=0) - ... qml.Rot(x[0], x[1], x[2], wires=0) - ... return qml.expval(qml.PauliZ(0)) - >>> weights = tf.Variable([0.2, 0.5, 0.1]) - >>> with tf.GradientTape() as tape: - ... res = circuit(weights) - >>> print(tape.gradient(res, weights)) - tf.Tensor([-2.2526717e-01 -1.0086454e+00 1.3877788e-17], shape=(3,), dtype=float32) - - Autograph mode will also work when using classical backpropagation: - - >>> @tf.function - ... def cost(weights): - ... return tf.reduce_sum(circuit(weights)**3) - 1 - >>> with tf.GradientTape() as tape: - ... res = cost(weights) - >>> print(tape.gradient(res, weights)) - tf.Tensor([-3.5471588e-01 -1.5882589e+00 3.4694470e-17], shape=(3,), dtype=float32) - - There are a couple of things to keep in mind when using the ``"backprop"`` - differentiation method for QNodes: - - * You must use the ``"tf"`` interface for classical backpropagation, as TensorFlow is - used as the device backend. - - * Only exact expectation values, variances, and probabilities are differentiable. - When instantiating the device with ``analytic=False``, differentiating QNode - outputs will result in ``None``. - - - If you wish to use a different machine-learning interface, or prefer to calculate quantum - gradients using the ``parameter-shift`` or ``finite-diff`` differentiation methods, - consider using the ``default.qubit`` device instead. - - - Args: - wires (int, Iterable[Number, str]): Number of subsystems represented by the device, - or iterable that contains unique labels for the subsystems as numbers (i.e., ``[-1, 0, 2]``) - or strings (``['ancilla', 'q1', 'q2']``). Default 1 if not specified. - shots (None, int): How many times the circuit should be evaluated (or sampled) to estimate - the expectation values. Defaults to ``None`` if not specified, which means - that the device returns analytical results. - If ``shots > 0`` is used, the ``diff_method="backprop"`` - QNode differentiation method is not supported and it is recommended to consider - switching device to ``default.qubit`` and using ``diff_method="parameter-shift"``. - """ - - name = "Default qubit (TensorFlow) PennyLane plugin" - short_name = "default.qubit.tf" - - parametric_ops = { - "PhaseShift": tf_ops.PhaseShift, - "ControlledPhaseShift": tf_ops.ControlledPhaseShift, - "RX": tf_ops.RX, - "RY": tf_ops.RY, - "RZ": tf_ops.RZ, - "Rot": tf_ops.Rot, - "MultiRZ": tf_ops.MultiRZ, - "CRX": tf_ops.CRX, - "CRY": tf_ops.CRY, - "CRZ": tf_ops.CRZ, - "CRot": tf_ops.CRot, - "SingleExcitation": tf_ops.SingleExcitation, - "SingleExcitationPlus": tf_ops.SingleExcitationPlus, - "SingleExcitationMinus": tf_ops.SingleExcitationMinus, - } - - C_DTYPE = tf.complex128 - R_DTYPE = tf.float64 - _asarray = staticmethod(tf.convert_to_tensor) - _dot = staticmethod(lambda x, y: tf.tensordot(x, y, axes=1)) - _abs = staticmethod(tf.abs) - _reduce_sum = staticmethod(tf.reduce_sum) - _reshape = staticmethod(tf.reshape) - _flatten = staticmethod(lambda tensor: tf.reshape(tensor, [-1])) - _gather = staticmethod(tf.gather) - _einsum = staticmethod(tf.einsum) - _cast = staticmethod(tf.cast) - _transpose = staticmethod(tf.transpose) - _tensordot = staticmethod(tf.tensordot) - _conj = staticmethod(tf.math.conj) - _imag = staticmethod(tf.math.imag) - _roll = staticmethod(tf.roll) - _stack = staticmethod(tf.stack) - - @staticmethod - def _asarray(array, dtype=None): - try: - res = tf.convert_to_tensor(array, dtype=dtype) - except InvalidArgumentError: - res = tf.concat([tf.reshape(i, [-1]) for i in array], axis=0) - - if dtype is not None: - res = tf.cast(res, dtype=dtype) - - return res - - def __init__(self, wires, *, shots=None): - super().__init__(wires, shots=shots, cache=0) - - # prevent using special apply method for this gate due to slowdown in TF implementation - del self._apply_ops["CZ"] - - # Versions of TF before 2.3.0 do not support using the special apply methods as they - # raise an error when calculating the gradient. For versions of TF after 2.3.0, - # special apply methods are also not supported when using more than 8 wires due to - # limitations with TF slicing. - if not SUPPORTS_APPLY_OPS or self.num_wires > 8: - self._apply_ops = {} - - @classmethod - def capabilities(cls): - capabilities = super().capabilities().copy() - capabilities.update( - passthru_interface="tf", - supports_reversible_diff=False, - ) - return capabilities - - @staticmethod - def _scatter(indices, array, new_dimensions): - indices = np.expand_dims(indices, 1) - return tf.scatter_nd(indices, array, new_dimensions) - - def _get_unitary_matrix(self, unitary): - """Return the matrix representing a unitary operation. - - Args: - unitary (~.Operation): a PennyLane unitary operation - - Returns: - tf.Tensor[complex] or array[complex]: Returns a 2D matrix representation of - the unitary in the computational basis, or, in the case of a diagonal unitary, - a 1D array representing the matrix diagonal. For non-parametric unitaries, - the return type will be a ``np.ndarray``. For parametric unitaries, a ``tf.Tensor`` - object will be returned. - """ - op_name = unitary.name.split(".inv")[0] - - if op_name in self.parametric_ops: - if op_name == "MultiRZ": - mat = self.parametric_ops[op_name](*unitary.parameters, len(unitary.wires)) - else: - mat = self.parametric_ops[op_name](*unitary.parameters) - - if unitary.inverse: - mat = self._transpose(self._conj(mat)) - - return mat - - if isinstance(unitary, DiagonalOperation): - return unitary.eigvals - - return unitary.matrix +# Copyright 2018-2020 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. +"""This module contains a TensorFlow implementation of the :class:`~.DefaultQubit` +reference plugin. +""" +import numpy as np +import semantic_version + +from pennylane.operation import DiagonalOperation + +try: + import tensorflow as tf + + if tf.__version__[0] == "1": + raise ImportError("default.qubit.tf device requires TensorFlow>=2.0") + + from tensorflow.python.framework.errors_impl import InvalidArgumentError + + SUPPORTS_APPLY_OPS = semantic_version.match(">=2.3.0", tf.__version__) + +except ImportError as e: + raise ImportError("default.qubit.tf device requires TensorFlow>=2.0") from e + + +# With TF 2.1+, the legacy tf.einsum was renamed to _einsum_v1, while +# the replacement tf.einsum introduced the bug. This try-except block +# will dynamically patch TensorFlow versions where _einsum_v1 exists, to make it the +# default einsum implementation. +# +# For more details, see https://github.com/tensorflow/tensorflow/issues/37307 +try: + from tensorflow.python.ops.special_math_ops import _einsum_v1 + + tf.einsum = _einsum_v1 +except ImportError: + pass + +from . import DefaultQubit +from . import tf_ops + + +class DefaultQubitTF(DefaultQubit): + """Simulator plugin based on ``"default.qubit"``, written using TensorFlow. + + **Short name:** ``default.qubit.tf`` + + This device provides a pure-state qubit simulator written using TensorFlow. + As a result, it supports classical backpropagation as a means to compute the Jacobian. This can + be faster than the parameter-shift rule for analytic quantum gradients + when the number of parameters to be optimized is large. + + To use this device, you will need to install TensorFlow: + + .. code-block:: console + + pip install tensorflow>=2.0 + + **Example** + + The ``default.qubit.tf`` is designed to be used with end-to-end classical backpropagation + (``diff_method="backprop"``) with the TensorFlow interface. This is the default method + of differentiation when creating a QNode with this device. + + Using this method, the created QNode is a 'white-box', and is + tightly integrated with your TensorFlow computation: + + >>> dev = qml.device("default.qubit.tf", wires=1) + >>> @qml.qnode(dev, interface="tf", diff_method="backprop") + ... def circuit(x): + ... qml.RX(x[1], wires=0) + ... qml.Rot(x[0], x[1], x[2], wires=0) + ... return qml.expval(qml.PauliZ(0)) + >>> weights = tf.Variable([0.2, 0.5, 0.1]) + >>> with tf.GradientTape() as tape: + ... res = circuit(weights) + >>> print(tape.gradient(res, weights)) + tf.Tensor([-2.2526717e-01 -1.0086454e+00 1.3877788e-17], shape=(3,), dtype=float32) + + Autograph mode will also work when using classical backpropagation: + + >>> @tf.function + ... def cost(weights): + ... return tf.reduce_sum(circuit(weights)**3) - 1 + >>> with tf.GradientTape() as tape: + ... res = cost(weights) + >>> print(tape.gradient(res, weights)) + tf.Tensor([-3.5471588e-01 -1.5882589e+00 3.4694470e-17], shape=(3,), dtype=float32) + + There are a couple of things to keep in mind when using the ``"backprop"`` + differentiation method for QNodes: + + * You must use the ``"tf"`` interface for classical backpropagation, as TensorFlow is + used as the device backend. + + * Only exact expectation values, variances, and probabilities are differentiable. + When instantiating the device with ``analytic=False``, differentiating QNode + outputs will result in ``None``. + + + If you wish to use a different machine-learning interface, or prefer to calculate quantum + gradients using the ``parameter-shift`` or ``finite-diff`` differentiation methods, + consider using the ``default.qubit`` device instead. + + + Args: + wires (int, Iterable[Number, str]): Number of subsystems represented by the device, + or iterable that contains unique labels for the subsystems as numbers (i.e., ``[-1, 0, 2]``) + or strings (``['ancilla', 'q1', 'q2']``). Default 1 if not specified. + shots (None, int): How many times the circuit should be evaluated (or sampled) to estimate + the expectation values. Defaults to ``None`` if not specified, which means + that the device returns analytical results. + If ``shots > 0`` is used, the ``diff_method="backprop"`` + QNode differentiation method is not supported and it is recommended to consider + switching device to ``default.qubit`` and using ``diff_method="parameter-shift"``. + """ + + name = "Default qubit (TensorFlow) PennyLane plugin" + short_name = "default.qubit.tf" + + parametric_ops = { + "PhaseShift": tf_ops.PhaseShift, + "ControlledPhaseShift": tf_ops.ControlledPhaseShift, + "RX": tf_ops.RX, + "RY": tf_ops.RY, + "RZ": tf_ops.RZ, + "Rot": tf_ops.Rot, + "MultiRZ": tf_ops.MultiRZ, + "CRX": tf_ops.CRX, + "CRY": tf_ops.CRY, + "CRZ": tf_ops.CRZ, + "CRot": tf_ops.CRot, + "SingleExcitation": tf_ops.SingleExcitation, + "SingleExcitationPlus": tf_ops.SingleExcitationPlus, + "SingleExcitationMinus": tf_ops.SingleExcitationMinus, + "DoubleExcitation": tf_ops.DoubleExcitation, + "DoubleExcitationPlus": tf_ops.DoubleExcitationPlus, + "DoubleExcitationMinus": tf_ops.DoubleExcitationMinus, + } + + C_DTYPE = tf.complex128 + R_DTYPE = tf.float64 + _asarray = staticmethod(tf.convert_to_tensor) + _dot = staticmethod(lambda x, y: tf.tensordot(x, y, axes=1)) + _abs = staticmethod(tf.abs) + _reduce_sum = staticmethod(tf.reduce_sum) + _reshape = staticmethod(tf.reshape) + _flatten = staticmethod(lambda tensor: tf.reshape(tensor, [-1])) + _gather = staticmethod(tf.gather) + _einsum = staticmethod(tf.einsum) + _cast = staticmethod(tf.cast) + _transpose = staticmethod(tf.transpose) + _tensordot = staticmethod(tf.tensordot) + _conj = staticmethod(tf.math.conj) + _imag = staticmethod(tf.math.imag) + _roll = staticmethod(tf.roll) + _stack = staticmethod(tf.stack) + + @staticmethod + def _asarray(array, dtype=None): + try: + res = tf.convert_to_tensor(array, dtype=dtype) + except InvalidArgumentError: + res = tf.concat([tf.reshape(i, [-1]) for i in array], axis=0) + + if dtype is not None: + res = tf.cast(res, dtype=dtype) + + return res + + def __init__(self, wires, *, shots=None): + super().__init__(wires, shots=shots, cache=0) + + # prevent using special apply method for this gate due to slowdown in TF implementation + del self._apply_ops["CZ"] + + # Versions of TF before 2.3.0 do not support using the special apply methods as they + # raise an error when calculating the gradient. For versions of TF after 2.3.0, + # special apply methods are also not supported when using more than 8 wires due to + # limitations with TF slicing. + if not SUPPORTS_APPLY_OPS or self.num_wires > 8: + self._apply_ops = {} + + @classmethod + def capabilities(cls): + capabilities = super().capabilities().copy() + capabilities.update( + passthru_interface="tf", + supports_reversible_diff=False, + ) + return capabilities + + @staticmethod + def _scatter(indices, array, new_dimensions): + indices = np.expand_dims(indices, 1) + return tf.scatter_nd(indices, array, new_dimensions) + + def _get_unitary_matrix(self, unitary): + """Return the matrix representing a unitary operation. + + Args: + unitary (~.Operation): a PennyLane unitary operation + + Returns: + tf.Tensor[complex] or array[complex]: Returns a 2D matrix representation of + the unitary in the computational basis, or, in the case of a diagonal unitary, + a 1D array representing the matrix diagonal. For non-parametric unitaries, + the return type will be a ``np.ndarray``. For parametric unitaries, a ``tf.Tensor`` + object will be returned. + """ + op_name = unitary.name.split(".inv")[0] + + if op_name in self.parametric_ops: + if op_name == "MultiRZ": + mat = self.parametric_ops[op_name](*unitary.parameters, len(unitary.wires)) + else: + mat = self.parametric_ops[op_name](*unitary.parameters) + + if unitary.inverse: + mat = self._transpose(self._conj(mat)) + + return mat + + if isinstance(unitary, DiagonalOperation): + return unitary.eigvals + + return unitary.matrix diff --git a/pennylane/devices/jax_ops.py b/pennylane/devices/jax_ops.py index 489c413d7a5..beed57f3fd8 100644 --- a/pennylane/devices/jax_ops.py +++ b/pennylane/devices/jax_ops.py @@ -1,226 +1,328 @@ -# Copyright 2018-2020 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. -r""" -Utility functions and numerical implementations of quantum operations for JAX-based devices. - -""" -import jax.numpy as jnp -from pennylane.utils import pauli_eigs - -C_DTYPE = jnp.complex64 # Use lower precision for better speed on JAX. -R_DTYPE = jnp.float32 - -I = jnp.array([[1, 0], [0, 1]], dtype=C_DTYPE) -X = jnp.array([[0, 1], [1, 0]], dtype=C_DTYPE) -Y = jnp.array([[0j, -1j], [1j, 0j]], dtype=C_DTYPE) -Z = jnp.array([[1, 0], [0, -1]], dtype=C_DTYPE) - -II = jnp.eye(4, dtype=C_DTYPE) -ZZ = jnp.array(jnp.kron(Z, Z), dtype=C_DTYPE) - -IX = jnp.array(jnp.kron(I, X), dtype=C_DTYPE) -IY = jnp.array(jnp.kron(I, Y), dtype=C_DTYPE) -IZ = jnp.array(jnp.kron(I, Z), dtype=C_DTYPE) - -ZI = jnp.array(jnp.kron(Z, I), dtype=C_DTYPE) -ZX = jnp.array(jnp.kron(Z, X), dtype=C_DTYPE) -ZY = jnp.array(jnp.kron(Z, Y), dtype=C_DTYPE) - - -def PhaseShift(phi): - r"""One-qubit phase shift. - - Args: - phi (float): phase shift angle - - Returns: - array[complex]: diagonal part of the phase shift matrix - """ - return jnp.array([1.0, jnp.exp(1j * phi)]) - - -def ControlledPhaseShift(phi): - r"""Two-qubit controlled phase shift. - - Args: - phi (float): phase shift angle - - Returns: - array[complex]: diagonal part of the controlled phase shift matrix - """ - return jnp.array([1.0, 1.0, 1.0, jnp.exp(1j * phi)]) - - -def RX(theta): - r"""One-qubit rotation about the x axis. - - Args: - theta (float): rotation angle - - Returns: - array[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_x \theta/2}` - """ - return jnp.cos(theta / 2) * I + 1j * jnp.sin(-theta / 2) * X - - -def RY(theta): - r"""One-qubit rotation about the y axis. - - Args: - theta (float): rotation angle - - Returns: - array[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_y \theta/2}` - """ - return jnp.cos(theta / 2) * I + 1j * jnp.sin(-theta / 2) * Y - - -def RZ(theta): - r"""One-qubit rotation about the z axis. - - Args: - theta (float): rotation angle - - Returns: - array[complex]: the diagonal part of the rotation matrix :math:`e^{-i \sigma_z \theta/2}` - """ - p = jnp.exp(-0.5j * theta) - return jnp.array([p, jnp.conj(p)]) - - -def Rot(a, b, c): - r"""Arbitrary one-qubit rotation using three Euler angles. - - Args: - a,b,c (float): rotation angles - - Returns: - array[complex]: unitary 2x2 rotation matrix ``rz(c) @ ry(b) @ rz(a)`` - """ - return jnp.diag(RZ(c)) @ RY(b) @ jnp.diag(RZ(a)) - - -def CRX(theta): - r"""Two-qubit controlled rotation about the x axis. - - Args: - theta (float): rotation angle - - Returns: - array[complex]: unitary 4x4 rotation matrix - :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_x(\theta)` - """ - return ( - jnp.cos(theta / 4) ** 2 * II - - 1j * jnp.sin(theta / 2) / 2 * IX - + jnp.sin(theta / 4) ** 2 * ZI - + 1j * jnp.sin(theta / 2) / 2 * ZX - ) - - -def CRY(theta): - r"""Two-qubit controlled rotation about the y axis. - - Args: - theta (float): rotation angle - Returns: - array[complex]: unitary 4x4 rotation matrix :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_y(\theta)` - """ - return ( - jnp.cos(theta / 4) ** 2 * II - - 1j * jnp.sin(theta / 2) / 2 * IY - + jnp.sin(theta / 4) ** 2 * ZI - + 1j * jnp.sin(theta / 2) / 2 * ZY - ) - - -def CRZ(theta): - r"""Two-qubit controlled rotation about the z axis. - - Args: - theta (float): rotation angle - Returns: - array[complex]: diagonal part of the 4x4 rotation matrix - :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_z(\theta)` - """ - p = jnp.exp(-0.5j * theta) - return jnp.array([1.0, 1.0, p, jnp.conj(p)]) - - -def CRot(a, b, c): - r"""Arbitrary two-qubit controlled rotation using three Euler angles. - - Args: - a,b,c (float): rotation angles - Returns: - array[complex]: unitary 4x4 rotation matrix - :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R(a,b,c)` - """ - return jnp.diag(CRZ(c)) @ (CRY(b) @ jnp.diag(CRZ(a))) - - -def MultiRZ(theta, n): - r"""Arbitrary multi Z rotation. - - Args: - theta (float): rotation angle :math:`\theta` - wires (Sequence[int] or int): the wires the operation acts on - Returns: - array[complex]: diagonal part of the multi-qubit rotation matrix - """ - return jnp.exp(-1j * theta / 2 * pauli_eigs(n)) - - -def SingleExcitation(phi): - r"""Single excitation rotation. - - Args: - phi (float): rotation angle - - Returns: - jnp.Tensor[float]: Single excitation rotation matrix - """ - c = jnp.cos(phi / 2) - s = jnp.sin(phi / 2) - return jnp.array([[1, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, 1]]) - - -def SingleExcitationPlus(phi): - r"""Single excitation rotation with positive phase-shift outside the rotation subspace. - - Args: - phi (float): rotation angle - - Returns: - jnp.Tensor[complex]: Single excitation rotation matrix with positive phase-shift - """ - c = jnp.cos(phi / 2) - s = jnp.sin(phi / 2) - e = jnp.exp(1j * phi / 2) - return jnp.array([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) - - -def SingleExcitationMinus(phi): - r"""Single excitation rotation with negative phase-shift outside the rotation subspace. - - Args: - phi (float): rotation angle - - Returns: - tf.Tensor[complex]: Single excitation rotation matrix with negative phase-shift - """ - c = jnp.cos(phi / 2) - s = jnp.sin(phi / 2) - e = jnp.exp(-1j * phi / 2) - return jnp.array([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) +# Copyright 2018-2020 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. +r""" +Utility functions and numerical implementations of quantum operations for JAX-based devices. + +""" +import jax.numpy as jnp +from pennylane.utils import pauli_eigs + +C_DTYPE = jnp.complex64 # Use lower precision for better speed on JAX. +R_DTYPE = jnp.float32 + +I = jnp.array([[1, 0], [0, 1]], dtype=C_DTYPE) +X = jnp.array([[0, 1], [1, 0]], dtype=C_DTYPE) +Y = jnp.array([[0j, -1j], [1j, 0j]], dtype=C_DTYPE) +Z = jnp.array([[1, 0], [0, -1]], dtype=C_DTYPE) + +II = jnp.eye(4, dtype=C_DTYPE) +ZZ = jnp.array(jnp.kron(Z, Z), dtype=C_DTYPE) + +IX = jnp.array(jnp.kron(I, X), dtype=C_DTYPE) +IY = jnp.array(jnp.kron(I, Y), dtype=C_DTYPE) +IZ = jnp.array(jnp.kron(I, Z), dtype=C_DTYPE) + +ZI = jnp.array(jnp.kron(Z, I), dtype=C_DTYPE) +ZX = jnp.array(jnp.kron(Z, X), dtype=C_DTYPE) +ZY = jnp.array(jnp.kron(Z, Y), dtype=C_DTYPE) + + +def PhaseShift(phi): + r"""One-qubit phase shift. + + Args: + phi (float): phase shift angle + + Returns: + array[complex]: diagonal part of the phase shift matrix + """ + return jnp.array([1.0, jnp.exp(1j * phi)]) + + +def ControlledPhaseShift(phi): + r"""Two-qubit controlled phase shift. + + Args: + phi (float): phase shift angle + + Returns: + array[complex]: diagonal part of the controlled phase shift matrix + """ + return jnp.array([1.0, 1.0, 1.0, jnp.exp(1j * phi)]) + + +def RX(theta): + r"""One-qubit rotation about the x axis. + + Args: + theta (float): rotation angle + + Returns: + array[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_x \theta/2}` + """ + return jnp.cos(theta / 2) * I + 1j * jnp.sin(-theta / 2) * X + + +def RY(theta): + r"""One-qubit rotation about the y axis. + + Args: + theta (float): rotation angle + + Returns: + array[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_y \theta/2}` + """ + return jnp.cos(theta / 2) * I + 1j * jnp.sin(-theta / 2) * Y + + +def RZ(theta): + r"""One-qubit rotation about the z axis. + + Args: + theta (float): rotation angle + + Returns: + array[complex]: the diagonal part of the rotation matrix :math:`e^{-i \sigma_z \theta/2}` + """ + p = jnp.exp(-0.5j * theta) + return jnp.array([p, jnp.conj(p)]) + + +def Rot(a, b, c): + r"""Arbitrary one-qubit rotation using three Euler angles. + + Args: + a,b,c (float): rotation angles + + Returns: + array[complex]: unitary 2x2 rotation matrix ``rz(c) @ ry(b) @ rz(a)`` + """ + return jnp.diag(RZ(c)) @ RY(b) @ jnp.diag(RZ(a)) + + +def CRX(theta): + r"""Two-qubit controlled rotation about the x axis. + + Args: + theta (float): rotation angle + + Returns: + array[complex]: unitary 4x4 rotation matrix + :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_x(\theta)` + """ + return ( + jnp.cos(theta / 4) ** 2 * II + - 1j * jnp.sin(theta / 2) / 2 * IX + + jnp.sin(theta / 4) ** 2 * ZI + + 1j * jnp.sin(theta / 2) / 2 * ZX + ) + + +def CRY(theta): + r"""Two-qubit controlled rotation about the y axis. + + Args: + theta (float): rotation angle + Returns: + array[complex]: unitary 4x4 rotation matrix :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_y(\theta)` + """ + return ( + jnp.cos(theta / 4) ** 2 * II + - 1j * jnp.sin(theta / 2) / 2 * IY + + jnp.sin(theta / 4) ** 2 * ZI + + 1j * jnp.sin(theta / 2) / 2 * ZY + ) + + +def CRZ(theta): + r"""Two-qubit controlled rotation about the z axis. + + Args: + theta (float): rotation angle + Returns: + array[complex]: diagonal part of the 4x4 rotation matrix + :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_z(\theta)` + """ + p = jnp.exp(-0.5j * theta) + return jnp.array([1.0, 1.0, p, jnp.conj(p)]) + + +def CRot(a, b, c): + r"""Arbitrary two-qubit controlled rotation using three Euler angles. + + Args: + a,b,c (float): rotation angles + Returns: + array[complex]: unitary 4x4 rotation matrix + :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R(a,b,c)` + """ + return jnp.diag(CRZ(c)) @ (CRY(b) @ jnp.diag(CRZ(a))) + + +def MultiRZ(theta, n): + r"""Arbitrary multi Z rotation. + + Args: + theta (float): rotation angle :math:`\theta` + wires (Sequence[int] or int): the wires the operation acts on + Returns: + array[complex]: diagonal part of the multi-qubit rotation matrix + """ + return jnp.exp(-1j * theta / 2 * pauli_eigs(n)) + + +def SingleExcitation(phi): + r"""Single excitation rotation. + + Args: + phi (float): rotation angle + + Returns: + jnp.Tensor[float]: Single excitation rotation matrix + """ + c = jnp.cos(phi / 2) + s = jnp.sin(phi / 2) + return jnp.array([[1, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, 1]]) + + +def SingleExcitationPlus(phi): + r"""Single excitation rotation with positive phase-shift outside the rotation subspace. + + Args: + phi (float): rotation angle + + Returns: + jnp.Tensor[complex]: Single excitation rotation matrix with positive phase-shift + """ + c = jnp.cos(phi / 2) + s = jnp.sin(phi / 2) + e = jnp.exp(1j * phi / 2) + return jnp.array([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) + + +def SingleExcitationMinus(phi): + r"""Single excitation rotation with negative phase-shift outside the rotation subspace. + + Args: + phi (float): rotation angle + + Returns: + tf.Tensor[complex]: Single excitation rotation matrix with negative phase-shift + """ + c = jnp.cos(phi / 2) + s = jnp.sin(phi / 2) + e = jnp.exp(-1j * phi / 2) + return jnp.array([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) + + +def DoubleExcitation(phi): + r"""Double excitation rotation. + + Args: + phi (float): rotation angle + Returns: + jnp.Tensor[float]: Double excitation rotation matrix + + """ + c = jnp.cos(phi / 2) + s = jnp.sin(phi / 2) + + U = [ + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, c, 0, 0, 0, 0, 0, 0, 0, 0, -s, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, s, 0, 0, 0, 0, 0, 0, 0, 0, c, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], + ] + + return jnp.array(U) + + +def DoubleExcitationPlus(phi): + r"""Double excitation rotation with positive phase-shift. + + Args: + phi (float): rotation angle + Returns: + jnp.Tensor[complex]: rotation matrix + """ + c = jnp.cos(phi / 2) + s = jnp.sin(phi / 2) + e = jnp.exp(1j * phi / 2) + + U = [ + [e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, c, 0, 0, 0, 0, 0, 0, 0, 0, -s, 0, 0, 0], + [0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0], + [0, 0, 0, s, 0, 0, 0, 0, 0, 0, 0, 0, c, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e], + ] + + return jnp.array(U) + + +def DoubleExcitationMinus(phi): + r"""Double excitation rotation with negative phase-shift. + + Args: + phi (float): rotation angle + Returns: + jnp.Tensor[complex]: rotation matrix + """ + c = jnp.cos(phi / 2) + s = jnp.sin(phi / 2) + e = jnp.exp(-1j * phi / 2) + + U = [ + [e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, c, 0, 0, 0, 0, 0, 0, 0, 0, -s, 0, 0, 0], + [0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0], + [0, 0, 0, s, 0, 0, 0, 0, 0, 0, 0, 0, c, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e], + ] + + return jnp.array(U) diff --git a/pennylane/devices/tests/test_gates.py b/pennylane/devices/tests/test_gates.py index 242292c48f9..17835b10e1f 100755 --- a/pennylane/devices/tests/test_gates.py +++ b/pennylane/devices/tests/test_gates.py @@ -71,6 +71,9 @@ "SingleExcitation": qml.SingleExcitation(0, wires=[0, 1]), "SingleExcitationPlus": qml.SingleExcitationPlus(0, wires=[0, 1]), "SingleExcitationMinus": qml.SingleExcitationMinus(0, wires=[0, 1]), + "DoubleExcitation": qml.DoubleExcitation(0, wires=[0, 1, 2, 3]), + "DoubleExcitationPlus": qml.DoubleExcitationPlus(0, wires=[0, 1, 2, 3]), + "DoubleExcitationMinus": qml.DoubleExcitationMinus(0, wires=[0, 1, 2, 3]), } all_ops = ops.keys() @@ -201,7 +204,7 @@ class TestSupportedGates: @pytest.mark.parametrize("operation", all_ops) def test_supported_gates_can_be_implemented(self, device_kwargs, operation): """Test that the device can implement all its supported gates.""" - device_kwargs["wires"] = 3 # maximum size of current gates + device_kwargs["wires"] = 4 # maximum size of current gates dev = qml.device(**device_kwargs) assert hasattr(dev, "operations") @@ -218,7 +221,7 @@ def circuit(): def test_inverse_gates_can_be_implemented(self, device_kwargs, operation): """Test that the device can implement the inverse of all its supported gates. This test is skipped for devices that do not support inverse operations.""" - device_kwargs["wires"] = 3 + device_kwargs["wires"] = 4 dev = qml.device(**device_kwargs) supports_inv = ( "supports_inverse_operations" in dev.capabilities() diff --git a/pennylane/devices/tf_ops.py b/pennylane/devices/tf_ops.py index 8ef7afd1c7c..9b569593570 100644 --- a/pennylane/devices/tf_ops.py +++ b/pennylane/devices/tf_ops.py @@ -1,241 +1,346 @@ -# Copyright 2018-2020 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. -r""" -Utility functions and numerical implementations of quantum operations TensorFlow devices. -""" -import tensorflow as tf -from numpy import kron -from pennylane.utils import pauli_eigs - -C_DTYPE = tf.complex128 -R_DTYPE = tf.float64 - -I = tf.constant([[1, 0], [0, 1]], dtype=C_DTYPE) -X = tf.constant([[0, 1], [1, 0]], dtype=C_DTYPE) -Y = tf.constant([[0j, -1j], [1j, 0j]], dtype=C_DTYPE) -Z = tf.constant([[1, 0], [0, -1]], dtype=C_DTYPE) - -II = tf.eye(4, dtype=C_DTYPE) -ZZ = tf.constant(kron(Z, Z), dtype=C_DTYPE) - -IX = tf.constant(kron(I, X), dtype=C_DTYPE) -IY = tf.constant(kron(I, Y), dtype=C_DTYPE) -IZ = tf.constant(kron(I, Z), dtype=C_DTYPE) - -ZI = tf.constant(kron(Z, I), dtype=C_DTYPE) -ZX = tf.constant(kron(Z, X), dtype=C_DTYPE) -ZY = tf.constant(kron(Z, Y), dtype=C_DTYPE) - - -def PhaseShift(phi): - r"""One-qubit phase shift. - - Args: - phi (float): phase shift angle - - Returns: - tf.Tensor[complex]: diagonal part of the phase shift matrix - """ - phi = tf.cast(phi, dtype=C_DTYPE) - return tf.convert_to_tensor([1.0, tf.exp(1j * phi)]) - - -def ControlledPhaseShift(phi): - r"""Two-qubit controlled phase shift. - - Args: - phi (float): phase shift angle - - Returns: - tf.Tensor[complex]: diagonal part of the controlled phase shift matrix - """ - phi = tf.cast(phi, dtype=C_DTYPE) - return tf.convert_to_tensor([1.0, 1.0, 1.0, tf.exp(1j * phi)]) - - -def RX(theta): - r"""One-qubit rotation about the x axis. - - Args: - theta (float): rotation angle - - Returns: - tf.Tensor[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_x \theta/2}` - """ - theta = tf.cast(theta, dtype=C_DTYPE) - return tf.cos(theta / 2) * I + 1j * tf.sin(-theta / 2) * X - - -def RY(theta): - r"""One-qubit rotation about the y axis. - - Args: - theta (float): rotation angle - - Returns: - tf.Tensor[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_y \theta/2}` - """ - theta = tf.cast(theta, dtype=C_DTYPE) - return tf.cos(theta / 2) * I + 1j * tf.sin(-theta / 2) * Y - - -def RZ(theta): - r"""One-qubit rotation about the z axis. - - Args: - theta (float): rotation angle - - Returns: - tf.Tensor[complex]: the diagonal part of the rotation matrix :math:`e^{-i \sigma_z \theta/2}` - """ - theta = tf.cast(theta, dtype=C_DTYPE) - p = tf.exp(-0.5j * theta) - return tf.convert_to_tensor([p, tf.math.conj(p)]) - - -def Rot(a, b, c): - r"""Arbitrary one-qubit rotation using three Euler angles. - - Args: - a,b,c (float): rotation angles - - Returns: - tf.Tensor[complex]: unitary 2x2 rotation matrix ``rz(c) @ ry(b) @ rz(a)`` - """ - return tf.linalg.diag(RZ(c)) @ RY(b) @ tf.linalg.diag(RZ(a)) - - -def MultiRZ(theta, n): - r"""Arbitrary multi Z rotation. - - Args: - theta (float): rotation angle - n (int): number of wires the rotation acts on - - Returns: - tf.Tensor[complex]: diagonal part of the MultiRZ matrix - """ - theta = tf.cast(theta, dtype=C_DTYPE) - multi_Z_rot_eigs = tf.exp(-1j * theta / 2 * pauli_eigs(n)) - return tf.convert_to_tensor(multi_Z_rot_eigs) - - -def CRX(theta): - r"""Two-qubit controlled rotation about the x axis. - - Args: - theta (float): rotation angle - - Returns: - tf.Tensor[complex]: unitary 4x4 rotation matrix - :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_x(\theta)` - """ - theta = tf.cast(theta, dtype=C_DTYPE) - return ( - tf.cos(theta / 4) ** 2 * II - - 1j * tf.sin(theta / 2) / 2 * IX - + tf.sin(theta / 4) ** 2 * ZI - + 1j * tf.sin(theta / 2) / 2 * ZX - ) - - -def CRY(theta): - r"""Two-qubit controlled rotation about the y axis. - - Args: - theta (float): rotation angle - Returns: - tf.Tensor[complex]: unitary 4x4 rotation matrix :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_y(\theta)` - """ - theta = tf.cast(theta, dtype=C_DTYPE) - return ( - tf.cos(theta / 4) ** 2 * II - - 1j * tf.sin(theta / 2) / 2 * IY - + tf.sin(theta / 4) ** 2 * ZI - + 1j * tf.sin(theta / 2) / 2 * ZY - ) - - -def CRZ(theta): - r"""Two-qubit controlled rotation about the z axis. - - Args: - theta (float): rotation angle - Returns: - tf.Tensor[complex]: diagonal part of the 4x4 rotation matrix - :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_z(\theta)` - """ - theta = tf.cast(theta, dtype=C_DTYPE) - p = tf.exp(-0.5j * theta) - return tf.convert_to_tensor([1.0, 1.0, p, tf.math.conj(p)]) - - -def CRot(a, b, c): - r"""Arbitrary two-qubit controlled rotation using three Euler angles. - - Args: - a,b,c (float): rotation angles - Returns: - tf.Tensor[complex]: unitary 4x4 rotation matrix - :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R(a,b,c)` - """ - return tf.linalg.diag(CRZ(c)) @ (CRY(b) @ tf.linalg.diag(CRZ(a))) - - -def SingleExcitation(phi): - r"""Single excitation rotation. - - Args: - phi (float): rotation angle - - Returns: - tf.Tensor[complex]: Single excitation rotation matrix - - """ - phi = tf.cast(phi, dtype=C_DTYPE) - c = tf.cos(phi / 2) - s = tf.sin(phi / 2) - return tf.convert_to_tensor([[1, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, 1]]) - - -def SingleExcitationPlus(phi): - r"""Single excitation rotation with positive phase-shift outside the rotation subspace. - - Args: - phi (float): rotation angle - - Returns: - tf.Tensor[complex]: Single excitation rotation matrix with positive phase-shift - """ - phi = tf.cast(phi, dtype=C_DTYPE) - c = tf.cos(phi / 2) - s = tf.sin(phi / 2) - e = tf.exp(1j * phi / 2) - return tf.convert_to_tensor([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) - - -def SingleExcitationMinus(phi): - r"""Single excitation rotation with negative phase-shift outside the rotation subspace. - - Args: - phi (float): rotation angle - - Returns: - tf.Tensor[complex]: Single excitation rotation matrix with negative phase-shift - """ - phi = tf.cast(phi, dtype=C_DTYPE) - c = tf.cos(phi / 2) - s = tf.sin(phi / 2) - e = tf.exp(-1j * phi / 2) - return tf.convert_to_tensor([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) +# Copyright 2018-2020 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. +r""" +Utility functions and numerical implementations of quantum operations TensorFlow devices. +""" +import tensorflow as tf +from numpy import kron +from pennylane.utils import pauli_eigs + +C_DTYPE = tf.complex128 +R_DTYPE = tf.float64 + +I = tf.constant([[1, 0], [0, 1]], dtype=C_DTYPE) +X = tf.constant([[0, 1], [1, 0]], dtype=C_DTYPE) +Y = tf.constant([[0j, -1j], [1j, 0j]], dtype=C_DTYPE) +Z = tf.constant([[1, 0], [0, -1]], dtype=C_DTYPE) + +II = tf.eye(4, dtype=C_DTYPE) +ZZ = tf.constant(kron(Z, Z), dtype=C_DTYPE) + +IX = tf.constant(kron(I, X), dtype=C_DTYPE) +IY = tf.constant(kron(I, Y), dtype=C_DTYPE) +IZ = tf.constant(kron(I, Z), dtype=C_DTYPE) + +ZI = tf.constant(kron(Z, I), dtype=C_DTYPE) +ZX = tf.constant(kron(Z, X), dtype=C_DTYPE) +ZY = tf.constant(kron(Z, Y), dtype=C_DTYPE) + + +def PhaseShift(phi): + r"""One-qubit phase shift. + + Args: + phi (float): phase shift angle + + Returns: + tf.Tensor[complex]: diagonal part of the phase shift matrix + """ + phi = tf.cast(phi, dtype=C_DTYPE) + return tf.convert_to_tensor([1.0, tf.exp(1j * phi)]) + + +def ControlledPhaseShift(phi): + r"""Two-qubit controlled phase shift. + + Args: + phi (float): phase shift angle + + Returns: + tf.Tensor[complex]: diagonal part of the controlled phase shift matrix + """ + phi = tf.cast(phi, dtype=C_DTYPE) + return tf.convert_to_tensor([1.0, 1.0, 1.0, tf.exp(1j * phi)]) + + +def RX(theta): + r"""One-qubit rotation about the x axis. + + Args: + theta (float): rotation angle + + Returns: + tf.Tensor[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_x \theta/2}` + """ + theta = tf.cast(theta, dtype=C_DTYPE) + return tf.cos(theta / 2) * I + 1j * tf.sin(-theta / 2) * X + + +def RY(theta): + r"""One-qubit rotation about the y axis. + + Args: + theta (float): rotation angle + + Returns: + tf.Tensor[complex]: unitary 2x2 rotation matrix :math:`e^{-i \sigma_y \theta/2}` + """ + theta = tf.cast(theta, dtype=C_DTYPE) + return tf.cos(theta / 2) * I + 1j * tf.sin(-theta / 2) * Y + + +def RZ(theta): + r"""One-qubit rotation about the z axis. + + Args: + theta (float): rotation angle + + Returns: + tf.Tensor[complex]: the diagonal part of the rotation matrix :math:`e^{-i \sigma_z \theta/2}` + """ + theta = tf.cast(theta, dtype=C_DTYPE) + p = tf.exp(-0.5j * theta) + return tf.convert_to_tensor([p, tf.math.conj(p)]) + + +def Rot(a, b, c): + r"""Arbitrary one-qubit rotation using three Euler angles. + + Args: + a,b,c (float): rotation angles + + Returns: + tf.Tensor[complex]: unitary 2x2 rotation matrix ``rz(c) @ ry(b) @ rz(a)`` + """ + return tf.linalg.diag(RZ(c)) @ RY(b) @ tf.linalg.diag(RZ(a)) + + +def MultiRZ(theta, n): + r"""Arbitrary multi Z rotation. + + Args: + theta (float): rotation angle + n (int): number of wires the rotation acts on + + Returns: + tf.Tensor[complex]: diagonal part of the MultiRZ matrix + """ + theta = tf.cast(theta, dtype=C_DTYPE) + multi_Z_rot_eigs = tf.exp(-1j * theta / 2 * pauli_eigs(n)) + return tf.convert_to_tensor(multi_Z_rot_eigs) + + +def CRX(theta): + r"""Two-qubit controlled rotation about the x axis. + + Args: + theta (float): rotation angle + + Returns: + tf.Tensor[complex]: unitary 4x4 rotation matrix + :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_x(\theta)` + """ + theta = tf.cast(theta, dtype=C_DTYPE) + return ( + tf.cos(theta / 4) ** 2 * II + - 1j * tf.sin(theta / 2) / 2 * IX + + tf.sin(theta / 4) ** 2 * ZI + + 1j * tf.sin(theta / 2) / 2 * ZX + ) + + +def CRY(theta): + r"""Two-qubit controlled rotation about the y axis. + + Args: + theta (float): rotation angle + Returns: + tf.Tensor[complex]: unitary 4x4 rotation matrix :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_y(\theta)` + """ + theta = tf.cast(theta, dtype=C_DTYPE) + return ( + tf.cos(theta / 4) ** 2 * II + - 1j * tf.sin(theta / 2) / 2 * IY + + tf.sin(theta / 4) ** 2 * ZI + + 1j * tf.sin(theta / 2) / 2 * ZY + ) + + +def CRZ(theta): + r"""Two-qubit controlled rotation about the z axis. + + Args: + theta (float): rotation angle + Returns: + tf.Tensor[complex]: diagonal part of the 4x4 rotation matrix + :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R_z(\theta)` + """ + theta = tf.cast(theta, dtype=C_DTYPE) + p = tf.exp(-0.5j * theta) + return tf.convert_to_tensor([1.0, 1.0, p, tf.math.conj(p)]) + + +def CRot(a, b, c): + r"""Arbitrary two-qubit controlled rotation using three Euler angles. + + Args: + a,b,c (float): rotation angles + Returns: + tf.Tensor[complex]: unitary 4x4 rotation matrix + :math:`|0\rangle\langle 0|\otimes \mathbb{I}+|1\rangle\langle 1|\otimes R(a,b,c)` + """ + return tf.linalg.diag(CRZ(c)) @ (CRY(b) @ tf.linalg.diag(CRZ(a))) + + +def SingleExcitation(phi): + r"""Single excitation rotation. + + Args: + phi (float): rotation angle + + Returns: + tf.Tensor[complex]: Single excitation rotation matrix + """ + phi = tf.cast(phi, dtype=C_DTYPE) + c = tf.cos(phi / 2) + s = tf.sin(phi / 2) + + return tf.convert_to_tensor([[1, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, 1]]) + + +def SingleExcitationPlus(phi): + r"""Single excitation rotation with positive phase-shift outside the rotation subspace. + + Args: + phi (float): rotation angle + + Returns: + tf.Tensor[complex]: Single excitation rotation matrix with positive phase-shift + """ + phi = tf.cast(phi, dtype=C_DTYPE) + c = tf.cos(phi / 2) + s = tf.sin(phi / 2) + e = tf.exp(1j * phi / 2) + return tf.convert_to_tensor([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) + + +def SingleExcitationMinus(phi): + r"""Single excitation rotation with negative phase-shift outside the rotation subspace. + + Args: + phi (float): rotation angle + + Returns: + tf.Tensor[complex]: Single excitation rotation matrix with negative phase-shift + """ + phi = tf.cast(phi, dtype=C_DTYPE) + c = tf.cos(phi / 2) + s = tf.sin(phi / 2) + e = tf.exp(-1j * phi / 2) + return tf.convert_to_tensor([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]]) + + +def DoubleExcitation(phi): + r"""Double excitation rotation. + + Args: + phi (float): rotation angle + Returns: + tf.Tensor[complex]: Double excitation rotation matrix + """ + + phi = tf.cast(phi, dtype=C_DTYPE) + c = tf.cos(phi / 2) + s = tf.sin(phi / 2) + + U = [ + [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, c, 0, 0, 0, 0, 0, 0, 0, 0, -s, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, s, 0, 0, 0, 0, 0, 0, 0, 0, c, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], + ] + + return tf.convert_to_tensor(U) + + +def DoubleExcitationPlus(phi): + r"""Double excitation rotation with positive phase-shift. + + Args: + phi (float): rotation angle + Returns: + tf.Tensor[complex]: rotation matrix + """ + phi = tf.cast(phi, dtype=C_DTYPE) + c = tf.cos(phi / 2) + s = tf.sin(phi / 2) + e = tf.exp(1j * phi / 2) + + U = [ + [e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, c, 0, 0, 0, 0, 0, 0, 0, 0, -s, 0, 0, 0], + [0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0], + [0, 0, 0, s, 0, 0, 0, 0, 0, 0, 0, 0, c, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e], + ] + + return tf.convert_to_tensor(U) + + +def DoubleExcitationMinus(phi): + r"""Double excitation rotation with negative phase-shift. + + Args: + phi (float): rotation angle + Returns: + tf.Tensor[complex]: rotation matrix + """ + phi = tf.cast(phi, dtype=C_DTYPE) + c = tf.cos(phi / 2) + s = tf.sin(phi / 2) + e = tf.exp(-1j * phi / 2) + + U = [ + [e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, c, 0, 0, 0, 0, 0, 0, 0, 0, -s, 0, 0, 0], + [0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0, 0, 0], + [0, 0, 0, s, 0, 0, 0, 0, 0, 0, 0, 0, c, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, e], + ] + + return tf.convert_to_tensor(U) diff --git a/pennylane/ops/qubit.py b/pennylane/ops/qubit.py index b65e38ac738..1ed673519ea 100644 --- a/pennylane/ops/qubit.py +++ b/pennylane/ops/qubit.py @@ -2113,6 +2113,199 @@ def adjoint(self, do_queue=False): return QFT(wires=self.wires, do_queue=do_queue).inv() +# ============================================================================= +# Quantum chemistry +# ============================================================================= + + +class DoubleExcitation(Operation): + r"""DoubleExcitation(phi, wires) + Double excitation rotation. + + This operation performs an :math:`SO(2)` rotation in the two-dimensional subspace :math:`\{ + |1100\rangle,|0011\rangle\}`. More precisely, it performs the transformation + + .. math:: + + &|0011\rangle \rightarrow \cos(\phi) |0011\rangle - \sin(\phi) |1100\rangle\\ + &|1100\rangle \rightarrow \cos(\phi) |1100\rangle + \sin(\phi) |0011\rangle, + + while leaving all other basis states unchanged. + + The name originates from the occupation-number representation of fermionic wavefunctions, where + the transformation from :math:`|1100\rangle` to :math:`|0011\rangle` is interpreted as + "exciting" two particles from the first pair of qubits to the second pair of qubits. + + **Details:** + + * Number of wires: 4 + * Number of parameters: 1 + * Gradient recipe: Obtained from its decomposition in terms of the + :class:`~.DoubleExcitationPlus` and :class:`~.DoubleExcitationMinus` operations + + Args: + phi (float): rotation angle :math:`\phi` + wires (Sequence[int]): the wires the operation acts on + + **Example** + + The following circuit performs the transformation :math:`|1100\rangle\rightarrow \cos( + \phi/2)|1100\rangle -\sin(\phi/2)|0011\rangle)`: + + .. code-block:: + + @qml.qnode(dev) + def circuit(phi): + qml.PauliX(wires=0) + qml.PauliX(wires=1) + qml.DoubleExcitation(phi, wires=[0, 1, 2, 3]) + """ + + num_params = 1 + num_wires = 4 + par_domain = "R" + grad_method = "A" + + G = np.zeros((16, 16), dtype=np.complex64) + G[3, 12] = -1j # 3 (dec) = 0011 (bin) + G[12, 3] = 1j # 12 (dec) = 1100 (bin) + generator = [G, -1 / 2] + + @classmethod + def _matrix(cls, *params): + theta = params[0] + c = math.cos(theta / 2) + s = math.sin(theta / 2) + + U = np.eye(16) + U[3, 3] = c # 3 (dec) = 0011 (bin) + U[3, 12] = -s # 12 (dec) = 1100 (bin) + U[12, 3] = s + U[12, 12] = c + + return U + + @staticmethod + def decomposition(theta, wires): + decomp_ops = [ + DoubleExcitationPlus(theta / 2, wires=wires), + DoubleExcitationMinus(theta / 2, wires=wires), + ] + return decomp_ops + + +class DoubleExcitationPlus(Operation): + r"""DoubleExcitationPlus(phi, wires) + Double excitation rotation with positive phase-shift outside the rotation subspace. + + This operation performs an :math:`SO(2)` rotation in the two-dimensional subspace :math:`\{ + |1100\rangle,|0011\rangle\}` while applying a phase-shift on other states. More precisely, + it performs the transformation + + .. math:: + + &|0011\rangle \rightarrow \cos(\phi) |0011\rangle - \sin(\phi) |1100\rangle\\ + &|1100\rangle \rightarrow \cos(\phi) |1100\rangle + \sin(\phi) |0011\rangle\\ + &|x\rangle \rightarrow e^{i\phi} |x\rangle, + + for all other basis states :math:`|x\rangle`. + + **Details:** + + * Number of wires: 4 + * Number of parameters: 1 + * Gradient recipe: :math:`\frac{d}{d\phi}f(U_+(\phi)) = \frac{1}{2}\left[f(U_+(\phi+\pi/2)) - f(U_+(\phi-\pi/2))\right]` + where :math:`f` is an expectation value depending on :math:`U_+(\phi)` + + Args: + phi (float): rotation angle :math:`\phi` + wires (Sequence[int]): the wires the operation acts on + """ + + num_params = 1 + num_wires = 4 + par_domain = "R" + grad_method = "A" + + G = -1 * np.eye(16, dtype=np.complex64) + G[3, 3] = 0 + G[12, 12] = 0 + G[3, 12] = -1j # 3 (dec) = 0011 (bin) + G[12, 3] = 1j # 12 (dec) = 1100 (bin) + generator = [G, -1 / 2] + + @classmethod + def _matrix(cls, *params): + theta = params[0] + c = math.cos(theta / 2) + s = math.sin(theta / 2) + e = cmath.exp(1j * theta / 2) + + U = e * np.eye(16, dtype=np.complex64) + U[3, 3] = c # 3 (dec) = 0011 (bin) + U[3, 12] = -s # 12 (dec) = 1100 (bin) + U[12, 3] = s + U[12, 12] = c + + return U + + +class DoubleExcitationMinus(Operation): + r"""DoubleExcitationMinus(phi, wires) + Double excitation rotation with negative phase-shift outside the rotation subspace. + + This operation performs an :math:`SO(2)` rotation in the two-dimensional subspace :math:`\{ + |1100\rangle,|0011\rangle\}` while applying a phase-shift on other states. More precisely, + it performs the transformation + + .. math:: + + &|0011\rangle \rightarrow \cos(\phi) |0011\rangle - \sin(\phi) |1100\rangle\\ + &|1100\rangle \rightarrow \cos(\phi) |1100\rangle + \sin(\phi) |0011\rangle\\ + &|x\rangle \rightarrow e^{-i\phi} |x\rangle, + + for all other basis states :math:`|x\rangle`. + + **Details:** + + * Number of wires: 4 + * Number of parameters: 1 + * Gradient recipe: :math:`\frac{d}{d\phi}f(U_-(\phi)) = \frac{1}{2}\left[f(U_-(\phi+\pi/2)) - f(U_-(\phi-\pi/2))\right]` + where :math:`f` is an expectation value depending on :math:`U_-(\phi)` + + Args: + phi (float): rotation angle :math:`\phi` + wires (Sequence[int]): the wires the operation acts on + """ + + num_params = 1 + num_wires = 4 + par_domain = "R" + grad_method = "A" + + G = np.eye(16, dtype=np.complex64) + G[3, 3] = 0 + G[12, 12] = 0 + G[3, 12] = -1j # 3 (dec) = 0011 (bin) + G[12, 3] = 1j # 12 (dec) = 1100 (bin) + generator = [G, -1 / 2] + + @classmethod + def _matrix(cls, *params): + theta = params[0] + c = math.cos(theta / 2) + s = math.sin(theta / 2) + e = cmath.exp(-1j * theta / 2) + + U = e * np.eye(16, dtype=np.complex64) + U[3, 3] = c # 3 (dec) = 0011 (bin) + U[3, 12] = -s # 12 (dec) = 1100 (bin) + U[12, 3] = s + U[12, 12] = c + + return U + + # ============================================================================= # State preparation # ============================================================================= @@ -2318,6 +2511,9 @@ def diagonalizing_gates(self): "SingleExcitation", "SingleExcitationPlus", "SingleExcitationMinus", + "DoubleExcitation", + "DoubleExcitationPlus", + "DoubleExcitationMinus", } diff --git a/tests/gate_data.py b/tests/gate_data.py index 44836637da6..e5b24c2bfbb 100644 --- a/tests/gate_data.py +++ b/tests/gate_data.py @@ -233,8 +233,8 @@ def ControlledPhaseShift(phi): array: the two-wire controlled-phase matrix """ return np.diag([1, 1, 1, np.exp(1j * phi)]) - - + + def SingleExcitation(phi): r"""Single excitation rotation. @@ -251,12 +251,12 @@ def SingleExcitation(phi): def SingleExcitationPlus(phi): r"""Single excitation rotation with positive phase shift. - + Args: - phi (float): rotation angle + phi (float): rotation angle Returns: - array: the two-qubit matrix describing the operation + array: the two-qubit Givens rotation describing the single excitation operation """ return np.array([[np.exp(1j*phi/2), 0, 0, 0], [0, np.cos(phi/2), -np.sin(phi/2), 0], @@ -275,3 +275,69 @@ def SingleExcitationMinus(phi): return np.array([[np.exp(-1j*phi/2), 0, 0, 0], [0, np.cos(phi/2), -np.sin(phi/2), 0], [0, np.sin(phi/2), np.cos(phi/2), 0], [0, 0, 0, np.exp(-1j*phi/2)]]) + + +def DoubleExcitation(phi): + r"""Double excitation rotation. + + Args: + phi (float): rotation angle + Returns: + array: the four-qubit Givens rotation describing the double excitation + """ + + c = math.cos(phi / 2) + s = math.sin(phi / 2) + + U = np.eye(16) + U[3, 3] = c # 3 (dec) = 0011 (bin) + U[3, 12] = -s # 12 (dec) = 1100 (bin) + U[12, 3] = s + U[12, 12] = c + + return U + + +def DoubleExcitationPlus(phi): + r"""Double excitation rotation with positive phase shift. + + Args: + phi (float): rotation angle + + Returns: + array: the four-qubit matrix describing the operation + """ + + c = math.cos(phi / 2) + s = math.sin(phi / 2) + e = cmath.exp(1j * phi / 2) + + U = e * np.eye(16, dtype=np.complex64) + U[3, 3] = c # 3 (dec) = 0011 (bin) + U[3, 12] = -s # 12 (dec) = 1100 (bin) + U[12, 3] = s + U[12, 12] = c + + return U + + +def DoubleExcitationMinus(phi): + r"""Double excitation rotation with negative phase shift. + + Args: + phi (float): rotation angle + Returns: + array: the four-qubit matrix describing the operation + """ + + c = math.cos(phi / 2) + s = math.sin(phi / 2) + e = cmath.exp(-1j * phi / 2) + + U = e * np.eye(16, dtype=np.complex64) + U[3, 3] = c # 3 (dec) = 0011 (bin) + U[3, 12] = -s # 12 (dec) = 1100 (bin) + U[12, 3] = s + U[12, 12] = c + + return U diff --git a/tests/ops/test_qubit_ops.py b/tests/ops/test_qubit_ops.py index 2098f87f9fb..009f74e0cbb 100644 --- a/tests/ops/test_qubit_ops.py +++ b/tests/ops/test_qubit_ops.py @@ -25,7 +25,8 @@ from pennylane.wires import Wires from gate_data import I, X, Y, Z, H, CNOT, SWAP, CZ, S, T, CSWAP, Toffoli, QFT, \ - ControlledPhaseShift, SingleExcitation, SingleExcitationPlus, SingleExcitationMinus + ControlledPhaseShift, SingleExcitation, SingleExcitationPlus, SingleExcitationMinus, \ + DoubleExcitation, DoubleExcitationPlus, DoubleExcitationMinus # Standard observables, their matrix representation, and eigenvalues @@ -1085,6 +1086,214 @@ def circuit(phi): ] +class TestDoubleExcitation: + @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4]) + def test_double_excitation_matrix(self, phi): + """Tests that the DoubleExcitation operation calculates the correct matrix""" + op = qml.DoubleExcitation(phi, wires=[0, 1, 2, 3]) + res = op.matrix + exp = DoubleExcitation(phi) + assert np.allclose(res, exp) + + @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4]) + def test_double_excitation_decomp(self, phi): + """Tests that the DoubleExcitation operation calculates the correct decomposition""" + op = qml.DoubleExcitation(phi, wires=[0, 1, 2, 3]) + decomp = op.decomposition(phi, wires=[0, 1, 2, 3]) + + mats = [m.matrix for m in decomp] + decomposed_matrix = mats[0] @ mats[1] + exp = DoubleExcitation(phi) + + assert np.allclose(decomposed_matrix, exp) + + @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4]) + def test_double_excitation_generator(self, phi): + """Tests that the DoubleExcitation operation calculates the correct generator""" + op = qml.DoubleExcitation(phi, wires=[0, 1, 2, 3]) + g, a = op.generator + + res = expm(1j * a * g * phi) + exp = DoubleExcitation(phi) + + assert np.allclose(res, exp) + + @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4]) + def test_double_excitation_plus_matrix(self, phi): + """Tests that the DoubleExcitationPlus operation calculates the correct matrix""" + op = qml.DoubleExcitationPlus(phi, wires=[0, 1, 2, 3]) + res = op.matrix + exp = DoubleExcitationPlus(phi) + assert np.allclose(res, exp) + + @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4]) + def test_double_excitation_plus_generator(self, phi): + """Tests that the DoubleExcitationPlus operation calculates the correct generator""" + op = qml.DoubleExcitationPlus(phi, wires=[0, 1, 2, 3]) + g, a = op.generator + + res = expm(1j * a * g * phi) + exp = DoubleExcitationPlus(phi) + + assert np.allclose(res, exp) + + @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4]) + def test_double_excitation_minus_matrix(self, phi): + """Tests that the DoubleExcitationMinus operation calculates the correct matrix""" + op = qml.DoubleExcitationMinus(phi, wires=[0, 1, 2, 3]) + res = op.matrix + exp = DoubleExcitationMinus(phi) + assert np.allclose(res, exp) + + @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4]) + def test_double_excitation_minus_generator(self, phi): + """Tests that the DoubleExcitationMinus operation calculates the correct generator""" + op = qml.DoubleExcitationMinus(phi, wires=[0, 1, 2, 3]) + g, a = op.generator + + res = expm(1j * a * g * phi) + exp = DoubleExcitationMinus(phi) + + assert np.allclose(res, exp) + + @pytest.mark.parametrize("excitation", [qml.DoubleExcitation, qml.DoubleExcitationPlus, + qml.DoubleExcitationMinus]) + def test_autograd(self, excitation): + """Tests that operations are computed correctly using the + autograd interface""" + + pytest.importorskip("autograd") + + dev = qml.device('default.qubit.autograd', wires=4) + state = np.array([0, 0, 0, -1 / np.sqrt(2), 0, 0, 0, 0, 0, 0, 0, 0, 1 / np.sqrt(2), 0, 0, + 0]) + + @qml.qnode(dev) + def circuit(phi): + qml.PauliX(wires=0) + qml.PauliX(wires=1) + excitation(phi, wires=[0, 1, 2, 3]) + + return qml.state() + + assert np.allclose(state, circuit(np.pi / 2)) + + + @pytest.mark.parametrize("excitation", [qml.DoubleExcitation, qml.DoubleExcitationPlus, + qml.DoubleExcitationMinus]) + def test_tf(self, excitation): + """Tests that operations are computed correctly using the + tensorflow interface""" + + pytest.importorskip("tensorflow") + + dev = qml.device('default.qubit.tf', wires=4) + state = np.array([0, 0, 0, -1 / np.sqrt(2), 0, 0, 0, 0, 0, 0, 0, 0, 1 / np.sqrt(2), 0, 0, + 0]) + + @qml.qnode(dev) + def circuit(phi): + qml.PauliX(wires=0) + qml.PauliX(wires=1) + excitation(phi, wires=[0, 1, 2, 3]) + + return qml.state() + + assert np.allclose(state, circuit(np.pi / 2)) + + @pytest.mark.parametrize("excitation", [qml.DoubleExcitation, qml.DoubleExcitationPlus, + qml.DoubleExcitationMinus]) + def test_jax(self, excitation): + """Tests that operations are computed correctly using the + jax interface""" + + pytest.importorskip("jax") + + dev = qml.device('default.qubit.jax', wires=4) + state = np.array([0, 0, 0, -1 / np.sqrt(2), 0, 0, 0, 0, 0, 0, 0, 0, 1 / np.sqrt(2), 0, 0, + 0]) + + @qml.qnode(dev) + def circuit(phi): + qml.PauliX(wires=0) + qml.PauliX(wires=1) + excitation(phi, wires=[0, 1, 2, 3]) + + return qml.state() + + assert np.allclose(state, circuit(np.pi / 2)) + + @pytest.mark.parametrize(("excitation", "phi"), [(qml.DoubleExcitation, -0.1), + (qml.DoubleExcitationPlus, 0.2), + (qml.DoubleExcitationMinus, np.pi / 4)]) + def test_autograd_grad(self, excitation, phi): + """Tests that gradients are computed correctly using the + autograd interface""" + + pytest.importorskip("autograd") + + dev = qml.device('default.qubit.autograd', wires=4) + + @qml.qnode(dev) + def circuit(phi): + qml.PauliX(wires=0) + qml.PauliX(wires=1) + excitation(phi, wires=[0, 1, 2, 3]) + + return qml.expval(qml.PauliZ(0)) + + assert np.allclose(qml.grad(circuit)(phi), np.sin(phi)) + + @pytest.mark.parametrize("diff_method", ["parameter-shift", "backprop"]) + @pytest.mark.parametrize(("excitation", "phi"), [(qml.DoubleExcitation, -0.1), + (qml.DoubleExcitationPlus, 0.2), + (qml.DoubleExcitationMinus, np.pi / 4)]) + def test_tf_grad(self, excitation, phi, diff_method): + """Tests that gradients are computed correctly using the + tensorflow interface""" + + tf = pytest.importorskip("tensorflow") + dev = qml.device('default.qubit.tf', wires=4) + + @qml.qnode(dev, interface="tf", diff_method=diff_method) + def circuit(phi): + qml.PauliX(wires=0) + qml.PauliX(wires=1) + excitation(phi, wires=[0, 1, 2, 3]) + return qml.expval(qml.PauliZ(0)) + + phi_t = tf.Variable(phi, dtype=tf.float64) + with tf.GradientTape() as tape: + res = circuit(phi_t) + + grad = tape.gradient(res, phi_t) + assert np.allclose(grad, np.sin(phi)) + + @pytest.mark.parametrize("diff_method", ["parameter-shift", "backprop"]) + @pytest.mark.parametrize(("excitation", "phi"), [(qml.DoubleExcitation, -0.1), + (qml.DoubleExcitationPlus, 0.2), + (qml.DoubleExcitationMinus, np.pi / 4)]) + def test_jax_grad(self, excitation, phi, diff_method): + """Tests that gradients and operations are computed correctly using the + jax interface""" + + if diff_method=="parameter-shift": + pytest.skip("JAX support for the parameter-shift method is still TBD") + + jax = pytest.importorskip("jax") + + dev = qml.device('default.qubit.jax', wires=4) + + @qml.qnode(dev, interface="jax", diff_method=diff_method) + def circuit(phi): + qml.PauliX(wires=0) + qml.PauliX(wires=1) + excitation(phi, wires=[0, 1, 2, 3]) + return qml.expval(qml.PauliZ(0)) + + assert np.allclose(jax.grad(circuit)(phi), np.sin(phi)) + + class TestPauliRot: """Test the PauliRot operation.""" diff --git a/tests/tape/test_reversible.py b/tests/tape/test_reversible.py index ba25b6d215f..a004e968a2a 100644 --- a/tests/tape/test_reversible.py +++ b/tests/tape/test_reversible.py @@ -1,512 +1,515 @@ -# Copyright 2018-2020 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. -"""Unit tests for the qubit parameter-shift QubitParamShiftTape""" -import pytest -from pennylane import numpy as np - -import pennylane as qml -from pennylane.interfaces.autograd import AutogradInterface -from pennylane.tape import JacobianTape, ReversibleTape -from pennylane import QNode, qnode -from pennylane.measure import MeasurementProcess - - -thetas = np.linspace(-2 * np.pi, 2 * np.pi, 8) - - -class TestReversibleTape: - """Unit tests for the reversible tape""" - - def test_diff_circuit_construction(self, mocker): - """Test that the diff circuit is correctly constructed""" - dev = qml.device("default.qubit", wires=2) - - with ReversibleTape() as tape: - qml.PauliX(wires=0) - qml.RX(0.542, wires=0) - qml.RY(0.542, wires=0) - qml.expval(qml.PauliZ(0)) - - spy = mocker.spy(dev, "execute") - tape.jacobian(dev) - - tape0 = spy.call_args_list[0][0][0] - tape1 = spy.call_args_list[1][0][0] - tape2 = spy.call_args_list[2][0][0] - - assert tape0 is tape - - assert len(tape1.operations) == 4 - assert len(tape1.measurements) == 1 - assert tape1.operations[0].name == "QubitStateVector" - assert tape1.operations[1].name == "RY.inv" - assert tape1.operations[2].name == "PauliX" - assert tape1.operations[3].name == "RY" - - assert len(tape2.operations) == 2 - assert len(tape1.measurements) == 1 - assert tape1.operations[0].name == "QubitStateVector" - assert tape2.operations[1].name == "PauliY" - - def test_rot_diff_circuit_construction(self, mocker): - """Test that the diff circuit is correctly constructed for the Rot gate""" - dev = qml.device("default.qubit", wires=2) - - with ReversibleTape() as tape: - qml.PauliX(wires=0) - qml.Rot(0.1, 0.2, 0.3, wires=0) - qml.expval(qml.PauliZ(0)) - - spy = mocker.spy(dev, "execute") - tape.jacobian(dev) - - tape0 = spy.call_args_list[0][0][0] - tape1 = spy.call_args_list[1][0][0] - tape2 = spy.call_args_list[2][0][0] - tape3 = spy.call_args_list[3][0][0] - - assert tape0 is tape - - assert len(tape1.operations) == 6 - assert len(tape1.measurements) == 1 - assert tape1.operations[0].name == "QubitStateVector" - assert tape1.operations[1].name == "RZ.inv" - assert tape1.operations[2].name == "RY.inv" - assert tape1.operations[3].name == "PauliZ" - assert tape1.operations[4].name == "RY" - assert tape1.operations[5].name == "RZ" - - assert len(tape2.operations) == 4 - assert len(tape1.measurements) == 1 - assert tape1.operations[0].name == "QubitStateVector" - assert tape2.operations[1].name == "RZ.inv" - assert tape2.operations[2].name == "PauliY" - assert tape2.operations[3].name == "RZ" - - assert len(tape3.operations) == 2 - assert len(tape1.measurements) == 1 - assert tape1.operations[0].name == "QubitStateVector" - assert tape3.operations[1].name == "PauliZ" - - @pytest.mark.parametrize("op, name", [(qml.CRX, "CRX"), (qml.CRY, "CRY"), (qml.CRZ, "CRZ")]) - def test_controlled_rotation_gates_exception(self, op, name): - """Tests that an exception is raised when a controlled - rotation gate is used with the ReversibleTape.""" - # TODO: remove this test when this support is added - dev = qml.device("default.qubit", wires=2) - - with ReversibleTape() as tape: - qml.PauliX(wires=0) - op(0.542, wires=[0, 1]) - qml.expval(qml.PauliZ(0)) - - with pytest.raises(ValueError, match="The {} gate is not currently supported".format(name)): - tape.jacobian(dev) - - def test_var_exception(self): - """Tests that an exception is raised when variance - is used with the ReversibleTape.""" - # TODO: remove this test when this support is added - dev = qml.device("default.qubit", wires=2) - - with ReversibleTape() as tape: - qml.PauliX(wires=0) - qml.RX(0.542, wires=0) - qml.var(qml.PauliZ(0)) - - with pytest.raises(ValueError, match="Variance is not supported"): - tape.jacobian(dev) - - def test_probs_exception(self): - """Tests that an exception is raised when probability - is used with the ReversibleTape.""" - # TODO: remove this test when this support is added - dev = qml.device("default.qubit", wires=2) - - with ReversibleTape() as tape: - qml.PauliX(wires=0) - qml.RX(0.542, wires=0) - qml.probs(wires=[0, 1]) - - with pytest.raises(ValueError, match="Probability is not supported"): - tape.jacobian(dev) - - def test_phaseshift_exception(self): - """Tests that an exception is raised when a PhaseShift gate - is used with the ReversibleTape.""" - # TODO: remove this test when this support is added - dev = qml.device("default.qubit", wires=1) - - with ReversibleTape() as tape: - qml.PauliX(wires=0) - qml.PhaseShift(0.542, wires=0) - qml.expval(qml.PauliZ(0)) - - with pytest.raises(ValueError, match="The PhaseShift gate is not currently supported"): - tape.jacobian(dev) - - -class TestGradients: - """Jacobian integration tests for qubit expectations.""" - - @pytest.mark.parametrize("theta", np.linspace(-2 * np.pi, 2 * np.pi, 7)) - @pytest.mark.parametrize("G", [qml.RX, qml.RY, qml.RZ]) - def test_pauli_rotation_gradient(self, G, theta, tol): - """Tests that the automatic gradients of Pauli rotations are correct.""" - dev = qml.device("default.qubit", wires=1) - - with ReversibleTape() as tape: - qml.QubitStateVector(np.array([1.0, -1.0]) / np.sqrt(2), wires=0) - G(theta, wires=[0]) - qml.expval(qml.PauliZ(0)) - - tape.trainable_params = {1} - - autograd_val = tape.jacobian(dev, method="analytic") - - # compare to finite differences - numeric_val = tape.jacobian(dev, method="numeric") - assert np.allclose(autograd_val, numeric_val, atol=tol, rtol=0) - - @pytest.mark.parametrize("theta", np.linspace(-2 * np.pi, 2 * np.pi, 7)) - def test_Rot_gradient(self, theta, tol): - """Tests that the automatic gradient of a arbitrary Euler-angle-parameterized gate is correct.""" - dev = qml.device("default.qubit", wires=1) - params = np.array([theta, theta ** 3, np.sqrt(2) * theta]) - - with ReversibleTape() as tape: - qml.QubitStateVector(np.array([1.0, -1.0]) / np.sqrt(2), wires=0) - qml.Rot(*params, wires=[0]) - qml.expval(qml.PauliZ(0)) - - tape.trainable_params = {1, 2, 3} - - autograd_val = tape.jacobian(dev, method="analytic") - - # compare to finite differences - numeric_val = tape.jacobian(dev, method="numeric") - assert np.allclose(autograd_val, numeric_val, atol=tol, rtol=0) - - @pytest.mark.parametrize("par", [1, -2, 1.623, -0.051, 0]) # intergers, floats, zero - def test_ry_gradient(self, par, mocker, tol): - """Test that the gradient of the RY gate matches the exact analytic - formula. Further, make sure the correct gradient methods - are being called.""" - - with ReversibleTape() as tape: - qml.RY(par, wires=[0]) - qml.expval(qml.PauliX(0)) - - tape.trainable_params = {0} - - dev = qml.device("default.qubit", wires=1) - - spy_numeric = mocker.spy(tape, "numeric_pd") - spy_analytic = mocker.spy(tape, "analytic_pd") - - # gradients - exact = np.cos(par) - grad_F = tape.jacobian(dev, method="numeric") - - spy_numeric.assert_called() - spy_analytic.assert_not_called() - - spy_device = mocker.spy(tape, "execute_device") - grad_A = tape.jacobian(dev, method="analytic") - - spy_analytic.assert_called() - spy_device.assert_called_once() # check that the state was only pre-computed once - - # different methods must agree - assert np.allclose(grad_F, exact, atol=tol, rtol=0) - assert np.allclose(grad_A, exact, atol=tol, rtol=0) - - def test_rx_gradient(self, tol): - """Test that the gradient of the RX gate matches the known formula.""" - dev = qml.device("default.qubit", wires=2) - a = 0.7418 - - with ReversibleTape() as tape: - qml.RX(a, wires=0) - qml.expval(qml.PauliZ(0)) - - circuit_output = tape.execute(dev) - expected_output = np.cos(a) - assert np.allclose(circuit_output, expected_output, atol=tol, rtol=0) - - # circuit jacobians - circuit_jacobian = tape.jacobian(dev, method="analytic") - expected_jacobian = -np.sin(a) - assert np.allclose(circuit_jacobian, expected_jacobian, atol=tol, rtol=0) - - def test_multiple_rx_gradient(self, tol): - """Tests that the gradient of multiple RX gates in a circuit - yeilds the correct result.""" - dev = qml.device("default.qubit", wires=3) - params = np.array([np.pi, np.pi / 2, np.pi / 3]) - - with ReversibleTape() as tape: - qml.RX(params[0], wires=0) - qml.RX(params[1], wires=1) - qml.RX(params[2], wires=2) - - for idx in range(3): - qml.expval(qml.PauliZ(idx)) - - circuit_output = tape.execute(dev) - expected_output = np.cos(params) - assert np.allclose(circuit_output, expected_output, atol=tol, rtol=0) - - # circuit jacobians - circuit_jacobian = tape.jacobian(dev, method="analytic") - expected_jacobian = -np.diag(np.sin(params)) - assert np.allclose(circuit_jacobian, expected_jacobian, atol=tol, rtol=0) - - qubit_ops = [getattr(qml, name) for name in qml.ops._qubit__ops__] - analytic_qubit_ops = {cls for cls in qubit_ops if cls.grad_method == "A"} - analytic_qubit_ops = analytic_qubit_ops - { - qml.CRX, - qml.CRY, - qml.CRZ, - qml.CRot, - qml.PhaseShift, - qml.ControlledPhaseShift, - qml.PauliRot, - qml.MultiRZ, - qml.U1, - qml.U2, - qml.U3, - qml.SingleExcitation, - qml.SingleExcitationPlus, - qml.SingleExcitationMinus, - } - - @pytest.mark.parametrize("obs", [qml.PauliX, qml.PauliY]) - @pytest.mark.parametrize("op", analytic_qubit_ops) - def test_gradients(self, op, obs, mocker, tol): - """Tests that the gradients of circuits match between the - finite difference and analytic methods.""" - args = np.linspace(0.2, 0.5, op.num_params) - - with ReversibleTape() as tape: - qml.Hadamard(wires=0) - qml.RX(0.543, wires=0) - qml.CNOT(wires=[0, 1]) - - op(*args, wires=range(op.num_wires)) - - qml.Rot(1.3, -2.3, 0.5, wires=[0]) - qml.RZ(-0.5, wires=0) - qml.RY(0.5, wires=1) - qml.CNOT(wires=[0, 1]) - - qml.expval(obs(wires=0)) - qml.expval(qml.PauliZ(wires=1)) - - dev = qml.device("default.qubit", wires=2) - res = tape.execute(dev) - - tape._update_gradient_info() - tape.trainable_params = set(range(1, 1 + op.num_params)) - - # check that every parameter is analytic - for i in range(op.num_params): - assert tape._par_info[1 + i]["grad_method"][0] == "A" - - grad_F = tape.jacobian(dev, method="numeric") - - spy = mocker.spy(ReversibleTape, "analytic_pd") - spy_execute = mocker.spy(tape, "execute_device") - grad_A = tape.jacobian(dev, method="analytic") - spy.assert_called() - - # check that the execute device method has only been called - # once, for all parameters. - spy_execute.assert_called_once() - - assert np.allclose(grad_A, grad_F, atol=tol, rtol=0) - - def test_gradient_gate_with_multiple_parameters(self, tol): - """Tests that gates with multiple free parameters yield correct gradients.""" - x, y, z = [0.5, 0.3, -0.7] - - with ReversibleTape() as tape: - qml.RX(0.4, wires=[0]) - qml.Rot(x, y, z, wires=[0]) - qml.RY(-0.2, wires=[0]) - qml.expval(qml.PauliZ(0)) - - tape.trainable_params = {1, 2, 3} - - dev = qml.device("default.qubit", wires=1) - grad_A = tape.jacobian(dev, method="analytic") - grad_F = tape.jacobian(dev, method="numeric") - - # gradient has the correct shape and every element is nonzero - assert grad_A.shape == (1, 3) - assert np.count_nonzero(grad_A) == 3 - # the different methods agree - assert np.allclose(grad_A, grad_F, atol=tol, rtol=0) - - -class TestQNodeIntegration: - """Test QNode integration with the reversible method""" - - def test_qnode(self, mocker, tol): - """Test that specifying diff_method allows the reversible - method to be selected""" - args = np.array([0.54, 0.1, 0.5], requires_grad=True) - dev = qml.device("default.qubit", wires=2) - - def circuit(x, y, z): - qml.Hadamard(wires=0) - qml.RX(0.543, wires=0) - qml.CNOT(wires=[0, 1]) - - qml.Rot(x, y, z, wires=0) - - qml.Rot(1.3, -2.3, 0.5, wires=[0]) - qml.RZ(-0.5, wires=0) - qml.RY(0.5, wires=1) - qml.CNOT(wires=[0, 1]) - - return qml.expval(qml.PauliX(0) @ qml.PauliZ(1)) - - qnode1 = QNode(circuit, dev, diff_method="reversible") - spy = mocker.spy(ReversibleTape, "analytic_pd") - - grad_fn = qml.grad(qnode1) - grad_A = grad_fn(*args) - - spy.assert_called() - assert isinstance(qnode1.qtape, ReversibleTape) - - qnode2 = QNode(circuit, dev, diff_method="finite-diff") - grad_fn = qml.grad(qnode2) - grad_F = grad_fn(*args) - - assert not isinstance(qnode2.qtape, ReversibleTape) - assert np.allclose(grad_A, grad_F, atol=tol, rtol=0) - - @pytest.mark.parametrize("reused_p", thetas ** 3 / 19) - @pytest.mark.parametrize("other_p", thetas ** 2 / 1) - def test_fanout_multiple_params(self, reused_p, other_p, tol): - """Tests that the correct gradient is computed for qnodes which - use the same parameter in multiple gates.""" - - from gate_data import Rotx as Rx, Roty as Ry, Rotz as Rz - - def expZ(state): - return np.abs(state[0]) ** 2 - np.abs(state[1]) ** 2 - - dev = qml.device("default.qubit", wires=1) - extra_param = np.array(0.31, requires_grad=False) - - @qnode(dev) - def cost(p1, p2): - qml.RX(extra_param, wires=[0]) - qml.RY(p1, wires=[0]) - qml.RZ(p2, wires=[0]) - qml.RX(p1, wires=[0]) - return qml.expval(qml.PauliZ(0)) - - zero_state = np.array([1.0, 0.0]) - - # analytic gradient - grad_fn = qml.grad(cost) - grad_A = grad_fn(reused_p, other_p) - - # manual gradient - grad_true0 = ( - expZ( - Rx(reused_p) @ Rz(other_p) @ Ry(reused_p + np.pi / 2) @ Rx(extra_param) @ zero_state - ) - - expZ( - Rx(reused_p) @ Rz(other_p) @ Ry(reused_p - np.pi / 2) @ Rx(extra_param) @ zero_state - ) - ) / 2 - grad_true1 = ( - expZ( - Rx(reused_p + np.pi / 2) @ Rz(other_p) @ Ry(reused_p) @ Rx(extra_param) @ zero_state - ) - - expZ( - Rx(reused_p - np.pi / 2) @ Rz(other_p) @ Ry(reused_p) @ Rx(extra_param) @ zero_state - ) - ) / 2 - expected = grad_true0 + grad_true1 # product rule - - assert np.allclose(grad_A[0], expected, atol=tol, rtol=0) - - def test_gradient_repeated_gate_parameters(self, mocker, tol): - """Tests that repeated use of a free parameter in a - multi-parameter gate yield correct gradients.""" - dev = qml.device("default.qubit", wires=1) - params = np.array([0.8, 1.3], requires_grad=True) - - def circuit(params): - qml.RX(np.array(np.pi / 4, requires_grad=False), wires=[0]) - qml.Rot(params[1], params[0], 2 * params[0], wires=[0]) - return qml.expval(qml.PauliX(0)) - - spy_numeric = mocker.spy(JacobianTape, "numeric_pd") - spy_analytic = mocker.spy(ReversibleTape, "analytic_pd") - - cost = QNode(circuit, dev, diff_method="finite-diff") - grad_fn = qml.grad(cost) - grad_F = grad_fn(params) - - spy_numeric.assert_called() - spy_analytic.assert_not_called() - - cost = QNode(circuit, dev, diff_method="reversible") - grad_fn = qml.grad(cost) - grad_A = grad_fn(params) - - spy_analytic.assert_called() - - # the different methods agree - assert np.allclose(grad_A, grad_F, atol=tol, rtol=0) - - -class TestHelperFunctions: - """Tests for additional helper functions.""" - - one_qubit_vec1 = np.array([1, 1]) - one_qubit_vec2 = np.array([1, 1j]) - two_qubit_vec = np.array([1, 1, 1, -1]) - single_qubit_obs1 = qml.PauliZ(0) - single_qubit_obs2 = qml.PauliY(0) - two_qubit_obs = qml.Hermitian(np.eye(4), wires=[0, 1]) - - @pytest.mark.parametrize( - "wires, vec1, obs, vec2, expected", - [ - (1, one_qubit_vec1, single_qubit_obs1, one_qubit_vec1, 0), - (1, one_qubit_vec2, single_qubit_obs1, one_qubit_vec2, 0), - (1, one_qubit_vec1, single_qubit_obs1, one_qubit_vec2, 1 - 1j), - (1, one_qubit_vec2, single_qubit_obs1, one_qubit_vec1, 1 + 1j), - (1, one_qubit_vec1, single_qubit_obs2, one_qubit_vec1, 0), - (1, one_qubit_vec2, single_qubit_obs2, one_qubit_vec2, 2), - (1, one_qubit_vec1, single_qubit_obs2, one_qubit_vec2, 1 + 1j), - (1, one_qubit_vec2, single_qubit_obs2, one_qubit_vec1, 1 - 1j), - (2, two_qubit_vec, single_qubit_obs1, two_qubit_vec, 0), - (2, two_qubit_vec, single_qubit_obs2, two_qubit_vec, 0), - (2, two_qubit_vec, two_qubit_obs, two_qubit_vec, 4), - ], - ) - def test_matrix_elem(self, wires, vec1, obs, vec2, expected): - """Tests for the helper function _matrix_elem""" - tape = ReversibleTape() - res = tape._matrix_elem(vec1, obs, vec2, qml.wires.Wires(range(wires))) - assert res == expected +# Copyright 2018-2020 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. +"""Unit tests for the qubit parameter-shift QubitParamShiftTape""" +import pytest +from pennylane import numpy as np + +import pennylane as qml +from pennylane.interfaces.autograd import AutogradInterface +from pennylane.tape import JacobianTape, ReversibleTape +from pennylane import QNode, qnode +from pennylane.measure import MeasurementProcess + + +thetas = np.linspace(-2 * np.pi, 2 * np.pi, 8) + + +class TestReversibleTape: + """Unit tests for the reversible tape""" + + def test_diff_circuit_construction(self, mocker): + """Test that the diff circuit is correctly constructed""" + dev = qml.device("default.qubit", wires=2) + + with ReversibleTape() as tape: + qml.PauliX(wires=0) + qml.RX(0.542, wires=0) + qml.RY(0.542, wires=0) + qml.expval(qml.PauliZ(0)) + + spy = mocker.spy(dev, "execute") + tape.jacobian(dev) + + tape0 = spy.call_args_list[0][0][0] + tape1 = spy.call_args_list[1][0][0] + tape2 = spy.call_args_list[2][0][0] + + assert tape0 is tape + + assert len(tape1.operations) == 4 + assert len(tape1.measurements) == 1 + assert tape1.operations[0].name == "QubitStateVector" + assert tape1.operations[1].name == "RY.inv" + assert tape1.operations[2].name == "PauliX" + assert tape1.operations[3].name == "RY" + + assert len(tape2.operations) == 2 + assert len(tape1.measurements) == 1 + assert tape1.operations[0].name == "QubitStateVector" + assert tape2.operations[1].name == "PauliY" + + def test_rot_diff_circuit_construction(self, mocker): + """Test that the diff circuit is correctly constructed for the Rot gate""" + dev = qml.device("default.qubit", wires=2) + + with ReversibleTape() as tape: + qml.PauliX(wires=0) + qml.Rot(0.1, 0.2, 0.3, wires=0) + qml.expval(qml.PauliZ(0)) + + spy = mocker.spy(dev, "execute") + tape.jacobian(dev) + + tape0 = spy.call_args_list[0][0][0] + tape1 = spy.call_args_list[1][0][0] + tape2 = spy.call_args_list[2][0][0] + tape3 = spy.call_args_list[3][0][0] + + assert tape0 is tape + + assert len(tape1.operations) == 6 + assert len(tape1.measurements) == 1 + assert tape1.operations[0].name == "QubitStateVector" + assert tape1.operations[1].name == "RZ.inv" + assert tape1.operations[2].name == "RY.inv" + assert tape1.operations[3].name == "PauliZ" + assert tape1.operations[4].name == "RY" + assert tape1.operations[5].name == "RZ" + + assert len(tape2.operations) == 4 + assert len(tape1.measurements) == 1 + assert tape1.operations[0].name == "QubitStateVector" + assert tape2.operations[1].name == "RZ.inv" + assert tape2.operations[2].name == "PauliY" + assert tape2.operations[3].name == "RZ" + + assert len(tape3.operations) == 2 + assert len(tape1.measurements) == 1 + assert tape1.operations[0].name == "QubitStateVector" + assert tape3.operations[1].name == "PauliZ" + + @pytest.mark.parametrize("op, name", [(qml.CRX, "CRX"), (qml.CRY, "CRY"), (qml.CRZ, "CRZ")]) + def test_controlled_rotation_gates_exception(self, op, name): + """Tests that an exception is raised when a controlled + rotation gate is used with the ReversibleTape.""" + # TODO: remove this test when this support is added + dev = qml.device("default.qubit", wires=2) + + with ReversibleTape() as tape: + qml.PauliX(wires=0) + op(0.542, wires=[0, 1]) + qml.expval(qml.PauliZ(0)) + + with pytest.raises(ValueError, match="The {} gate is not currently supported".format(name)): + tape.jacobian(dev) + + def test_var_exception(self): + """Tests that an exception is raised when variance + is used with the ReversibleTape.""" + # TODO: remove this test when this support is added + dev = qml.device("default.qubit", wires=2) + + with ReversibleTape() as tape: + qml.PauliX(wires=0) + qml.RX(0.542, wires=0) + qml.var(qml.PauliZ(0)) + + with pytest.raises(ValueError, match="Variance is not supported"): + tape.jacobian(dev) + + def test_probs_exception(self): + """Tests that an exception is raised when probability + is used with the ReversibleTape.""" + # TODO: remove this test when this support is added + dev = qml.device("default.qubit", wires=2) + + with ReversibleTape() as tape: + qml.PauliX(wires=0) + qml.RX(0.542, wires=0) + qml.probs(wires=[0, 1]) + + with pytest.raises(ValueError, match="Probability is not supported"): + tape.jacobian(dev) + + def test_phaseshift_exception(self): + """Tests that an exception is raised when a PhaseShift gate + is used with the ReversibleTape.""" + # TODO: remove this test when this support is added + dev = qml.device("default.qubit", wires=1) + + with ReversibleTape() as tape: + qml.PauliX(wires=0) + qml.PhaseShift(0.542, wires=0) + qml.expval(qml.PauliZ(0)) + + with pytest.raises(ValueError, match="The PhaseShift gate is not currently supported"): + tape.jacobian(dev) + + +class TestGradients: + """Jacobian integration tests for qubit expectations.""" + + @pytest.mark.parametrize("theta", np.linspace(-2 * np.pi, 2 * np.pi, 7)) + @pytest.mark.parametrize("G", [qml.RX, qml.RY, qml.RZ]) + def test_pauli_rotation_gradient(self, G, theta, tol): + """Tests that the automatic gradients of Pauli rotations are correct.""" + dev = qml.device("default.qubit", wires=1) + + with ReversibleTape() as tape: + qml.QubitStateVector(np.array([1.0, -1.0]) / np.sqrt(2), wires=0) + G(theta, wires=[0]) + qml.expval(qml.PauliZ(0)) + + tape.trainable_params = {1} + + autograd_val = tape.jacobian(dev, method="analytic") + + # compare to finite differences + numeric_val = tape.jacobian(dev, method="numeric") + assert np.allclose(autograd_val, numeric_val, atol=tol, rtol=0) + + @pytest.mark.parametrize("theta", np.linspace(-2 * np.pi, 2 * np.pi, 7)) + def test_Rot_gradient(self, theta, tol): + """Tests that the automatic gradient of a arbitrary Euler-angle-parameterized gate is correct.""" + dev = qml.device("default.qubit", wires=1) + params = np.array([theta, theta ** 3, np.sqrt(2) * theta]) + + with ReversibleTape() as tape: + qml.QubitStateVector(np.array([1.0, -1.0]) / np.sqrt(2), wires=0) + qml.Rot(*params, wires=[0]) + qml.expval(qml.PauliZ(0)) + + tape.trainable_params = {1, 2, 3} + + autograd_val = tape.jacobian(dev, method="analytic") + + # compare to finite differences + numeric_val = tape.jacobian(dev, method="numeric") + assert np.allclose(autograd_val, numeric_val, atol=tol, rtol=0) + + @pytest.mark.parametrize("par", [1, -2, 1.623, -0.051, 0]) # intergers, floats, zero + def test_ry_gradient(self, par, mocker, tol): + """Test that the gradient of the RY gate matches the exact analytic + formula. Further, make sure the correct gradient methods + are being called.""" + + with ReversibleTape() as tape: + qml.RY(par, wires=[0]) + qml.expval(qml.PauliX(0)) + + tape.trainable_params = {0} + + dev = qml.device("default.qubit", wires=1) + + spy_numeric = mocker.spy(tape, "numeric_pd") + spy_analytic = mocker.spy(tape, "analytic_pd") + + # gradients + exact = np.cos(par) + grad_F = tape.jacobian(dev, method="numeric") + + spy_numeric.assert_called() + spy_analytic.assert_not_called() + + spy_device = mocker.spy(tape, "execute_device") + grad_A = tape.jacobian(dev, method="analytic") + + spy_analytic.assert_called() + spy_device.assert_called_once() # check that the state was only pre-computed once + + # different methods must agree + assert np.allclose(grad_F, exact, atol=tol, rtol=0) + assert np.allclose(grad_A, exact, atol=tol, rtol=0) + + def test_rx_gradient(self, tol): + """Test that the gradient of the RX gate matches the known formula.""" + dev = qml.device("default.qubit", wires=2) + a = 0.7418 + + with ReversibleTape() as tape: + qml.RX(a, wires=0) + qml.expval(qml.PauliZ(0)) + + circuit_output = tape.execute(dev) + expected_output = np.cos(a) + assert np.allclose(circuit_output, expected_output, atol=tol, rtol=0) + + # circuit jacobians + circuit_jacobian = tape.jacobian(dev, method="analytic") + expected_jacobian = -np.sin(a) + assert np.allclose(circuit_jacobian, expected_jacobian, atol=tol, rtol=0) + + def test_multiple_rx_gradient(self, tol): + """Tests that the gradient of multiple RX gates in a circuit + yeilds the correct result.""" + dev = qml.device("default.qubit", wires=3) + params = np.array([np.pi, np.pi / 2, np.pi / 3]) + + with ReversibleTape() as tape: + qml.RX(params[0], wires=0) + qml.RX(params[1], wires=1) + qml.RX(params[2], wires=2) + + for idx in range(3): + qml.expval(qml.PauliZ(idx)) + + circuit_output = tape.execute(dev) + expected_output = np.cos(params) + assert np.allclose(circuit_output, expected_output, atol=tol, rtol=0) + + # circuit jacobians + circuit_jacobian = tape.jacobian(dev, method="analytic") + expected_jacobian = -np.diag(np.sin(params)) + assert np.allclose(circuit_jacobian, expected_jacobian, atol=tol, rtol=0) + + qubit_ops = [getattr(qml, name) for name in qml.ops._qubit__ops__] + analytic_qubit_ops = {cls for cls in qubit_ops if cls.grad_method == "A"} + analytic_qubit_ops = analytic_qubit_ops - { + qml.CRX, + qml.CRY, + qml.CRZ, + qml.CRot, + qml.PhaseShift, + qml.ControlledPhaseShift, + qml.PauliRot, + qml.MultiRZ, + qml.U1, + qml.U2, + qml.U3, + qml.SingleExcitation, + qml.SingleExcitationPlus, + qml.SingleExcitationMinus, + qml.DoubleExcitation, + qml.DoubleExcitationPlus, + qml.DoubleExcitationMinus, + } + + @pytest.mark.parametrize("obs", [qml.PauliX, qml.PauliY]) + @pytest.mark.parametrize("op", analytic_qubit_ops) + def test_gradients(self, op, obs, mocker, tol): + """Tests that the gradients of circuits match between the + finite difference and analytic methods.""" + args = np.linspace(0.2, 0.5, op.num_params) + + with ReversibleTape() as tape: + qml.Hadamard(wires=0) + qml.RX(0.543, wires=0) + qml.CNOT(wires=[0, 1]) + + op(*args, wires=range(op.num_wires)) + + qml.Rot(1.3, -2.3, 0.5, wires=[0]) + qml.RZ(-0.5, wires=0) + qml.RY(0.5, wires=1) + qml.CNOT(wires=[0, 1]) + + qml.expval(obs(wires=0)) + qml.expval(qml.PauliZ(wires=1)) + + dev = qml.device("default.qubit", wires=2) + res = tape.execute(dev) + + tape._update_gradient_info() + tape.trainable_params = set(range(1, 1 + op.num_params)) + + # check that every parameter is analytic + for i in range(op.num_params): + assert tape._par_info[1 + i]["grad_method"][0] == "A" + + grad_F = tape.jacobian(dev, method="numeric") + + spy = mocker.spy(ReversibleTape, "analytic_pd") + spy_execute = mocker.spy(tape, "execute_device") + grad_A = tape.jacobian(dev, method="analytic") + spy.assert_called() + + # check that the execute device method has only been called + # once, for all parameters. + spy_execute.assert_called_once() + + assert np.allclose(grad_A, grad_F, atol=tol, rtol=0) + + def test_gradient_gate_with_multiple_parameters(self, tol): + """Tests that gates with multiple free parameters yield correct gradients.""" + x, y, z = [0.5, 0.3, -0.7] + + with ReversibleTape() as tape: + qml.RX(0.4, wires=[0]) + qml.Rot(x, y, z, wires=[0]) + qml.RY(-0.2, wires=[0]) + qml.expval(qml.PauliZ(0)) + + tape.trainable_params = {1, 2, 3} + + dev = qml.device("default.qubit", wires=1) + grad_A = tape.jacobian(dev, method="analytic") + grad_F = tape.jacobian(dev, method="numeric") + + # gradient has the correct shape and every element is nonzero + assert grad_A.shape == (1, 3) + assert np.count_nonzero(grad_A) == 3 + # the different methods agree + assert np.allclose(grad_A, grad_F, atol=tol, rtol=0) + + +class TestQNodeIntegration: + """Test QNode integration with the reversible method""" + + def test_qnode(self, mocker, tol): + """Test that specifying diff_method allows the reversible + method to be selected""" + args = np.array([0.54, 0.1, 0.5], requires_grad=True) + dev = qml.device("default.qubit", wires=2) + + def circuit(x, y, z): + qml.Hadamard(wires=0) + qml.RX(0.543, wires=0) + qml.CNOT(wires=[0, 1]) + + qml.Rot(x, y, z, wires=0) + + qml.Rot(1.3, -2.3, 0.5, wires=[0]) + qml.RZ(-0.5, wires=0) + qml.RY(0.5, wires=1) + qml.CNOT(wires=[0, 1]) + + return qml.expval(qml.PauliX(0) @ qml.PauliZ(1)) + + qnode1 = QNode(circuit, dev, diff_method="reversible") + spy = mocker.spy(ReversibleTape, "analytic_pd") + + grad_fn = qml.grad(qnode1) + grad_A = grad_fn(*args) + + spy.assert_called() + assert isinstance(qnode1.qtape, ReversibleTape) + + qnode2 = QNode(circuit, dev, diff_method="finite-diff") + grad_fn = qml.grad(qnode2) + grad_F = grad_fn(*args) + + assert not isinstance(qnode2.qtape, ReversibleTape) + assert np.allclose(grad_A, grad_F, atol=tol, rtol=0) + + @pytest.mark.parametrize("reused_p", thetas ** 3 / 19) + @pytest.mark.parametrize("other_p", thetas ** 2 / 1) + def test_fanout_multiple_params(self, reused_p, other_p, tol): + """Tests that the correct gradient is computed for qnodes which + use the same parameter in multiple gates.""" + + from gate_data import Rotx as Rx, Roty as Ry, Rotz as Rz + + def expZ(state): + return np.abs(state[0]) ** 2 - np.abs(state[1]) ** 2 + + dev = qml.device("default.qubit", wires=1) + extra_param = np.array(0.31, requires_grad=False) + + @qnode(dev) + def cost(p1, p2): + qml.RX(extra_param, wires=[0]) + qml.RY(p1, wires=[0]) + qml.RZ(p2, wires=[0]) + qml.RX(p1, wires=[0]) + return qml.expval(qml.PauliZ(0)) + + zero_state = np.array([1.0, 0.0]) + + # analytic gradient + grad_fn = qml.grad(cost) + grad_A = grad_fn(reused_p, other_p) + + # manual gradient + grad_true0 = ( + expZ( + Rx(reused_p) @ Rz(other_p) @ Ry(reused_p + np.pi / 2) @ Rx(extra_param) @ zero_state + ) + - expZ( + Rx(reused_p) @ Rz(other_p) @ Ry(reused_p - np.pi / 2) @ Rx(extra_param) @ zero_state + ) + ) / 2 + grad_true1 = ( + expZ( + Rx(reused_p + np.pi / 2) @ Rz(other_p) @ Ry(reused_p) @ Rx(extra_param) @ zero_state + ) + - expZ( + Rx(reused_p - np.pi / 2) @ Rz(other_p) @ Ry(reused_p) @ Rx(extra_param) @ zero_state + ) + ) / 2 + expected = grad_true0 + grad_true1 # product rule + + assert np.allclose(grad_A[0], expected, atol=tol, rtol=0) + + def test_gradient_repeated_gate_parameters(self, mocker, tol): + """Tests that repeated use of a free parameter in a + multi-parameter gate yield correct gradients.""" + dev = qml.device("default.qubit", wires=1) + params = np.array([0.8, 1.3], requires_grad=True) + + def circuit(params): + qml.RX(np.array(np.pi / 4, requires_grad=False), wires=[0]) + qml.Rot(params[1], params[0], 2 * params[0], wires=[0]) + return qml.expval(qml.PauliX(0)) + + spy_numeric = mocker.spy(JacobianTape, "numeric_pd") + spy_analytic = mocker.spy(ReversibleTape, "analytic_pd") + + cost = QNode(circuit, dev, diff_method="finite-diff") + grad_fn = qml.grad(cost) + grad_F = grad_fn(params) + + spy_numeric.assert_called() + spy_analytic.assert_not_called() + + cost = QNode(circuit, dev, diff_method="reversible") + grad_fn = qml.grad(cost) + grad_A = grad_fn(params) + + spy_analytic.assert_called() + + # the different methods agree + assert np.allclose(grad_A, grad_F, atol=tol, rtol=0) + + +class TestHelperFunctions: + """Tests for additional helper functions.""" + + one_qubit_vec1 = np.array([1, 1]) + one_qubit_vec2 = np.array([1, 1j]) + two_qubit_vec = np.array([1, 1, 1, -1]) + single_qubit_obs1 = qml.PauliZ(0) + single_qubit_obs2 = qml.PauliY(0) + two_qubit_obs = qml.Hermitian(np.eye(4), wires=[0, 1]) + + @pytest.mark.parametrize( + "wires, vec1, obs, vec2, expected", + [ + (1, one_qubit_vec1, single_qubit_obs1, one_qubit_vec1, 0), + (1, one_qubit_vec2, single_qubit_obs1, one_qubit_vec2, 0), + (1, one_qubit_vec1, single_qubit_obs1, one_qubit_vec2, 1 - 1j), + (1, one_qubit_vec2, single_qubit_obs1, one_qubit_vec1, 1 + 1j), + (1, one_qubit_vec1, single_qubit_obs2, one_qubit_vec1, 0), + (1, one_qubit_vec2, single_qubit_obs2, one_qubit_vec2, 2), + (1, one_qubit_vec1, single_qubit_obs2, one_qubit_vec2, 1 + 1j), + (1, one_qubit_vec2, single_qubit_obs2, one_qubit_vec1, 1 - 1j), + (2, two_qubit_vec, single_qubit_obs1, two_qubit_vec, 0), + (2, two_qubit_vec, single_qubit_obs2, two_qubit_vec, 0), + (2, two_qubit_vec, two_qubit_obs, two_qubit_vec, 4), + ], + ) + def test_matrix_elem(self, wires, vec1, obs, vec2, expected): + """Tests for the helper function _matrix_elem""" + tape = ReversibleTape() + res = tape._matrix_elem(vec1, obs, vec2, qml.wires.Wires(range(wires))) + assert res == expected