Add a transform get_unitary_matrix to obtain matrix representations of circuits

[glassnotes]

As suggested in #1458, it would be useful to have a transform to obtain the matrix representation of a quantum function. The transform should look and behave like the following pseudocode (similar to the draw or specs transform):

def get_unitary_matrix(qnode, wire_order=None):

    @wraps(qnode)
    def wrapper(*args, **kwargs):
        qnode.construct(args, kwargs)
     
        mat = qml.math.eye(2 ** num_wires)

        for op in qnode.qtape.operations:
            next_matrix = # construct matrix from op
            mat = qml.math.dot(next_matrix, mat)
    
        return mat

    return wrapper

This transform should be able to process all operations, including controlled gates on non-adjacent wires or wires in reverse order.

Note: this was written as a QNode transform because wire order is important (and we could extract that from the device). However, it could also be written as a tape or quantum function transform where the wire order must be passed explicitly every time.

The functions compute_matrix_from_ops_one_qubit and compute_matrix_from_ops_two_qubit contain some simple starting code. Once the new transform is implemented, usage of these should be replaced by the transform, as well as matrix computation in other parts of the code (such as SingleExcitation and DoubleExcitation tests).

---

List of tasks for completion:

1. Create a quantum function transform that takes the device wires as an argument;
2. Create a QNode transform that calls the quantum function transform;
3. Replace the usage of compute_matrix_from_ops_one_qubit and compute_matrix_from_ops_two_qubit;
4. Update the tests for the decomposition of the following operations to a new transform:

• ControlledPhaseShift
• SingleExcitationPlus
• SingleExcitationMinus
• DoubleExcitation

[josh146]

However it could also be written as a tape or quantum function transform where the wire order must be passed explicitly every time.

Oh! Actually, should we have both? The QNode function can just call this qfunc transform and pass the device wires as an argument.

This feels more flexible, since then you can use it even if you don't have a device

[glassnotes]

Yes, great idea! 💯

[ingstra]

I'm interested to look into this :)

[glassnotes]

@ingstra fantastic 🌟 Please feel free to ask us any questions you might have!

[ingstra]

Which file should the transforms be defined in?

Oh and another question :)
Do I need to take into account the possibility of multiqubit gates with an arbitrary number of qubits involved in the gate? Single qubit gates will obviously not be so hard to do this for, the two-qubit gates should also be ok after some thinking, then in the list of operations I see Toffoli so three-qubit gates have to work as well. Right now I'm thinking to look at the list of wires for a gate (e.g. qml.CNOT(wires=[0,2]) and construct the appropriate matrix case by case given how many wires are given, kind of like what's done in compute_matrix_from_ops_two_qubit. Would it be ok to throw a NotImplementedError if a list of 4 wires was given for a gate, for example?

[glassnotes]

Hi @ingstra ,

Which file should the transforms be defined in?

You can put it in a new file in the pennylane/transforms directory, with the same name as the transform itself.

Do I need to take into account the possibility of multiqubit gates with an arbitrary number of qubits involved in the gate?

Yes we'd like to have the transform work for arbitrary-sized gates. All of the existing gates do have matrix representations available regardless of size that you can use as a starting point (for example, you can do qml.Toffoli(wires=[0, 1, 2]).matrix, or qml.MultiControlledX(control_wires=[0, 1, 2, 3, 4, 5, 6], wires=7).matrix), it's just that the ordering of those matrix representations is fixed, and doesn't take into account the ordering of the wires on the device (which is why compute_matrix_from_ops_two_qubit has a wire_order argument).

For non-controlled gates, one way to work around this would be to apply SWAPs (but this could get computationally expensive). For controlled gates, there is a trick you can use to construct them using the projectors |0><0| = np.array([[1, 0], [0, 0]]) and |1><1| = np.array([[0, 0], [0, 1]]). A CNOT can be written as |0><0| \otimes I + |1><1| \otimes X. You can use this representation to compute the matrices for, e.g., CNOT(wires=[0, 2]) where wires are non-adjacent or in a different order, and it generalizes to arbitrary controlled operations and arbitrary number of controls.

[ingstra]

Gotcha, thanks. I got something that works for single-qubit gates now, but ran into a bit of a side-issue. qml.S.matrix seem to not give the matrix? If i do print(qml.S.matrix) it prints <property object at 0x7f80d714aef0> instead of the matrix. Same for the T and SX gates. Those that gives the matrix (Pauli, Hadamard) have defined the variable matrix and the function _matrix returns cls.matrix. The ones that don't give the matrix has no matrix variable and the function _matrix returns the numpy array. This is defined in non_parametric_ops.py I have no idea why this is happening?

[antalszava]

Hi @ingstra, indeed! The reason is that matrix is intended to be a property, so the expected behaviour for an operation is to call query matrix on an instance, e.g., qml.S(wires=[0]).matrix. However, as you mentioned, some operations might also have a class attribute called matrix (hence qml.CNOT.matrix returns the valid matrix). Having that class attribute is not per se a convention, but rather something that was introduced as a convenience for some operation classes. This is great to know about nonetheless, as it would be good to standardize this behaviour.

