How to_unitary_mut should work

[bramathon]

Currently, to_unitary_mut requires the user to specify the dimension of the unitary.

This seem unnecessary and I'm not sure what it's used for. I suspect, however, that it's used to "lift" unitaries to a larger Hilbert space based on the qubit indices. This is not a good approach because it assumes that the space is for qubits 0 to n, in that order while in fact a Hilbert space could be defined on qubits 5, 3, 10, 33, for example. If lifting is required, a separate function should be provided and it should require the Hilbert space qubit order to be specified.

Rather to_unitary_mut should return the canonical matrix of the gate. There is no need to specify the dimension, it's a property of the gate and should produce a matrix of dimensions (2^n, 2^n).

Furthermore, I would like all Gate instructions in a program to have a unitary attached to them. For standard quil gates, this is achieved by to_unitary_mut, however, for defined gates in the program this is not the case. Rather the user needs to check the defined gates of the program to know what is the definition of an instruction. Ideally, this would not be necessary, the instruction would have the defgate associated with it allowing the following.

program = """DEFGATE BLARG AS MATRIX:
    -0.28705672572498936-0.11277532694223659i, -0.5998793485470028-0.3613248187346551i, -0.05259662196286458-0.6195952054027546i, 0.05949628004150796+0.15577186480089364i
    0.0042693221909431-0.37767648756583294i, 0.21434622020589333+0.0861879617069231i, -0.3702023023607123-0.3181064164720223i, -0.6471484586315557-0.3833105951903434i
    -0.3999605210419518+0.6267674838894352i, 0.16961228281142476+0.24120155674611812i, 0.3866424148211197-0.30283186464280754i, -0.344791688723465-0.012914391002983222i
    -0.407825281226674+0.2075934686261242i, 0.4307462151034142-0.42863464210242547i, -0.36409272342149845-0.03648196124031003i, 0.39962918525837543-0.35737319861073946i

BLARG 0 1"""

program = Program.parse(program)

program.body_instructions[0].to_unitary_mut()

[asaites]

Thank you for opening this issue. Indeed, the purpose of to_unitary is to "lift a Gate to the full n_qubits-qubit Hilbert space". The implementation appears to be a port of the PyQuil lifted_gate function, which is arguably a more appropriate name. The documentation there may give better context for your primary question.

The method was added to quil-py to support getting the matrix represented by a particular Program (provided that it consists only of Gate instructions with known parameters). That may help explain the n_qubits parameter.

In your comment, you said:

[I]t assumes that the space is for qubits 0 to n, in that order while in fact a Hilbert space could be defined on qubits 5, 3, 10, 33, for example.

The comments on those functions indicate it should work with arbitrary qubit ordering, and you can see there's work done to construct a permutation matrix and use it to modify the result. Unfortunately, here's what happens when trying the example described in the function's comments:

>>> from quil.instructions import Instruction
>>> Instruction.parse("CCNOT 20 15 10").inner().to_unitary_mut(3)
[...]
pyo3_runtime.PanicException: These arrays cover the same range.

The same basic problem exists with the function in PyQuil:

>>> from pyquil.gates import Gate
>>> from pyquil.simulation.tools import lifted_gate
>>> g = lifted_gate(Gate("CCNOT", [], [Qubit(n) for n in (20,15,10)], []), 3)
[...]
IndexError: index 0 is out of bounds for axis 0 with size 0

My guess is there's a bug in part of the code that generates the permutation matrix, and was ported when the feature was added to quil-rs. The same basic bug exists in PyQuil. As far as I can tell, it was even in the original implementation (which made it to PyQuil as part of this PR).

It's quite surprising the bug has been there nearly a decade without surfacing. My guess is that any uses of this method must be indexing qubits from $0..n$.

So at the very least, for this method/the PyQuil function to be valuable, this bug needs to be fixed. Doing so would likely need to happen with fixing #464.

---

In a related vein, it's become apparent that Gates really must be immutable for Python users (as we want them to be hashable), so if this particular method is to remain, it needs to become immutable. But we could get rid of it while still supporting the Program case described above. In either case, it likely makes sense to have some way to a matrix from a Gate, particularly those which take parameters -- or otherwise to have some way of easily obtaining a numpy array equivalent to the standard gates.

