Add argument assume_unitary to Operator.power. (#13319)

* Add `assume_unitary` kwarg to `Operator.power`

This is approximately a 3x performance improvement over the generalised
`fractional_matrix_power` method, but only works for unitary matrices.

* review comments

qiskit/circuit/gate.py

@@ -89,7 +89,9 @@ def power(self, exponent: float, annotated: bool = False):
         from qiskit.circuit.library.generalized_gates.unitary import UnitaryGate
 
         if not annotated:
-            return UnitaryGate(Operator(self).power(exponent), label=f"{self.name}^{exponent}")
+            return UnitaryGate(
+                Operator(self).power(exponent, assume_unitary=True), label=f"{self.name}^{exponent}"
+            )
         else:
             return AnnotatedOperation(self, PowerModifier(exponent))
 

qiskit/quantum_info/operators/operator.py

@@ -541,7 +541,9 @@ def compose(self, other: Operator, qargs: list | None = None, front: bool = Fals
         ret._op_shape = new_shape
         return ret
 
-    def power(self, n: float, branch_cut_rotation=cmath.pi * 1e-12) -> Operator:
+    def power(
+        self, n: float, branch_cut_rotation=cmath.pi * 1e-12, assume_unitary=False
+    ) -> Operator:
         """Return the matrix power of the operator.
 
         Non-integer powers of operators with an eigenvalue whose complex phase is :math:`\\pi` have
@@ -572,13 +574,20 @@ def power(self, n: float, branch_cut_rotation=cmath.pi * 1e-12) -> Operator:
                 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.
+            assume_unitary (bool): if ``True``, the operator is assumed to be unitary. In this case,
+                for fractional powers we employ a faster implementation based on Schur's decomposition.
 
         Returns:
             Operator: the resulting operator ``O ** n``.
 
         Raises:
             QiskitError: if the input and output dimensions of the operator
                          are not equal.
+
+        .. note::
+            It is only safe to set the argument ``assume_unitary`` to ``True`` when the operator
+            is unitary (or, more generally, normal). Otherwise, the function will return an
+            incorrect output.
         """
         if self.input_dims() != self.output_dims():
             raise QiskitError("Can only power with input_dims = output_dims.")
@@ -588,11 +597,23 @@ def power(self, n: float, branch_cut_rotation=cmath.pi * 1e-12) -> Operator:
         else:
             import scipy.linalg
 
-            ret._data = cmath.rect(
-                1, branch_cut_rotation * n
-            ) * scipy.linalg.fractional_matrix_power(
-                cmath.rect(1, -branch_cut_rotation) * self.data, n
-            )
+            if assume_unitary:
+                # Experimentally, for fractional powers this seems to be 3x faster than
+                # calling scipy.linalg.fractional_matrix_power(self.data, exponent)
+                decomposition, unitary = scipy.linalg.schur(
+                    cmath.rect(1, -branch_cut_rotation) * 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 = cmath.rect(1, branch_cut_rotation * n) * unitary_power
+            else:
+                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/operator-power-assume-unitary-0a2f97ea9de91b49.yaml

@@ -0,0 +1,6 @@
+---
+features_quantum_info:
+  - |
+    Added a new argument ``assume_unitary`` to :meth:`qiskit.quantum_info.Operator.power`.
+    When ``True``, we use a faster method based on Schur's decomposition to raise an
+    ``Operator`` to a fractional power.

test/python/quantum_info/operators/test_operator.py

@@ -522,8 +522,16 @@ def test_power_of_nonunitary(self):
         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):
+    @ddt.data(
+        (0.5, False),
+        (1.0 / 3.0, False),
+        (0.25, False),
+        (0.5, True),
+        (1.0 / 3.0, True),
+        (0.25, True),
+    )
+    @ddt.unpack
+    def test_root_stability(self, root, assume_unitary):
         """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)
@@ -540,17 +548,24 @@ def test_root_stability(self, root):
         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(matrix).power(root, assume_unitary=assume_unitary) for matrix in matrices],
             [Operator(single) for single in expected],
         )
 
-    def test_root_branch_cut(self):
+    @ddt.data(True, False)
+    def test_root_branch_cut(self, assume_unitary):
         """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()))
+        self.assertEqual(
+            z_op.power(0.5, branch_cut_rotation=1e-3, assume_unitary=assume_unitary),
+            Operator(library.SGate()),
+        )
+        self.assertEqual(
+            z_op.power(0.5, branch_cut_rotation=-1e-3, assume_unitary=assume_unitary),
+            Operator(library.SdgGate()),
+        )
 
     def test_expand(self):
         """Test expand method."""