From 847064880a600dad923374ed64ba0b293e3fa332 Mon Sep 17 00:00:00 2001 From: Jake Lishman Date: Tue, 22 Oct 2024 18:06:46 +0100 Subject: [PATCH 1/2] 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 test changes to the Kraus and Stinespring modules are to cope with the two operators only being defined up to a global phase, which the test previously did not account for. The conversion to Kraus-operator form happens to work fine with OpenBLAS, but caused global-phase differences on macOS Accelerate. A previous version of this commit attempted to revert the Choi-to-Kraus conversion back to using `eigh` instead of the Schur decomposition, but the `eigh` instabilities noted in fdd5603af76 (gh-3884) (the time of Scipy 1.1 to 1.3, with OpenBLASes around 0.3.6) remain with Scipy 1.13/1.14 and OpenBLAS 0.3.27. Co-authored-by: Lev S. Bishop <18673315+levbishop@users.noreply.github.com> Co-authored-by: Alexander Ivrii --- constraints.txt | 4 -- .../operators/channel/transformations.py | 49 +++++++++++-------- qiskit/quantum_info/operators/operator.py | 21 +++++--- .../notes/scipy-1.14-951d1c245473aee9.yaml | 18 +++++++ .../operators/channel/test_kraus.py | 11 ++++- .../operators/channel/test_stinespring.py | 7 ++- .../symplectic/test_sparse_pauli_op.py | 2 +- .../quantum_info/operators/test_operator.py | 40 +++++++++++++-- 8 files changed, 110 insertions(+), 42 deletions(-) create mode 100644 releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml 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..79e665e1b38c 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 @@ -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: diff --git a/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml b/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml new file mode 100644 index 000000000000..939b8fe8dee6 --- /dev/null +++ b/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 `__. 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..f206da551169 100644 --- a/test/python/quantum_info/operators/test_operator.py +++ b/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 From 38ddaf8e553b43ac9268548c32a7c9b9c117a2a5 Mon Sep 17 00:00:00 2001 From: Jake Lishman Date: Thu, 24 Oct 2024 16:00:11 +0100 Subject: [PATCH 2/2] Expose `branch_cut_rotation` parameter in `Operator.power` The rotation used to stabilise matrix roots has an impact on which matrix is selected as the principal root. Exposing it to users to allow control makes the most sense. --- qiskit/quantum_info/operators/operator.py | 44 +++++++++++++------ .../notes/scipy-1.14-951d1c245473aee9.yaml | 9 +++- .../quantum_info/operators/test_operator.py | 9 ++++ 3 files changed, 48 insertions(+), 14 deletions(-) diff --git a/qiskit/quantum_info/operators/operator.py b/qiskit/quantum_info/operators/operator.py index 79e665e1b38c..fc3dd9ee187a 100644 --- a/qiskit/quantum_info/operators/operator.py +++ b/qiskit/quantum_info/operators/operator.py @@ -541,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``. @@ -562,18 +588,10 @@ def power(self, n: float) -> Operator: else: import scipy.linalg - # 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 + 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 diff --git a/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml b/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml index 939b8fe8dee6..41f4b8790286 100644 --- a/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml +++ b/releasenotes/notes/scipy-1.14-951d1c245473aee9.yaml @@ -13,6 +13,13 @@ fixes: :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. + 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/test_operator.py b/test/python/quantum_info/operators/test_operator.py index f206da551169..d9423d0ec141 100644 --- a/test/python/quantum_info/operators/test_operator.py +++ b/test/python/quantum_info/operators/test_operator.py @@ -27,6 +27,7 @@ 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 @@ -543,6 +544,14 @@ def test_root_stability(self, root): [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