`YGate.power(1/2)` behaves differently with scipy 1.14 on macOS

[t-imamichi]

Environment

• Qiskit version: 1.2.4
• Python version: 3.12.7
• Operating system: macOS 15.0.1

What is happening?

YGate.power(1/2) behaves differently with scipy 1.14 on macOS compared with outcomes with other version and OS.

How can we reproduce the issue?

import scipy
from qiskit.circuit.library import YGate

print('scipy', scipy.__version__)
print('Y^0.5', YGate().power(1 / 2))
print('Y^0.5 †', YGate().power(1 / 2).adjoint())

macOS + qiskit 1.2.4 + scipy 1.14.1

scipy 1.14.1
Y^0.5 Instruction(name='unitary', num_qubits=1, num_clbits=0, params=[array([[ 0.5-0.5j,  0.5-0.5j],
       [-0.5+0.5j,  0.5-0.5j]])])
Y^0.5 † Instruction(name='unitary', num_qubits=1, num_clbits=0, params=[array([[ 0.5+0.5j, -0.5-0.5j],
       [ 0.5+0.5j,  0.5+0.5j]])])

macOS + qiskit 1.2.4 + scipy 1.14.0

scipy 1.14.0
Y^0.5 Instruction(name='unitary', num_qubits=1, num_clbits=0, params=[array([[ 0.5-0.5j,  0.5-0.5j],
       [-0.5+0.5j,  0.5-0.5j]])])
Y^0.5 † Instruction(name='unitary', num_qubits=1, num_clbits=0, params=[array([[ 0.5+0.5j, -0.5-0.5j],
       [ 0.5+0.5j,  0.5+0.5j]])])

macOS + qiskit 1.2.4 + scipy 1.13.1

scipy 1.13.1
Y^0.5 Instruction(name='unitary', num_qubits=1, num_clbits=0, params=[array([[ 0.5+0.5j, -0.5-0.5j],
       [ 0.5+0.5j,  0.5+0.5j]])])
Y^0.5 † Instruction(name='unitary', num_qubits=1, num_clbits=0, params=[array([[ 0.5-0.5j,  0.5-0.5j],
       [-0.5+0.5j,  0.5-0.5j]])])

RedHat 9.4 + qiskit 1.2.4 + scipy 1.14.1

scipy 1.14.1
Y^0.5 Instruction(name='unitary', num_qubits=1, num_clbits=0, params=[array([[ 0.5+0.5j, -0.5-0.5j],
       [ 0.5+0.5j,  0.5+0.5j]])])
Y^0.5 † Instruction(name='unitary', num_qubits=1, num_clbits=0, params=[array([[ 0.5-0.5j,  0.5-0.5j],
       [-0.5+0.5j,  0.5-0.5j]])])

What should happen?

Consistent outcomes regardless of scipy version

Any suggestions?

pin scipy version, e.g., scipy<1.14 for macOS?

[jakelishman]

Hmm, this is kind of weird - I don't think anything should have changed between the 1.2 series and main to change this. At the end of the day, this code should just defer to Scipy to do the decomposition.

Are you using the same versions of Scipy between the macOS 1.2.4 and main tests? Are they both running on the same processor? I can't reproduce on Intel macOS with scipy 1.13.1

[jakelishman]

Ah, I see you just updated to say it's a scipy 1.14 problem - thanks.

[t-imamichi]

I updated the description just now. I confirmed that it is due to scipy 1.14.

[jakelishman]

Ok, so Scipy 1.14 started using Accelerate as its default BLAS on macOS. This is quite possibly the difference - it used OpenBLAS before. I need to look a little bit into where something's going wrong, because we want to return the principal square root. I'm hoping we weren't just doing that by chance previously.

[t-imamichi]

Thank you for investigating this issue so quickly. Using Accelerate is a big change.

[jakelishman]

Ok, so digging a little into the Operator code, I feel like we might actually have some problems going on here, and Accelerate's behaviour is not ideal for us, but not wrong.

Let's handle the immediate problem first: at the end of the day, the trouble is that Accelerate is finding that the eigenvalues of $Y$ to be $(1, -(1 + \text{eps}\times i))$ where eps is ~$10^{-32}$. The OpenBLAS implementation found them to be 1 and -1 exactly. This unfortunately means that when we then take powers of the diagonal matrix, Accelerate's negative eigenvalue happens to fall on the other side of the principal-root branch cut.

Now, an additional problem: 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:

>>> 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 ]]

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'll need to fix Operator.power for non-unitary operators, but that's a separate issue to this one.

Unfortunately for us, scipy.linalg.fractional_matrix_power has the same trouble with the branch cut between OpenBLAS and Accelerate (unsurprisingly - it's using a Schur-based method internally, and I wouldn't be surprised if the difference comes in in the QR step, which is the base of basically everything). I don't immediately have a way to mitigate this branch cut effect that's actually going to be reliable - it's an annoying common problem in computational linear algebra, but I'm not sure I've got enough experience to see the best way out of the situation.

[jakelishman]

Fwiw, while we're seeing the problem with Accelerate here with the default Scipy 1.14 wheels, in principle it can happen with any Scipy version if the user switches out the BLAS implementation it's built against. So pinning Scipy would just unfortunately squashing a single symptom.

If we can't come up with a fix to this that at least makes it less surprisingly, I guess we'll at least need to heavily document in Operator.power and Gate.power that the results of fractional powers may well fall afoul of branch-cuts in the complex plane depending on the BLAS version, available SIMD instructions, etc.