Fix compatibility issues with SciPy 1.14

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,39 @@ 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)
-            # If we are converting a zero matrix, we need to return a Kraus set
-            # with a single zero-element Kraus matrix
+        # Ideally we'd use `eigh`, but `scipy.linalg.eigh` has stability problems on macOS (at a
+        # minimum from SciPy 1.1 to 1.13 with the bundled OpenBLAS, or ~0.3.6 before they started
+        # bundling one in).  The Schur form of a Hermitian matrix is guaranteed diagonal:
+        #
+        #       H = U T U+   for upper-triangular T.
+        #   => H+ = U T+ U+
+        #   =>  T = T+       because H = H+, and thus T cannot have super-diagonal elements.
+        #
+        # So the eigenvalues are on the diagonal, therefore the basis-transformation matrix must be
+        # a spanning set of the eigenspace.
+        triangular, vecs = scipy.linalg.schur(data)
+        values = triangular.diagonal().real
+        # If we're not a CP map, fall-through back to the generalization 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.
+    # 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
@@ -540,11 +541,37 @@ def compose(self, other: Operator, qargs: list | None = None, front: bool = Fals
         ret._op_shape = new_shape
         return ret
 
-    def power(self, n: float) -> Operator:
+    def power(self, n: float, branch_cut_rotation=cmath.pi * 1e-12) -> Operator:
         """Return the matrix power of the operator.
 
+        Non-integer powers of operators with an eigenvalue whose complex phase is :math:`\\pi` have
+        a branch cut in the complex plane, which makes the calculation of the principal root around
+        this cut subject to precision / differences in BLAS implementation.  For example, the square
+        root of Pauli Y can return the :math:`\\pi/2` or :math:`\\-pi/2` Y rotation depending on
+        whether the -1 eigenvalue is found as ``complex(-1, tiny)`` or ``complex(-1, -tiny)``. Such
+        eigenvalues are really common in quantum information, so this function first phase-rotates
+        the input matrix to shift the branch cut to a far less common point.  The underlying
+        numerical precision issues around the branch-cut point remain, if an operator has an
+        eigenvalue close to this phase.  The magnitude of this rotation can be controlled with the
+        ``branch_cut_rotation`` parameter.
+
+        The choice of ``branch_cut_rotation`` affects the principal root that is found.  For
+        example, the square root of :class:`.ZGate` will be calculated as either :class:`.SGate` or
+        :class:`.SdgGate` depending on which way the rotation is done::
+
+            from qiskit.circuit import library
+            from qiskit.quantum_info import Operator
+
+            z_op = Operator(library.ZGate())
+            assert z_op.power(0.5, branch_cut_rotation=1e-3) == Operator(library.SGate())
+            assert z_op.power(0.5, branch_cut_rotation=-1e-3) == Operator(library.SdgGate())
+
         Args:
             n (float): the power to raise the matrix to.
+            branch_cut_rotation (float): The rotation angle to apply to the branch cut in the
+                complex plane.  This shifts the branch cut away from the common point of :math:`-1`,
+                but can cause a different root to be selected as the principal root.  The rotation
+                is anticlockwise, following the standard convention for complex phase.
 
         Returns:
             Operator: the resulting operator ``O ** n``.
@@ -561,13 +588,11 @@ 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
+            ret._data = cmath.rect(
+                1, branch_cut_rotation * n
+            ) * scipy.linalg.fractional_matrix_power(
+                cmath.rect(1, -branch_cut_rotation) * self.data, n
+            )
         return ret
 
     def tensor(self, other: Operator) -> Operator:

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

@@ -0,0 +1,25 @@
+---
+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.  Use the new ``branch_cut_rotation`` parameter to have more control over this.
+
+    See `#13305 <https://github.com/Qiskit/qiskit/issues/13305>`__.
+features_quantum_info:
+  - |
+    The method :meth:`.Operator.power` has a new parameter ``branch_cut_rotation``.  This can be
+    used to shift the branch-cut point of the root around, which can affect which matrix is chosen
+    as the principal root.  By default, it is set to a small positive rotation to make roots of
+    operators with a real-negative eigenvalue (like Pauli operators) more stable against numerical
+    precision differences.

test/python/quantum_info/operators/channel/test_kraus.py

@@ -19,7 +19,7 @@
 
 from qiskit import QiskitError
 from qiskit.quantum_info.states import DensityMatrix
-from qiskit.quantum_info import Kraus
+from qiskit.quantum_info import Kraus, Operator
 from .channel_test_case import ChannelTestCase
 
 
@@ -68,7 +68,14 @@ def test_circuit_init(self):
         circuit, target = self.simple_circuit_no_measure()
         op = Kraus(circuit)
         target = Kraus(target)
-        self.assertEqual(op, target)
+        # The given separable circuit should only have a single Kraus operator.
+        self.assertEqual(len(op.data), 1)
+        self.assertEqual(len(target.data), 1)
+        kraus_op = Operator(op.data[0])
+        kraus_target = Operator(target.data[0])
+        # THe Kraus representation is not unique, but for a single operator, the only gauge freedom
+        # is the global phase.
+        self.assertTrue(kraus_op.equiv(kraus_target))
 
     def test_circuit_init_except(self):
         """Test initialization from circuit with measure raises exception."""

test/python/quantum_info/operators/channel/test_stinespring.py

@@ -19,7 +19,7 @@
 
 from qiskit import QiskitError
 from qiskit.quantum_info.states import DensityMatrix
-from qiskit.quantum_info import Stinespring
+from qiskit.quantum_info import Stinespring, Operator
 from .channel_test_case import ChannelTestCase
 
 
@@ -61,7 +61,10 @@ def test_circuit_init(self):
         circuit, target = self.simple_circuit_no_measure()
         op = Stinespring(circuit)
         target = Stinespring(target)
-        self.assertEqual(op, target)
+        # If the Stinespring is CPTP (and it should be), it's defined in terms of a single
+        # rectangular operator, which has global-phase gauge freedom.
+        self.assertTrue(op.is_cptp())
+        self.assertTrue(Operator(op.data).equiv(Operator(target.data)))
 
     def test_circuit_init_except(self):
         """Test initialization from circuit with measure raises exception."""

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,12 +20,14 @@
 
 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
+from qiskit.circuit import library
 from qiskit.circuit.library import HGate, CHGate, CXGate, QFT
 from qiskit.transpiler import CouplingMap
 from qiskit.transpiler.layout import Layout, TranspileLayout
@@ -97,7 +99,7 @@ def simple_circuit_with_measure(self):
         return circ
 
 
-@ddt
+@ddt.ddt
 class TestOperator(OperatorTestCase):
     """Tests for Operator linear operator class."""
 
@@ -290,7 +292,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 +497,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 +515,43 @@ 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_root_branch_cut(self):
+        """Test that we can choose where the branch cut appears in the root."""
+        z_op = Operator(library.ZGate())
+        # Depending on the direction we move the branch cut, we should be able to select the root to
+        # be either of the two valid options.
+        self.assertEqual(z_op.power(0.5, branch_cut_rotation=1e-3), Operator(library.SGate()))
+        self.assertEqual(z_op.power(0.5, branch_cut_rotation=-1e-3), Operator(library.SdgGate()))
+
     def test_expand(self):
         """Test expand method."""
         mat1 = self.UX