Adds elementary gate decomposition for single-qubit `QubitUnitary` #1427

glassnotes · 2021-06-22T20:18:23Z

Context: Currently, there is no decomposition/synthesis for the arbitrary QubitUnitary operation.

Description of the Change: Adds the math needed to express arbitrary single-qubit matrices as a Rot gate (or an RZ gate, if they are diagonal).

Benefits: This should allow QubitUnitary to run on more devices. Furthermore, I think that by using the decomposition we can make QubitUnitary fully differentiable, at least for the single-qubit case (see example below).

Possible Drawbacks: None

Related GitHub Issues: None

Example:

import pennylane as qml
from pennylane import numpy as np

dev = qml.device('default.qubit', wires=1)

def qfunc(x):
    qml.Hadamard(wires=0)

    # Create a unitary matrix based on an input parameter
    Ux = np.array([[1, 0], [0, np.exp(1j * x)]])
    ops = qml.QubitUnitary.decomposition(Ux, wires=0)
    for op in ops:
        op.queue()

    return qml.expval(qml.PauliX(0))

>>> original_qnode = qml.QNode(qfunc, dev, diff_method="adjoint")
>>> print(qml.grad(original_qnode)(0.2))
-0.19866933079506127

github-actions · 2021-06-22T20:18:44Z

Hello. You may have forgotten to update the changelog!
Please edit .github/CHANGELOG.md with:

A one-to-two sentence description of the change. You may include a small working example for new features.
A link back to this PR.
Your name (or GitHub username) in the contributors section.

josh146

@glassnotes this is awesome! A couple of things that crossed my mind:

This means that we can now support qml.Hermitian (arbitrary Hermitian matrix observables) directly on hardware, since it currently diagonalizes down to a QubitUnitary rotation 🙌
I love that it results in differentiability 🤯
In your code block in the PR comment, I don't think you need the op.queue()? Calling .decomposition should automatically queue the ops, so they might be being queued twice.
Part of me wonder if this feels like a compilation transform? Becuase there will be use cases where users want to force a unitary decomposition on a device that natively supports it. But I don't know? it would add overhead compared to simply calling QubitUnitary.decomposition() manually.
Don't forget to add interface tests! To ensure it works with different tensor input
Along the same lines, gradient tests should also be added 🙂

pennylane/ops/qubit.py

josh146 · 2021-06-23T08:26:15Z

pennylane/ops/qubit.py

+            # if desired. If the top left element is 0, can only use the off-diagonal elements
+            if qml.math.isclose(U[0, 0], 0):
+                phi = 1j * qml.math.log(U[0, 1] / U[1, 0])
+                omega = -phi - 2 * qml.math.angle(U[1, 0])


You may need to be careful here:

phi is complex

2 * angle(...) will be real (I think float64)

and TensorFlow (maybe also PyTorch?) will not automatically cast the latter to a complex number before doing the subtraction.

I think the following should work:

Suggested change

omega = -phi - 2 * qml.math.angle(U[1, 0])

omega = -phi - qml.math.cast_like(2 * qml.math.angle(U[1, 0]), phi)

(Note: qml.math.cast(..., "complex128") would also work above, but cast_like is safer, in that it will support both complex64 and complex128 depending on whatever phi is)

Oh thanks, yes this looks safer! I originally had the conversion to real in the line of phi rather than in the return value, do you think it'd be better to convert phi first so that the subtraction ends up real? Or is it better to just do the casting anyways?

I think you'd need the casting regardless, potentially you could have float64 subtract float32, which would also cause an error

josh146 · 2021-06-23T08:27:33Z

pennylane/ops/qubit.py

+                omega = -phi - 2 * qml.math.angle(U[1, 0])
+            else:
+                omega = 1j * qml.math.log(qml.math.tan(theta / 2) * U[0, 0] / U[1, 0])
+                phi = -omega - 2 * qml.math.angle(U[0, 0])


Suggested change

phi = -omega - 2 * qml.math.angle(U[0, 0])

phi = -omega - qml.math.cast_like(2 * qml.math.angle(U[0, 0]), omega)

tests/ops/test_qubit_ops.py

Co-authored-by: Josh Izaac <josh146@gmail.com>

glassnotes · 2021-06-23T13:19:55Z

Part of me wonder if this feels like a compilation transform? Becuase there will be use cases where users want to force a unitary decomposition on a device that natively supports it. But I don't know? it would add overhead compared to simply calling QubitUnitary.decomposition() manually.

I am wondering this too; so, for additional context, I'd like to do this for the 2-qubit case as well. That's going to make the decomposition method extremely hairy. If we use transforms, we can make things more flexible (e.g., we can have different transforms for ZYZ versus ZXZ decompositions), and better separate the 1- and 2-qubit cases from each other. Then we can have something that passes through the circuit, finds QubitUnitarys, and depending on the number of qubits calls the appropriate transform.

josh146 · 2021-06-23T13:28:31Z

I am wondering this too; so, for additional context, I'd like to do this for the 2-qubit case as well.

Oh! Maybe the answer is... both?

We have transforms to make it flexible and keep the code organized.
The QubitUnitary.decomposition() method calls the transforms as needed, in order to support devices where QubitUnitary is not defined.

glassnotes · 2021-06-24T12:36:42Z

* We have transforms to make it flexible and keep the code organized.

* The `QubitUnitary.decomposition()` method calls the transforms as needed, in order to support devices where `QubitUnitary` is not defined.

I like it! I'll give it a go and we can see what it looks like. After deciding on that I'll make the rest of the tests.

…ngle.

josh146

I really like the new organization! Especially the fact that the decompositions can be done three ways:

Automatically if the device does not support it
Manually by calling qml.QubitUnitary.decomposition(U)
Manually on an existing qfunc via a transform

josh146 · 2021-06-28T09:00:02Z

pennylane/ops/qubit.py

+    @staticmethod
+    def decomposition(U, wires):
+        # Decompose arbitrary single-qubit unitaries as the form RZ RY RZ
+        if qml.math.shape(U)[0] == 2:
+            wire = Wires(wires)[0]
+            decomp_ops = qml.transforms.decompositions._zyz_decomposition(U, wire)
+            return decomp_ops
+
+        return NotImplementedError("Decompositions only supported for single-qubit unitaries")


This is really nice! In particular, I like how the logic is in the transforms package in its own module. This will make it much easier to modify and test for new contributors, rather than looking through the already large qubit.py file.

Yes, thanks for the suggestion to separate it out! 😁 It will also make it much cleaner to add different decompositions, one for 2-qubit unitaries, etc. while keeping the qubit ops file tidy.

pennylane/transforms/__init__.py