GlobalPhaseGate generates an error when trying to generate Statevector

[diemilio]

Environment

• Qiskit Terra version: 0.24.0rc1
• Python version: 3.10
• Operating system: MacOS Monterrey

What is happening?

Statevector class generates a tuple index out of range error when a quantum circuit includes an instance of the recently-added GlobalPhaseGate.

How can we reproduce the issue?

from math import pi
from qiskit import QuantumCircuit
from qiskit.circuit.library import GlobalPhaseGate
from qiskit.quantum_info import Statevector

qc = QuantumCircuit(1)
qc.append(GlobalPhaseGate(pi/2))
print(Statevector(qc))

What should happen?

The issue seems to be related to how the operator used to evolve the initial state is being generated. This operator uses the instruction .to_matrix() method, which for the GlobalPhaseGate results in a matrix of shape (1,1), which is basically an array containing the global phase being added. For example:

>>> from qiskit.circuit.library import GlobalPhaseGate
>>> from math import pi
>>> GlobalPhaseGate(pi/3).to_matrix()
array([[0.5+0.8660254j]])

For this to work correctly, the operator needs to be a diagonal matrix with the global phase value of shape (2^n_qubits, 2^n_qubits) instead of just the value itself.

Any suggestions?

I would volunteer to help with this, but if I am not mistaken, the to_matrix() method is inherited from the Gate class, so I can't think of an easy solution other than treating GlobalPhaseGate as a special case.

[jakelishman]

Thanks for this report - there were bound to be some more issues with 0q operations that #8272 missed.

The matrix returned shouldn't be full shape (and can't - we don't pass that information to to_matrix). I think the return value as things stand is acceptable, since it's returning a 2d array of shape (2**n_qubits, 2**n_qubits). It feels more like a problem with the handling of scalar values within quantum_info, I think; the code that explodes here is because there's an operator that's of shape (1,), whereas it might want to be (1, 1). I haven't fully looked through the code / thought about things though.

edit: having looked a little further, OpShape.shape deliberately returns a one-element tuple when there's no qubits on the right-hand side so that it looks like a statevector. I suspect then that the proper fix is for it to detect the num_qargs_l = num_qargs_r = 0 case, and reduce that to a returned scalar shape of (). That will probably have some knock-on effects throughout the code-base, but it's probably the most correct thing for us to try.

[jakelishman]

To add on to my edit: the code in Statevector._evolve_operator would still need updating in that situation to handle the possibility of a scalar operator.

[diemilio]

Hi @jakelishman. Sorry, I just realized that saying that "GlobalPhaseGate might need to be treated as a special case" maybe caused some confusion. I didn't mean to say that the to_matrix() method should take as an input the size of the system, or that it shouldn't return something different than what it already does; sorry about that.

What I meant is that in Statevector._evolve_instruction the operator is generated using Operator._instruction_to_matrix(), which generates the matrix using the to_matrix() method of the instruction of interest. Since for the GlobalPhaseGate this is a (1, 1) matrix, the evolution is not working correctly because the operator should instead be of size (2^n, 2^n). So, somewhere along this chain, this needs to be fixed, perhaps exclusively for the GlobalPhaseGate.

One thing that it's worth pointing out, is that the correct matrix is generated when using the Operator class. For example:

>>> from math import pi
>>> from qiskit import QuantumCircuit
>>> from qiskit.quantum_info import Operator
>>> from qiskit.circuit.library import GlobalPhaseGate
>>> qc = QuantumCircuit(2)
>>> qc.append(GlobalPhaseGate(pi/3))
>>> print(Operator(qc))
Operator([[0.5+0.8660254j, 0. +0.j       , 0. +0.j       , 0. +0.j       ],
          [0. +0.j       , 0.5+0.8660254j, 0. +0.j       , 0. +0.j       ],
          [0. +0.j       , 0. +0.j       , 0.5+0.8660254j, 0. +0.j       ],
          [0. +0.j       , 0. +0.j       , 0. +0.j       , 0.5+0.8660254j]],
         input_dims=(2, 2), output_dims=(2, 2))

What I noticed is that when creating an Operator object, a matrix of the size of the system is first created, and then the overall operator is constructed using Operator.compose, which checks that the dimensions are compatible, generating the right output. So, something similar would have to be done when trying to generate a statevector.

[jakelishman]

My point is that the global phase gate is of the correct size already - it's a 0q operator, so it matches your definition, and Operator working correctly shows that. Statevector has slightly different handling of the OpShape, and that's what made me notice that the 0q case is being handled as "vector-like" instead of "operator-like". This coincidentally doesn't cause bugs with the construction of an Operator, but it is still part of the underlying issue.

[jakelishman]

I changed my mind again: I think actually treating the OpShape.shape as producing a (1, 1) value if both the left and right sizes are 1 (as opposed to producing a vector-like shape) is the cleanest way to go about this; the rest of the handling should "just work" if we do this.

If you'd like to make the PR, please feel free to, otherwise I can do it no trouble.

[diemilio]

@jakelishman I see what you mean now; sorry for the confusion.

I added a conditional for when both self.num_qargs_l and num_qargs_r are zero, and that fixed the problem. I could create a PR but not sure what's the most legible way to implement it (completely separate if statement for the 0q case?); so I think it would actually save you time if you do it instead of having to answer my questions :)

[jakelishman]

No worries, I've made #10031 that does the fix I think both you and I have described.