Adds elementary gate decomposition for single-qubit QubitUnitary

.github/CHANGELOG.md

@@ -2,6 +2,61 @@
 
 <h3>New features since last release</h3>
 
+* A decomposition has been added to ``QubitUnitary`` that makes the
+  single-qubit case fully differentiable in all interfaces. Furthermore,
+  a quantum function transform, ``unitary_to_rot()``, has been added to decompose all
+  single-qubit instances of ``QubitUnitary`` in a quantum circuit.
+  [(#1427)](https://github.com/PennyLaneAI/pennylane/pull/1427)
+
+  Instances of ``QubitUnitary`` may now be decomposed directly to ``Rot``
+  operations, or ``RZ`` operations if the input matrix is diagonal. For
+  example, let
+
+  ```python
+  >>> U = np.array([
+      [-0.28829348-0.78829734j,  0.30364367+0.45085995j],
+      [ 0.53396245-0.10177564j,  0.76279558-0.35024096j]
+  ])
+  ```
+
+  Then, we can compute the decomposition as:
+
+  ```pycon
+  >>> qml.QubitUnitary.decomposition(U, wires=0)
+  [Rot(-0.24209530281458358, 1.1493817777199102, 1.733058145303424, wires=[0])]
+  ```
+
+  We can also apply the transform directly to a quantum function, and compute the
+  gradients of parameters used to construct the unitary matrices.
+
+  ```python
+  def qfunc_with_qubit_unitary(angles):
+      z, x = angles[0], angles[1]
+
+      Z_mat = np.array([[np.exp(-1j * z / 2), 0.0], [0.0, np.exp(1j * z / 2)]])
+
+      c = np.cos(x / 2)
+      s = np.sin(x / 2) * 1j
+      X_mat = np.array([[c, -s], [-s, c]])
+
+      qml.Hadamard(wires="a")
+      qml.QubitUnitary(Z_mat, wires="a")
+      qml.QubitUnitary(X_mat, wires="b")
+      qml.CNOT(wires=["b", "a"])
+      return qml.expval(qml.PauliX(wires="a"))
+  ```
+
+  ```pycon
+  >>> dev = qml.device("default.qubit", wires=["a", "b"])
+  >>> transformed_qfunc = qml.transforms.unitary_to_rot(qfunc_with_qubit_unitary)
+  >>> transformed_qnode = qml.QNode(transformed_qfunc, dev)
+  >>> input = np.array([0.3, 0.4], requires_grad=True)
+  >>> transformed_qnode(input)
+  tensor(0.95533649, requires_grad=True)
+  >>> qml.grad(transformed_qnode)(input)
+  array([-0.29552021,  0.        ])
+  ```
+
 * The new ``qml.apply`` function can be used to add operations that might have
   already been instantiated elsewhere to the QNode and other queuing contexts:
   [(#1433)](https://github.com/PennyLaneAI/pennylane/pull/1433)

pennylane/ops/qubit.py

@@ -2159,6 +2159,17 @@ def _matrix(cls, *params):
 
         return U
 
+    @staticmethod
+    def decomposition(U, wires):
+        # Decomposes arbitrary single-qubit unitaries as Rot gates (RZ - RY - RZ format),
+        # or a single RZ for diagonal matrices.
+        if qml.math.shape(U) == (2, 2):
+            wire = Wires(wires)[0]
+            decomp_ops = qml.transforms.decompositions.zyz_decomposition(U, wire)
+            return decomp_ops
+
+        raise NotImplementedError("Decompositions only supported for single-qubit unitaries")
+
     def adjoint(self):
         return QubitUnitary(qml.math.T(qml.math.conj(self.matrix)), wires=self.wires)
 

pennylane/transforms/__init__.py

@@ -45,10 +45,19 @@
 
     ~adjoint
     ~ctrl
+    ~transforms.unitary_to_rot
     ~transforms.invisible
     ~apply_controlled_Q
     ~quantum_monte_carlo
 
+There are also utility functions and decompositions available that assist with
+both transforms, and decompositions within the larger PennyLane codebase.
+
+.. autosummary::
+    :toctree: api
+
+    ~transforms.zyz_decomposition
+
 Transforms that act on tapes
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -75,14 +84,17 @@
     ~qfunc_transform
     ~transforms.make_tape
 """
+# Import the decorators first to prevent circular imports when used in other transforms
+from .qfunc_transforms import make_tape, single_tape_transform, qfunc_transform
 from .adjoint import adjoint
 from .classical_jacobian import classical_jacobian
 from .control import ControlledOperation, ctrl
+from .decompositions import zyz_decomposition
 from .draw import draw
 from .hamiltonian_expand import hamiltonian_expand
 from .invisible import invisible
 from .measurement_grouping import measurement_grouping
 from .metric_tensor import metric_tensor, metric_tensor_tape
 from .specs import specs
-from .qfunc_transforms import make_tape, single_tape_transform, qfunc_transform
 from .qmc import apply_controlled_Q, quantum_monte_carlo
+from .unitary_to_rot import unitary_to_rot

pennylane/transforms/decompositions/__init__.py

@@ -0,0 +1,18 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+r"""This module contains decompositions for (numerically-specified) arbitrary
+unitary operations into sequences of elementary operations.
+"""
+
+from .single_qubit_unitary import zyz_decomposition

pennylane/transforms/decompositions/single_qubit_unitary.py

@@ -0,0 +1,105 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Contains transforms and helpers functions for decomposing arbitrary unitary
+operations into elementary gates.
+"""
+import pennylane as qml
+from pennylane import math
+
+
+def _convert_to_su2(U):
+    r"""Check unitarity of a matrix and convert it to :math:`SU(2)` if possible.
+
+    Args:
+        U (array[complex]): A matrix, presumed to be :math:`2 \times 2` and unitary.
+
+    Returns:
+        array[complex]: A :math:`2 \times 2` matrix in :math:`SU(2)` that is
+        equivalent to U up to a global phase.
+    """
+    # Check unitarity
+    if not math.allclose(math.dot(U, math.T(math.conj(U))), math.eye(2), atol=1e-7):
+        raise ValueError("Operator must be unitary.")
+
+    # Compute the determinant
+    det = U[0, 0] * U[1, 1] - U[0, 1] * U[1, 0]
+
+    # Convert to SU(2) if it's not close to 1
+    if not math.allclose(det, [1.0]):
+        exp_angle = -1j * math.cast_like(math.angle(det), 1j) / 2
+        U = math.cast_like(U, exp_angle) * math.exp(exp_angle)
+
+    return U
+
+
+def zyz_decomposition(U, wire):
+    r"""Recover the decomposition of a single-qubit matrix :math:`U` in terms of
+    elementary operations.
+
+    Diagonal operations will be converted to a single :class:`.RZ` gate, while non-diagonal
+    operations will be converted to a :class:`.Rot` gate that implements the original operation
+    up to a global phase in the form :math:`RZ(\omega) RY(\theta) RZ(\phi)`.
+
+    Args:
+        U (tensor): A 2 x 2 unitary matrix.
+        wire (Union[Wires, Sequence[int] or int]): The wire on which to apply the operation.
+
+    Returns:
+        list[qml.Operation]: A ``Rot`` gate on the specified wire that implements ``U``
+        up to a global phase, or an equivalent ``RZ`` gate if ``U`` is diagonal.
+
+    **Example**
+
+    Suppose we would like to apply the following unitary operation:
+
+    .. code-block:: python3
+
+        U = np.array([
+            [-0.28829348-0.78829734j,  0.30364367+0.45085995j],
+            [ 0.53396245-0.10177564j,  0.76279558-0.35024096j]
+        ])
+
+    For PennyLane devices that cannot natively implement ``QubitUnitary``, we
+    can instead recover a ``Rot`` gate that implements the same operation, up
+    to a global phase:
+
+    >>> decomp = zyz_decomposition(U, 0)
+    >>> decomp
+    [Rot(-0.24209529417800013, 1.14938178234275, 1.7330581433950871, wires=[0])]
+    """
+    U = _convert_to_su2(U)
+
+    # Check if the matrix is diagonal; only need to check one corner.
+    # If it is diagonal, we don't need a full Rot, just return an RZ.
+    if math.allclose(U[0, 1], [0.0]):
+        omega = 2 * math.angle(U[1, 1])
+        return [qml.RZ(omega, wires=wire)]
+
+    # If not diagonal, compute the angle of the RY
+    cos2_theta_over_2 = math.abs(U[0, 0] * U[1, 1])
+    theta = 2 * math.arccos(math.sqrt(cos2_theta_over_2))
+
+    # If the top left element is 0, can only use the off-diagonal elements We
+    # have to be very careful with the math here to ensure things that get
+    # multiplied together are of the correct type in the different interfaces.
+    if math.allclose(U[0, 0], [0.0]):
+        phi = 1j * math.log(U[0, 1] / U[1, 0])
+        omega = -phi - math.cast_like(2 * math.angle(U[1, 0]), phi)
+    else:
+        el_division = U[0, 0] / U[1, 0]
+        tan_part = math.cast_like(math.tan(theta / 2), el_division)
+        omega = 1j * math.log(tan_part * el_division)
+        phi = -omega - math.cast_like(2 * math.angle(U[0, 0]), omega)
+
+    return [qml.Rot(math.real(phi), math.real(theta), math.real(omega), wires=wire)]

pennylane/transforms/unitary_to_rot.py

@@ -0,0 +1,82 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+A transform for decomposing arbitrary single-qubit QubitUnitary gates into elementary gates.
+"""
+import pennylane as qml
+from pennylane.transforms import qfunc_transform
+from pennylane.transforms.decompositions import zyz_decomposition
+
+
+@qfunc_transform
+def unitary_to_rot(tape):
+    """Quantum function transform to decomposes all instances of single-qubit :class:`~.QubitUnitary`
+    operations to parametrized single-qubit operations.
+
+    Diagonal operations will be converted to a single :class:`.RZ` gate, while non-diagonal
+    operations will be converted to a :class:`.Rot` gate that implements the original operation
+    up to a global phase.
+
+    Args:
+        qfunc (function): a quantum function
+
+    **Example**
+
+    Suppose we would like to apply the following unitary operation:
+
+    .. code-block:: python3
+
+        U = np.array([
+            [-0.17111489+0.58564875j, -0.69352236-0.38309524j],
+            [ 0.25053735+0.75164238j,  0.60700543-0.06171855j]
+        ])
+
+    The ``unitary_to_rot`` transform enables us to decompose
+    such numerical operations (as well as unitaries that may be defined by parameters
+    within the QNode, and instantiated therein), while preserving differentiability.
+
+
+    .. code-block:: python3
+
+        def qfunc():
+            qml.QubitUnitary(U, wires=0)
+            return qml.expval(qml.PauliZ(0))
+
+    The original circuit is:
+
+    >>> dev = qml.device('default.qubit', wires=1)
+    >>> qnode = qml.QNode(qfunc, dev)
+    >>> print(qml.draw(qnode)())
+     0: ──U0──┤ ⟨Z⟩
+    U0 =
+    [[-0.17111489+0.58564875j -0.69352236-0.38309524j]
+     [ 0.25053735+0.75164238j  0.60700543-0.06171855j]]
+
+    We can use the transform to decompose the gate:
+
+    >>> transformed_qfunc = unitary_to_rot(qfunc)
+    >>> transformed_qnode = qml.QNode(transformed_qfunc, dev)
+    >>> print(qml.draw(transformed_qnode)())
+     0: ──Rot(-1.35, 1.83, -0.606)──┤ ⟨Z⟩
+
+    """
+    for op in tape.operations + tape.measurements:
+        if isinstance(op, qml.QubitUnitary):
+            # Only act on single-qubit unitary operations
+            if qml.math.shape(op.parameters[0]) != (2, 2):
+                continue
+
+            zyz_decomposition(op.parameters[0], op.wires[0])
+        else:
+            qml.apply(op)

tests/circuit_graph/test_qasm.py

@@ -227,12 +227,12 @@ def test_basis_state_initialization_decomposition(self):
     def test_unsupported_gate(self):
         """Test an exception is raised if an unsupported operation is
         applied."""
-        U = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
-
         with qml.tape.QuantumTape() as circuit:
-            qml.S(wires=0), qml.QubitUnitary(U, wires=[0, 1])
+            qml.S(wires=0), qml.DoubleExcitationPlus(0.5, wires=[0, 1, 2, 3])
 
-        with pytest.raises(ValueError, match="QubitUnitary not supported by the QASM serializer"):
+        with pytest.raises(
+            ValueError, match="DoubleExcitationPlus not supported by the QASM serializer"
+        ):
             res = circuit.to_openqasm()
 
     def test_rotations(self):
@@ -578,18 +578,19 @@ def qnode(state=None):
     def test_unsupported_gate(self):
         """Test an exception is raised if an unsupported operation is
         applied."""
-        U = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
-        dev = qml.device("default.qubit", wires=1)
+        dev = qml.device("default.qubit", wires=4)
 
         @qml.qnode(dev)
         def qnode():
             qml.S(wires=0)
-            qml.QubitUnitary(U, wires=0)
+            qml.DoubleExcitationPlus(0.5, wires=[0, 1, 2, 3])
             return qml.expval(qml.PauliZ(0))
 
         qnode()
 
-        with pytest.raises(ValueError, match="QubitUnitary not supported by the QASM serializer"):
+        with pytest.raises(
+            ValueError, match="DoubleExcitationPlus not supported by the QASM serializer"
+        ):
             qnode.qtape.to_openqasm()
 
     def test_rotations(self):