Fix compatibility issues with SciPy 1.14

The main change here is that SciPy 1.14 on macOS now uses the system
Accelerate rather than a bundled OpenBLAS by default.  These have
different characteristics for several LAPACK drivers, which caused
numerical instability in our test suite.  Fundamentally, these problems
existed before; it was always possible to switch out the BLAS/LAPACK
implementation that SciPy used, but in practice, the vast majority of
users (and our CI) use the system defaults.

The modification to `Operator.power` to shift the branch cut was
suggested by Lev.  As a side-effect of how it's implemented, it fixes
an issue with `Operator.power` on non-unitary matrices, which Sasha had
been looking at.

The modification to the Choi-to-Kraus conversion reverts back from a
Schur-based decomposition to an `eigh`-based one.  This was noted in a
code comment that it was causing instability, and tracing the git
history through to fdd5603 (gh-3884) suggests that this was around
the time of Scipy 1.1 to 1.3, with OpenBLASes around 0.3.6.  The bundled
version of OpenBLAS with SciPy 1.13 was 0.3.27, so several releases on.
If `eigh` problems re-appear, we can look at changing the decomposition
back to something else, but the current `eigh` seems to be more stable
than the Schur form for a wider range of inputs now.

Co-authored-by: Lev S. Bishop <[email protected]>
Co-authored-by: Alexander Ivrii <[email protected]>

constraints.txt

@@ -3,10 +3,6 @@
 # https://github.com/Qiskit/qiskit-terra/issues/10345 for current details.
 scipy<1.11; python_version<'3.12'
 
-# Temporary pin to avoid CI issues caused by scipy 1.14.0
-# See https://github.com/Qiskit/qiskit/issues/12655 for current details.
-scipy==1.13.1; python_version=='3.12'
-
 # z3-solver from 4.12.3 onwards upped the minimum macOS API version for its
 # wheels to 11.7. The Azure VM images contain pre-built CPythons, of which at
 # least CPython 3.8 was compiled for an older macOS, so does not match a

qiskit/quantum_info/operators/channel/transformations.py

@@ -220,32 +220,28 @@ def _kraus_to_choi(data):
 
 def _choi_to_kraus(data, input_dim, output_dim, atol=ATOL_DEFAULT):
     """Transform Choi representation to Kraus representation."""
-    from scipy import linalg as la
+    import scipy.linalg
 
     # Check if hermitian matrix
     if is_hermitian_matrix(data, atol=atol):
-        # Get eigen-decomposition of Choi-matrix
-        # This should be a call to la.eigh, but there is an OpenBlas
-        # threading issue that is causing segfaults.
-        # Need schur here since la.eig does not
-        # guarantee orthogonality in degenerate subspaces
-        w, v = la.schur(data, output="complex")
-        w = w.diagonal().real
-        # Check eigenvalues are non-negative
-        if len(w[w < -atol]) == 0:
-            # CP-map Kraus representation
-            kraus = []
-            for val, vec in zip(w, v.T):
-                if abs(val) > atol:
-                    k = np.sqrt(val) * vec.reshape((output_dim, input_dim), order="F")
-                    kraus.append(k)
+        values, vecs = scipy.linalg.eigh(data)
+        # If we're not a CP map, fall-through back to the general handling.  Since we needed to get
+        # the eigenvalues anyway, we can do the CP check manually rather than deferring to a
+        # separate re-calculation.
+        if all(values >= -atol):
+            kraus = [
+                math.sqrt(value) * vec.reshape((output_dim, input_dim), order="F")
+                for value, vec in zip(values, vecs.T)
+                if abs(value) > atol
+            ]
             # If we are converting a zero matrix, we need to return a Kraus set
             # with a single zero-element Kraus matrix
             if not kraus:
-                kraus.append(np.zeros((output_dim, input_dim), dtype=complex))
+                kraus = [np.zeros((output_dim, input_dim), dtype=complex)]
             return kraus, None
-    # Non-CP-map generalized Kraus representation
-    mat_u, svals, mat_vh = la.svd(data)
+        # Fall through to the general non-CP form.
+    # Non-CP-map generalized Kraus representation.
+    mat_u, svals, mat_vh = scipy.linalg.svd(data)
     kraus_l = []
     kraus_r = []
     for val, vec_l, vec_r in zip(svals, mat_u.T, mat_vh.conj()):

qiskit/quantum_info/operators/operator.py

@@ -16,6 +16,7 @@
 
 from __future__ import annotations
 
+import cmath
 import copy as _copy
 import re
 from numbers import Number
@@ -561,13 +562,19 @@ def power(self, n: float) -> Operator:
         else:
             import scipy.linalg
 
-            # Experimentally, for fractional powers this seems to be 3x faster than
-            # calling scipy.linalg.fractional_matrix_power(self.data, n)
-            decomposition, unitary = scipy.linalg.schur(self.data, output="complex")
-            decomposition_diagonal = decomposition.diagonal()
-            decomposition_power = [pow(element, n) for element in decomposition_diagonal]
-            unitary_power = unitary @ np.diag(decomposition_power) @ unitary.conj().T
-            ret._data = unitary_power
+            # Non-integer powers of operators with an eigenvalue of -1 have a branch cut in the
+            # complex plane, which makes the calculation of the principal root around this cut
+            # finnicky and subject to precision / differences in BLAS implementation.  For example,
+            # the square root of Pauli Y can return the pi/2 or -pi/2 Y rotation depending on
+            # whether the -1 eigenvalue is found as `complex(-1, tiny)` or `complex(-1, -tiny)`.
+            # Such -1 eigenvalues are really common, so we first phase the matrix slightly to move
+            # common eigenvalues away from the branch-cut point of the power function.  The exact
+            # value of the epsilon doesn't matter much, but shouldn't coincide with any "nice"
+            # eigenvalues we expect to encounter. This isn't perfect, just pragmatic.
+            eps = cmath.pi * 1e-3
+            ret._data = cmath.rect(1, eps * n) * scipy.linalg.fractional_matrix_power(
+                cmath.rect(1, -eps) * self.data, n
+            )
         return ret
 
     def tensor(self, other: Operator) -> Operator:

releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml

@@ -0,0 +1,18 @@
+---
+fixes:
+  - |
+    Fixed :meth:`.Operator.power` when called with non-integer powers on a matrix whose Schur form
+    is not diagonal (for example, most non-unitary matrices).
+  - |
+    :meth:`.Operator.power` will now more reliably return the expected principal value from a
+    fractional matrix power of a unitary matrix with a :math:`-1` eigenvalue.  This is tricky in
+    general, because floating-point rounding effects can cause a matrix to _truly_ have an eigenvalue
+    on the negative side of the branch cut (even if its exact mathematical relation would not), and
+    imprecision in various BLAS calls can falsely find the wrong side of the branch cut.
+
+    :meth:`.Operator.power` now shifts the branch-cut location for matrix powers to be a small
+    complex rotation away from :math:`-1`.  This does not solve the problem, it just shifts it to a
+    place where it is far less likely to be noticeable for the types of operators that usually
+    appear.
+
+    See `#13305 <https://github.com/Qiskit/qiskit/issues/13305>`__.

test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py

@@ -271,7 +271,7 @@ def test_to_matrix_zero(self):
 
         zero_sparse = zero.to_matrix(sparse=True)
         self.assertIsInstance(zero_sparse, scipy.sparse.csr_matrix)
-        np.testing.assert_array_equal(zero_sparse.A, zero_numpy)
+        np.testing.assert_array_equal(zero_sparse.todense(), zero_numpy)
 
     def test_to_matrix_parallel_vs_serial(self):
         """Parallel execution should produce the same results as serial execution up to

test/python/quantum_info/operators/test_operator.py

@@ -20,9 +20,10 @@
 
 from test import combine
 import numpy as np
-from ddt import ddt
+import ddt
 from numpy.testing import assert_allclose
-import scipy.linalg as la
+import scipy.stats
+import scipy.linalg
 
 from qiskit import QiskitError
 from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
@@ -97,7 +98,7 @@ def simple_circuit_with_measure(self):
         return circ
 
 
-@ddt
+@ddt.ddt
 class TestOperator(OperatorTestCase):
     """Tests for Operator linear operator class."""
 
@@ -290,7 +291,7 @@ def test_copy(self):
     def test_is_unitary(self):
         """Test is_unitary method."""
         # X-90 rotation
-        X90 = la.expm(-1j * 0.5 * np.pi * np.array([[0, 1], [1, 0]]) / 2)
+        X90 = scipy.linalg.expm(-1j * 0.5 * np.pi * np.array([[0, 1], [1, 0]]) / 2)
         self.assertTrue(Operator(X90).is_unitary())
         # Non-unitary should return false
         self.assertFalse(Operator([[1, 0], [0, 0]]).is_unitary())
@@ -495,7 +496,7 @@ def test_compose_front_subsystem(self):
 
     def test_power(self):
         """Test power method."""
-        X90 = la.expm(-1j * 0.5 * np.pi * np.array([[0, 1], [1, 0]]) / 2)
+        X90 = scipy.linalg.expm(-1j * 0.5 * np.pi * np.array([[0, 1], [1, 0]]) / 2)
         op = Operator(X90)
         self.assertEqual(op.power(2), Operator([[0, -1j], [-1j, 0]]))
         self.assertEqual(op.power(4), Operator(-1 * np.eye(2)))
@@ -513,6 +514,35 @@ def test_floating_point_power(self):
 
         self.assertEqual(op.power(0.25), expected_op)
 
+    def test_power_of_nonunitary(self):
+        """Test power method for a non-unitary matrix."""
+        data = [[1, 1], [0, -1]]
+        powered = Operator(data).power(0.5)
+        expected = Operator([[1.0 + 0.0j, 0.5 - 0.5j], [0.0 + 0.0j, 0.0 + 1.0j]])
+        assert_allclose(powered.data, expected.data)
+
+    @ddt.data(0.5, 1.0 / 3.0, 0.25)
+    def test_root_stability(self, root):
+        """Test that the root of operators that have eigenvalues that are -1 up to floating-point
+        imprecision stably choose the positive side of the principal-root branch cut."""
+        rng = np.random.default_rng(2024_10_22)
+
+        eigenvalues = np.array([1.0, -1.0], dtype=complex)
+        # We have the eigenvalues exactly, so this will safely find the principal root.
+        root_eigenvalues = eigenvalues**root
+
+        # Now, we can construct a bunch of Haar-random unitaries with our chosen eigenvalues.  Since
+        # we already know their eigenvalue decompositions exactly (up to floating-point precision in
+        # the matrix multiplications), we can also compute the principal values of the roots of the
+        # complete matrices.
+        bases = scipy.stats.unitary_group.rvs(2, size=16, random_state=rng)
+        matrices = [basis.conj().T @ np.diag(eigenvalues) @ basis for basis in bases]
+        expected = [basis.conj().T @ np.diag(root_eigenvalues) @ basis for basis in bases]
+        self.assertEqual(
+            [Operator(matrix).power(root) for matrix in matrices],
+            [Operator(single) for single in expected],
+        )
+
     def test_expand(self):
         """Test expand method."""
         mat1 = self.UX