diff --git a/constraints.txt b/constraints.txt index cef0c15114b7..42e11f9e1460 100644 --- a/constraints.txt +++ b/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 diff --git a/qiskit/quantum_info/operators/channel/transformations.py b/qiskit/quantum_info/operators/channel/transformations.py index 8f429cad8cea..657ee62703ec 100644 --- a/qiskit/quantum_info/operators/channel/transformations.py +++ b/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()): diff --git a/qiskit/quantum_info/operators/operator.py b/qiskit/quantum_info/operators/operator.py index 42593626f2cc..fc3dd9ee187a 100644 --- a/qiskit/quantum_info/operators/operator.py +++ b/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: diff --git a/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml b/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml new file mode 100644 index 000000000000..41f4b8790286 --- /dev/null +++ b/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 `__. +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. diff --git a/test/python/quantum_info/operators/channel/test_kraus.py b/test/python/quantum_info/operators/channel/test_kraus.py index 5d50ee9b4759..3b75b2dd614b 100644 --- a/test/python/quantum_info/operators/channel/test_kraus.py +++ b/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.""" diff --git a/test/python/quantum_info/operators/channel/test_stinespring.py b/test/python/quantum_info/operators/channel/test_stinespring.py index 9bcc886a026c..693e85d7c1cc 100644 --- a/test/python/quantum_info/operators/channel/test_stinespring.py +++ b/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.""" diff --git a/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py b/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py index dedd84279a8d..65f19eb8e44c 100644 --- a/test/python/quantum_info/operators/symplectic/test_sparse_pauli_op.py +++ b/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 diff --git a/test/python/quantum_info/operators/test_operator.py b/test/python/quantum_info/operators/test_operator.py index 725c46576a9d..d9423d0ec141 100644 --- a/test/python/quantum_info/operators/test_operator.py +++ b/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