[bramathon]

Yes, the pyquil state-vector machinery, which is where lifted_gate comes from, assumes qubits indexed from 0 to n. I'm not sure this is a bug, I think it was made that way on purpose.

[bramathon]

My point was a little more subtle though. If you have a program

CCNOT 20 15 10
CCNOT 15 10 20

and you used this method, both gates would return the same unitary and if you multiplied those unitary to fund the program unitary you would get the wrong answer. The lifting operation is missing information - the order of the qubits in the Hilbert space. If you force your qubits to be indexed 0 to n, that problem is resolved.

[asaites]

The lifting operation is missing information - the order of the qubits in the Hilbert space.

If I understand you and the code correctly, the current code intends to use the indexes from the Gate.qubits to build a permutation matrix, which it applies to the gate matrix. Does this work as you expect (note the order 2 0 1 vs 2 1 0)?

>>> Instruction.parse("CCNOT 2 0 1").inner().to_unitary_mut(3)
array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j]])
>>> Instruction.parse("CCNOT 2 1 0").inner().to_unitary_mut(3)
array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j]])

I think that works as intended, but using higher numbered qubits will only work if n_qubits is at least 1+max(q.inner() for q in gate.qubits):

>>> Instruction.parse("CNOT 1 2").inner().to_unitary_mut(2)
[...]
pyo3_runtime.PanicException: ndarray: inputs 0 × 0 and 4 × 4 are not compatible for matrix multiplication

I'm pretty sure the specific error is a red herring, and that this actually failed due to #464. Smaller numbers fail within the permutation code:

>>> Instruction.parse("CNOT 1 2").inner().to_unitary_mut(1)
[...]
pyo3_runtime.PanicException: These arrays cover the same range.
>>> Instruction.parse("CNOT 1 2").inner().to_unitary_mut(0)
[...]
pyo3_runtime.PanicException: These arrays cover the same range.
>>>
>>> lifted_gate(Gate._from_rs_gate(Instruction.parse("CNOT 1 2").inner()), 3)
array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]])
>>> lifted_gate(Gate._from_rs_gate(Instruction.parse("CNOT 1 2").inner()), 2)
[...]
ValueError: Permutation SWAP index not valid
>>> lifted_gate(Gate._from_rs_gate(Instruction.parse("CNOT 1 2").inner()), 1)
[...]
IndexError: index 0 is out of bounds for axis 0 with size 0

[antalsz]

@bramathon Can you clarify what you're using to_unitary_mut for?

[bramathon]

We need to know the unitary for each gate in order to perform simulation and compilation passes. We currently use a python function to do this, but it only supports matrix based gates.

[bramathon]

Does this work as you expect (note the order 2 0 1 vs 2 1 0)?

This works as expected (now that I understand it's a lifted_gate) - the CCNOT unitaries are lifted to the Hilbert space with the qubit order 0 1 2. However, it ought to work by simply returning the canonical CCNOT unitary, regardless of what the target qubits are:

        [1, 0, 0, 0, 0, 0, 0, 0],
        [0, 1, 0, 0, 0, 0, 0, 0],
        [0, 0, 1, 0, 0, 0, 0, 0],
        [0, 0, 0, 1, 0, 0, 0, 0],
        [0, 0, 0, 0, 1, 0, 0, 0],
        [0, 0, 0, 0, 0, 1, 0, 0],
        [0, 0, 0, 0, 0, 0, 0, 1],
        [0, 0, 0, 0, 0, 0, 1, 0],

This is necessary because users don't always use qubits 0...n, nor do we always want the state to be in that order. Rather the order of the Hilbert space should be specified. CCNOT should simply represent the CCNOT matrix and the simulator or compiler will apply it how it sees fit.

[bramathon]

I think the best way to achieve this is as follows:

1. Programs always include defgates. The std gate set is included, and these gates cannot be redefined.
2. Each Gate has a unitary property which is callable.
3. Gates are created with the unitary property set to raise an error.
4. When a gate is added to a program, the unitary property is set according to the defgates of the program.
5. When a defgate is added to a program, the unitary property of all gates in the program is also modified.
6. The unitary method of parametric gates accepts parameters.