From 38ddaf8e553b43ac9268548c32a7c9b9c117a2a5 Mon Sep 17 00:00:00 2001 From: Jake Lishman Date: Thu, 24 Oct 2024 16:00:11 +0100 Subject: [PATCH] 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