Reorder Pauli terms before Trotterization

qiskit/synthesis/evolution/lie_trotter.py

@@ -22,7 +22,7 @@
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
 from qiskit.utils.deprecation import deprecate_arg
 
-from .product_formula import ProductFormula
+from .product_formula import ProductFormula, reorder_paulis
 
 
 class LieTrotter(ProductFormula):
@@ -78,6 +78,7 @@ def __init__(
             | None
         ) = None,
         wrap: bool = False,
+        reorder: bool = False,
     ) -> None:
         """
         Args:
@@ -97,12 +98,25 @@ def __init__(
                 built.
             wrap: Whether to wrap the atomic evolutions into custom gate objects. This only takes
                 effect when ``atomic_evolution is None``.
+            reorder: Whether to allow reordering the terms of the operator to
+                potentially yield a shallower evolution circuit. Not relevant
+                when synthesizing operator with a single term.
         """
-        super().__init__(1, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
+        super().__init__(
+            1,
+            reps,
+            insert_barriers,
+            cx_structure,
+            atomic_evolution,
+            wrap,
+            reorder=reorder,
+        )
 
     def synthesize(self, evolution):
         # get operators and time to evolve
         operators = evolution.operator
+        if self.reorder:
+            operators = reorder_paulis(operators)
         time = evolution.time
 
         # construct the evolution circuit

qiskit/synthesis/evolution/product_formula.py

@@ -16,9 +16,12 @@
 
 import inspect
 from collections.abc import Callable
+from collections import defaultdict
+from itertools import combinations
 from typing import Any
 from functools import partial
 import numpy as np
+import rustworkx as rx
 from qiskit.circuit.parameterexpression import ParameterExpression
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.quantum_info import SparsePauliOp, Pauli
@@ -60,6 +63,7 @@ def __init__(
             | None
         ) = None,
         wrap: bool = False,
+        reorder: bool = False,
     ) -> None:
         """
         Args:
@@ -80,11 +84,15 @@ def __init__(
                 built.
             wrap: Whether to wrap the atomic evolutions into custom gate objects. This only takes
                 effect when ``atomic_evolution is None``.
+            reorder: Whether to allow reordering the terms of the operator to
+                potentially yield a shallower evolution circuit. Not relevant
+                when synthesizing operator with a single term.
         """
         super().__init__()
         self.order = order
         self.reps = reps
         self.insert_barriers = insert_barriers
+        self.reorder = reorder
 
         # user-provided atomic evolution, stored for serialization
         self._atomic_evolution = atomic_evolution
@@ -129,6 +137,7 @@ def settings(self) -> dict[str, Any]:
             "insert_barriers": self.insert_barriers,
             "cx_structure": self._cx_structure,
             "wrap": self._wrap,
+            "reorder": self.reorder,
         }
 
 
@@ -177,6 +186,64 @@ def evolve_pauli(
         _multi_qubit_evolution(output, pauli, time, cx_structure, wrap)
 
 
+def reorder_paulis(
+    operators: SparsePauliOp | list[SparsePauliOp],
+) -> SparsePauliOp | list[SparsePauliOp]:
+    r"""
+    Creates an equivalent operator by reordering terms in order to yield a
+    shallower circuit after evolution synthesis. The original operator remains
+    unchanged.
+
+    This method works in three steps. First, a graph is constructed, where the
+    nodes are the terms of the operator and where two nodes are connected if
+    their term acts on the same qubit (for example, the terms $IXX$ and $IYI$
+    would be connected, but not $IXX$ and $YII$). Then, the graph is colored.
+    Two terms with the same color thus do not act on the same qubit, and in
+    particular, their evolution subcircuits can be run in parallel in the
+    greater evolution circuit of ``operator``. Finally, a new
+    :class:`~qiskit.quantum_info.SparsePauliOp` is created where terms of the
+    same color are grouped together.
+
+    If the input is in fact a list of
+    :class:`~qiskit.quantum_info.SparsePauliOp`, then the terms of all operators
+    will be coalesced and reordered into a single
+    :class:`~qiskit.quantum_info.SparsePauliOp`.
+
+    Args:
+        operators: The operator or list of operators whose terms to reorder.
+    """
+    if not isinstance(operators, list):
+        operators = [operators]
+    # Do nothing in trivial cases
+    if not (
+        # operators is a list of > 1 SparsePauliOp
+        len(operators) > 1
+        # operators has a single SparsePauliOp, which has > 1 terms
+        or (len(operators) == 1 and len(operators[0]) > 1)
+    ):
+        return operators
+    graph = rx.PyGraph()
+    graph.add_nodes_from(sum((o.to_list() for o in operators), []))
+    indexed_nodes = list(enumerate(graph.nodes()))
+    for (idx1, term1), (idx2, term2) in combinations(indexed_nodes, 2):
+        # Add an edge between term1 and term2 if they touch the same qubit
+        for a, b in zip(term1[0], term2[0]):
+            if not (a == "I" or b == "I"):
+                graph.add_edge(idx1, idx2, None)
+                break
+    # TODO: The graph is likely not very big, so an exact coloring could be
+    # computed in a reasonable time.
+    coloring = rx.graph_greedy_color(graph)
+    indices_by_color = defaultdict(list)
+    for term, color in coloring.items():
+        indices_by_color[color].append(term)
+    terms: list[tuple[str, complex]] = []
+    for color, indices in indices_by_color.items():
+        for i in indices:
+            terms.append(graph.nodes()[i])
+    return SparsePauliOp.from_list(terms)
+
+
 def _single_qubit_evolution(output, pauli, time, wrap):
     dest = QuantumCircuit(1) if wrap else output
     # Note that all phases are removed from the pauli label and are only in the coefficients.

qiskit/synthesis/evolution/suzuki_trotter.py

@@ -24,7 +24,7 @@
 from qiskit.utils.deprecation import deprecate_arg
 
 
-from .product_formula import ProductFormula
+from .product_formula import ProductFormula, reorder_paulis
 
 
 class SuzukiTrotter(ProductFormula):
@@ -82,6 +82,7 @@ def __init__(
             | None
         ) = None,
         wrap: bool = False,
+        reorder: bool = False,
     ) -> None:
         """
         Args:
@@ -101,6 +102,9 @@ def __init__(
                 built.
             wrap: Whether to wrap the atomic evolutions into custom gate objects. This only takes
                 effect when ``atomic_evolution is None``.
+            reorder: Whether to allow reordering the terms of the operator to
+                potentially yield a shallower evolution circuit. Not relevant
+                when synthesizing operator with a single term.
         Raises:
             ValueError: If order is not even
         """
