Operator.power does not handle non-unitary inputs correctly

[jakelishman]

Environment

• Qiskit version: 1.2.4
• Python version: 3.10
• Operating system: macOS

What is happening?

The way Operator.power works, it's strongly assuming that the internal matrix is unitary, because it's assuming that it can calculate fractional matrix powers by taking the Schur decomposition taking the matrix power of the Schur form then applying the unitary transformation (ok so far), except it assumes that the matrix power of the Schur form is the power of its diagonal. The Schur form is upper-triangular, but not necessarily diagonal (roughly, I think it's diagonal in the case that the input is a scaled unitary). So our Operator.power code fails for non-unitary operators.

How can we reproduce the issue?

>>> from qiskit.quantum_info import Operator
>>> import scipy.linalg
>>> # Not unitary.
>>> data = [[1, 1], [0, -1]]
>>> print(Operator(data).power(0.5).data)
[[1.000000e+00+0.j 0.000000e+00+0.j]
 [0.000000e+00+0.j 6.123234e-17+1.j]]
>>> print(scipy.linalg.fractional_matrix_power(data, 0.5)))
[[1.00000000e+00+0.j  5.00000000e-01-0.5j]
 [0.00000000e+00+0.j  4.93096709e-17+1.j ]]

What should happen?

The two should give (approximately) equal results (though see #13305 for branch-cut problems).

Any suggestions?

The ability to do fractional matrix powers of Operator is relatively recent (#11534), and it seems to have pulled the code from Gate.power, which was empowered to make that assumption.

We probably should just use scipy.linalg.fractional_matrix_power for Operator.power, and specialise Gate.power to use a more explicit eigenvector/eigenvalue form, since it's allowed to assume that its input was unitary.