For now, could it work to use the matrix property of an instance for every operation or would that cause any challenges?

[ingstra]

I see! Thanks. Yes it works :)

[antalszava]

Created an issue to clear this up:
#1595

[ingstra]

Status update: I now have a quantum function transform that works for (has been tested on) single-qubit gates, two-qubit controlled gates, and Toffoli with adjacent or nonadjacent wires. My first idea was actually using the |0><0| \otimes I + |1><1| \otimes X kind of representation for controlled gates, but I ended up just doing SWAPs for everything. Now I need to test with MultiControlledX, and fix the other things on the to-do list. I'm just wondering if you'd like to see what I have so far, before I start editing previously existing code (tasks 3 & 4)? Maybe you can see something I missed or have other suggestions. I actually have one question for you guys (posed in a comment in the code ^_^)

One thing I was briefly thinking about was checking the is_symmetric_over_all_wires and is_symmetric_over_control_wires to do something to avoid doing superfluous SWAPs, but eh...

[ingstra]

I did some more testing today, and it turned out that it was just a fluke everything worked yesterday, but now I think it works (famous last words). Anyway, I started replacing compute_matrix_from_ops_one_qubit and compute_matrix_from_ops_two_qubit, and found that my method got a different result compared to compute_matrix_from_ops_two_qubit in function test_single_qubit_fusion_multiple_qubits in test_single_qubit_fusion.py. So I manually constructed the matrix, and it looks like I'm right and compute_matrix_from_ops_two_qubit is wrong o_o

[ingstra]

Sry for spamming but here's a question: I now have a quantum function transform named get_unitary_matrix. What should the QNode transform be named? Something different? The same? I learned from a quick Google search that Python doesn't support function overloading but I guess some type of workarounds exist.

[glassnotes]

Hi @ingstra , thanks for all your work on this! 🚀 I think you are correct that there is an issue in the current. compute_matrix_from_ops_two_qubit function. I'm pretty sure I know the cause (I added a swap, but did not swap back! The test cases it was used in must not have been affected). I will confirm tomorrow, and also provide some guidance for your other questions 🙂

BTW you can make a PR whenever you like, happy to quickly review it in an intermediate state (sometimes it's helpful to have just for pointing out where you might have questions).

[ingstra]

Thanks! I created a PR. In get_unitary matrix I copied some code from adjoint.py because it felt appropriate, but I don't really know how the tape stuff works so I don't know if all of it it is necessary/appropriate here.

Edit: Oh it seems to fail some tests. Interesting, it was fine when I just ran the tests locally. It seems to be some issue with the SequenceMatcher: TypeError: find_longest_match() missing 4 required positional arguments: 'alo', 'ahi', 'blo', and 'bhi'. Don't know why the error appears here but I don't get it locally?
Edit2: I found out why, it's the 3.7 test that fails, which makes sense since not having to have arguments was apparently implemented in 3.9 (and I use 3.9. That's what I get for not running the whole test suite...)

[glassnotes]

Thanks @ingstra , I will leave a few comments regarding your questions over in the PR!

[ingstra]

Right now this does not work with a circuit containing a controlled op generated by qml.ctrl, since it throws a NotImplementedError when get_unitary_matrix() tries to get op.matrix. Is this something I need to be concerned about at this point?

[ingstra]

I'm having some problems in test_qubit_ops.py for task 4, test_double_excitation_decomp: Something makes it not work when looping through the list of decompostitions. As a simple example, the following just prints the identity matrix two times. What is happening?

wires= [0,1]

list = [qml.RY(0.5,wires=[0]), qml.PauliX(wires=1)]

for op in list:
    def circ():
        op
        
    get_matrix = get_unitary_matrix(circ, wires)
    matrix = get_matrix()

    print(matrix)

[glassnotes]

Right now this does not work with a circuit containing a controlled op generated by qml.ctrl, since it throws a NotImplementedError when get_unitary_matrix() tries to get op.matrix. Is this something I need to be concerned about at this point?

In this case since the matrix is legit not available, I'd say don't worry about it. (This should maybe be written up as a separate issue, actually.)

I'm having some problems in test_qubit_ops.py for task 4, test_double_excitation_decomp: Something makes it not work when looping through the list of decompostitions.

I would need a bit more context, but one thing to note is that PennyLane operators are queued to the underlying tape upon construction, so you cannot apply an already-constructed operation as is done above. You'll need to use the qml.apply() function, so:

for op in list:
    def circ():
        qml.apply(op)

The identity is coming out because unless the operations are manually applied like this, nothing will be queued. This might solve the problem; let me know if it doesn't (seeing the full function would be helpful here).

[ingstra]

@glassnotes qml.apply() did the trick, thanks!

[glassnotes]

Closing as #1609 has been merged in 🎉