@@ -110,11 +114,21 @@ def __init__(
                 "Suzuki product formulae are symmetric and therefore only defined "
                 "for even orders."
             )
-        super().__init__(order, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
+        super().__init__(
+            order,
+            reps,
+            insert_barriers,
+            cx_structure,
+            atomic_evolution,
+            wrap,
+            reorder=reorder,
+        )
 
     def synthesize(self, evolution):
         # get operators and time to evolve
         operators = evolution.operator
+        if self.reorder:
+            operators = reorder_paulis(operators)
         time = evolution.time
 
         if not isinstance(operators, list):
@@ -150,6 +164,8 @@ def _recurse(order, time, pauli_list):
                 order - 2, time=reduction * time, pauli_list=pauli_list
             )
             inner = SuzukiTrotter._recurse(
-                order - 2, time=(1 - 4 * reduction) * time, pauli_list=pauli_list
+                order - 2,
+                time=(1 - 4 * reduction) * time,
+                pauli_list=pauli_list,
             )
             return outer + inner + outer

test/python/circuit/library/test_evolution_gate.py

@@ -19,8 +19,16 @@
 
 from qiskit.circuit import QuantumCircuit, Parameter
 from qiskit.circuit.library import PauliEvolutionGate
-from qiskit.synthesis import LieTrotter, SuzukiTrotter, MatrixExponential, QDrift
-from qiskit.synthesis.evolution.product_formula import cnot_chain, diagonalizing_clifford
+from qiskit.synthesis import (
+    LieTrotter,
+    SuzukiTrotter,
+    MatrixExponential,
+    QDrift,
+)
+from qiskit.synthesis.evolution.product_formula import (
+    cnot_chain,
+    diagonalizing_clifford,
+)
 from qiskit.converters import circuit_to_dag
 from qiskit.quantum_info import Operator, SparsePauliOp, Pauli, Statevector
 from test import QiskitTestCase  # pylint: disable=wrong-import-order
@@ -133,6 +141,23 @@ def test_suzuki_trotter_manual(self):
 
         self.assertEqual(evo_gate.definition, expected)
 
+    def test_suzuki_trotter_reordered_manual(self):
+        """Test the evolution circuit of Suzuki Trotter against a manually constructed circuit."""
+        op = (X ^ X ^ I ^ I) + (I ^ Y ^ Y ^ I) + (I ^ I ^ Z ^ Z)
+        time, reps = 0.1, 1
+        evo_gate = PauliEvolutionGate(
+            op,
+            time,
+            synthesis=SuzukiTrotter(order=2, reps=reps, reorder=True),
+        )
+        expected = QuantumCircuit(4)
+        expected.ryy(time, 1, 2)
+        expected.rxx(time, 2, 3)
+        expected.rzz(2 * time, 0, 1)
+        expected.rxx(time, 2, 3)
+        expected.ryy(time, 1, 2)
+        self.assertEqual(evo_gate.definition, expected)
+
     @data(
         (X + Y, 0.5, 1, [(Pauli("X"), 0.5), (Pauli("X"), 0.5)]),
         (X, 0.238, 2, [(Pauli("X"), 0.238)]),
@@ -299,6 +324,23 @@ def test_lie_trotter_two_qubit_correct_order(self):
 
         self.assertTrue(Operator(lie_trotter).equiv(exact))
 
+    def test_lie_trotter_reordered_manual(self):
+        """Test the evolution circuit of Lie Trotter against a manually constructed circuit."""
+        op = (X ^ I ^ I ^ I) + (X ^ X ^ I ^ I) + (I ^ Y ^ Y ^ I) + (I ^ I ^ Z ^ Z)
+        time, reps = 0.1, 1
+        evo_gate = PauliEvolutionGate(
+            op,
+            time,
+            synthesis=LieTrotter(reps=reps, reorder=True),
+        )
+        # manually construct expected evolution
+        expected = QuantumCircuit(4)
+        expected.rxx(2 * time, 2, 3)
+        expected.rzz(2 * time, 0, 1)
+        expected.rx(2 * time, 3)
+        expected.ryy(2 * time, 1, 2)
+        self.assertEqual(evo_gate.definition, expected)
+
     def test_complex_op_raises(self):
         """Test an operator with complex coefficient raises an error."""
         with self.assertRaises(ValueError):
@@ -359,7 +401,9 @@ def atomic_evolution(pauli, time):
         reps = 4
         with self.assertWarns(PendingDeprecationWarning):
             evo_gate = PauliEvolutionGate(
-                op, time, synthesis=LieTrotter(reps=reps, atomic_evolution=atomic_evolution)
+                op,
+                time,
+                synthesis=LieTrotter(reps=reps, atomic_evolution=atomic_evolution),
             )
         decomposed = evo_gate.definition.decompose()
         self.assertEqual(decomposed.count_ops()["cx"], reps * 3 * 4)