-
Notifications
You must be signed in to change notification settings - Fork 586
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Decomposition of arbitrary two-qubit unitaries. #1552
Decomposition of arbitrary two-qubit unitaries. #1552
Conversation
Hello. You may have forgotten to update the changelog!
|
@josh146 this is ready for a first look. I am just having trouble getting the differentiability working. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is an awesome new feature @glassnotes!
pennylane/ops/qubit/matrix_ops.py
Outdated
if qml.math.shape(U) == (4, 4): | ||
wires = Wires(wires) | ||
decomp_ops = qml.transforms.decompositions.two_qubit_decomposition(U, wires) | ||
return decomp_ops |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Nice 👍 This is super cool.
I must admit though, after writing the tracing ADR, everywhere I look I see things that won't be supported by tracing 😆 If we decide to go down that route, it means we will have to rethink a lot of the code base (or retrace any time argument shapes change).
|
||
# This gate E is called the "magic basis". It can be used to convert between | ||
# SO(4) and SU(2) x SU(2). For A in SO(4), E A E^\dag is in SU(2) x SU(2). | ||
E = qml.math.array([[1, 1j, 0, 0], [0, 0, 1j, 1], [0, 0, 1j, -1], [1, -1j, 0, 0]]) / np.sqrt(2) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have a feeling that qml.math.array()
simply dispatches to np.array
? Would it make more sense to just use NumPy here, so that it is clear that these are non-differentiable constants?
Then you would be able to do
E = np.array([[1, 1j, 0, 0], [0, 0, 1j, 1], [0, 0, 1j, -1], [1, -1j, 0, 0]]) / np.sqrt(2)
Et = E.T
Edag = Et.conj()
mat = qml.math.zeros((4, 4)) | ||
|
||
for row_idx in range(4): | ||
mat[row_idx, seq[row_idx]] = 1 | ||
|
||
return mat |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This can be vectorized for performance:
>>> seq = np.array([0, 1, 4, 3, 2])
>>> mat = np.identity(5)
>>> mat[:, seq]
array([[1., 0., 0., 0., 0.],
[0., 1., 0., 0., 0.],
[0., 0., 0., 0., 1.],
[0., 0., 0., 1., 0.],
[0., 0., 1., 0., 0.]])
I think this is also a case where you can just use plain NumPy?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh sweet - when streamlined like this, it doesn't even need to be in its own function. Thanks!
D_ops = zyz_decomposition(D, wires[1]) | ||
|
||
# Return the full decomposition | ||
return C_ops + D_ops + interior_decomp + A_ops + B_ops |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
🥇
assert check_matrix_equivalence(U, obtained_matrix, atol=1e-7) | ||
|
||
|
||
class TestTwoQubitUnitaryDecompositionInterfaces: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It might be the case that gradient tests, when added here, might expose issues (e.g., the array assignment)
assert qml.math.allclose(original_grad, transformed_grad, atol=1e-6) | ||
|
||
@pytest.mark.parametrize("U", test_u4_unitaries) | ||
def test_gradient_unitary_to_rot_torch_two_qubit(self, U): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
From looking at the CI, it seems that the (current) reasons Torch is failing is qml.math.linalg.multi_dot([Edag, U, E])
. It seems that Torch has no equivalent to this, so qml.math
is seeing a list input and simply attempting to dispatch to vanilla NumPy :(
This is probably a case of multi-dispatch - we need to look at the types of Edag, U, and E, and determine which interface to dispatch to based on that. So we have two options:
- Add
qml.math.multi_dot
tomulti_dispatch.py
- Chain
qml.math.dot
instead
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The first one seems like the nicer solution, though maybe something out of scope of this PR. For now I will chain dot
s just for the purpose of getting a differentiable version working, and then I will look into this 👍
Co-authored-by: Josh Izaac <josh146@gmail.com>
pennylane/ops/qubit/matrix_ops.py
Outdated
raise NotImplementedError("Decompositions only supported for single-qubit unitaries") | ||
if qml.math.shape(U) == (4, 4): | ||
wires = Wires(wires) | ||
decomp_ops = qml.transforms.decompositions.two_qubit_decomposition(U, wires) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@josh146 so doing this here actually messes up the behaviour when tape.expand()
is called, because of the order in which the operations are created (the centre part gets computed before the SU(2) stuff around). However if I call this with invisible
around it, nothing happens during expand
because nothing gets queue. It seems kind of weird to loop through the decomp_ops
and queue them here, is there a different workaround for this?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Huh, this is a really good question 🤔
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
if qml.math.shape(U) == (4, 4):
wires = Wires(wires)
decomp_ops = qml.transforms.invisible(
qml.transforms.two_qubit_decomposition
)(U, wires)
current_tape = qml.tape.get_active_tape()
if current_tape:
for op in decomp_ops:
qml.apply(op, current_tape)
return decomp_ops
Huh, this feels gross, but it works.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
My mind is going blank, all I can think about is
loop through the decomp_ops and queue them here
Also, it means that users won't be able to use qml.transforms.decompositions.two_qubit_decomposition(U, wires)
directly either, since the order will be wrong.
I suppose what you could do is call the utility functions invisibly, inside two_qubit_decomposition
, in order to generate a list of operators in the correct order. And then, right before returning, apply()
each op in order? This is exactly the same as what you suggest above, just 'inside' two_qubit_decomposition
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Huh, this feels gross, but it works.
😆
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
^ What you suggest is the nicest. It further simplifies unitary_to_rot
so that the 2-qubit case looks just like the 1-qubit case.
|
||
|
||
################################################################################### | ||
# Developer notes: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is the best info I can get right now; it doesn't seem like there is a straightforward fix so I will leave this to future work.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is super useful to have 💯
… of https://github.com/PennyLaneAI/pennylane into ch7103-implement-differentiable-two-qubit-qubitunitary
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
💯
@@ -88,7 +88,7 @@ def decomposition(U, wires): | |||
|
|||
if qml.math.shape(U) == (4, 4): | |||
wires = Wires(wires) | |||
decomp_ops = qml.transforms.decompositions.two_qubit_decomposition(U, wires) | |||
decomp_ops = qml.transforms.two_qubit_decomposition(U, wires) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
👨🍳
-╭U- = -A- | ||
-╰U- = -B- |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
oh! nice
-╭U- = -C--╭C--A- | ||
-╰U- = -D--╰X--B- |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
much easier to read!
# U depends on the input parameters | ||
qml.QubitUnitary(U, wires=["a", "b"]) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah, I think I get this now - U
cannot depend on trainable parameters, it must be constant
elif qml.math.shape(op.parameters[0]) == (4, 4): | ||
two_qubit_decomposition(op.parameters[0], op.wires) | ||
else: | ||
qml.apply(op) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is definitely covered in tests; the circuits run through the transform consist of more than just Never mind, I see the issue.QubitUnitary
operations.
Co-authored-by: Josh Izaac <josh146@gmail.com>
if not math.allclose(ev_p, ev_q): | ||
new_q_order = [] | ||
for _, eigval in enumerate(ev_p): | ||
are_close = [math.allclose(x, eigval) for x in ev_q] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I am not sure why these 3 lines are not showing as covered; the 2- and 3-CNOT test cases should all be going through this portion.
… of https://github.com/PennyLaneAI/pennylane into ch7103-implement-differentiable-two-qubit-qubitunitary
* sort qubit ops tests * remove imports * pr #1552 update * move into folder. update some docstrings * black and changelog * oops. * Apply suggestions from code review Co-authored-by: antalszava <antalszava@gmail.com> * respond to review * rename file * complete renaming Co-authored-by: antalszava <antalszava@gmail.com>
Context: Addition of compilation tools to PennyLane.
Description of the Change: Adds a decomposition method for arbitrary two-qubit unitary operations.
Benefits: Any two-qubit unitary in a circuit can now be decomposed into elementary operations, in all interfaces.
Possible Drawbacks: This method was very very challenging to implement. I was not able to get full differentiability (there are some notes in the
two_qubit_unitary.py
file about this). It is possible to take gradients when the two-qubit unitaries that are decomposed are constant (i.e., do not depend on the QNode parameters), but not when they do depend on the QNode parameters.Related GitHub Issues: None