Reorder Pauli terms before Trotterization

qiskit/circuit/library/pauli_evolution.py

@@ -110,17 +110,19 @@ class docstring for an example.
         else:
             operator = _to_sparse_pauli_op(operator)
 
-        if synthesis is None:
-            from qiskit.synthesis.evolution import LieTrotter
-
-            synthesis = LieTrotter()
-
         if label is None:
             label = _get_default_label(operator)
 
         num_qubits = operator[0].num_qubits if isinstance(operator, list) else operator.num_qubits
         super().__init__(name="PauliEvolution", num_qubits=num_qubits, params=[time], label=label)
         self.operator = operator
+
+        if synthesis is None:
+            # pylint: disable=cyclic-import
+            from qiskit.synthesis.evolution import LieTrotter
+
+            synthesis = LieTrotter()
+
         self.synthesis = synthesis
 
     @property

qiskit/synthesis/evolution/lie_trotter.py

@@ -60,6 +60,7 @@ def __init__(
             | None
         ) = None,
         wrap: bool = False,
+        preserve_order: bool = False,
     ) -> None:
         """
         Args:
@@ -79,8 +80,19 @@ def __init__(
                 built.
             wrap: Whether to wrap the atomic evolutions into custom gate objects. This only takes
                 effect when ``atomic_evolution is None``.
+            preserve_order: If ``False``, allows 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,
+            preserve_order=preserve_order,
+        )
 
     @property
     def settings(self) -> dict[str, Any]:

qiskit/synthesis/evolution/product_formula.py

@@ -15,9 +15,13 @@
 from __future__ import annotations
 
 import inspect
-from collections.abc import Callable
+import itertools
+from collections.abc import Callable, Sequence
+from collections import defaultdict
+from itertools import combinations
 import typing
 import numpy as np
+import rustworkx as rx
 from qiskit.circuit.parameterexpression import ParameterExpression
 from qiskit.circuit.quantumcircuit import QuantumCircuit, ParameterValueType
 from qiskit.quantum_info import SparsePauliOp, Pauli
@@ -29,6 +33,8 @@
 if typing.TYPE_CHECKING:
     from qiskit.circuit.library import PauliEvolutionGate
 
+SparsePauliLabel = typing.Tuple[str, list[int], ParameterValueType]
+
 
 class ProductFormula(EvolutionSynthesis):
     """Product formula base class for the decomposition of non-commuting operator exponentials.
@@ -63,6 +69,7 @@ def __init__(
             | None
         ) = None,
         wrap: bool = False,
+        preserve_order: bool = False,
     ) -> None:
         """
         Args:
@@ -84,11 +91,15 @@ def __init__(
             wrap: Whether to wrap the atomic evolutions into custom gate objects. Note that setting
                 this to ``True`` is slower than ``False``. This only takes effect when
                 ``atomic_evolution is None``.
+            preserve_order: If ``False``, allows 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.preserve_order = preserve_order
 
         # user-provided atomic evolution, stored for serialization
         self._atomic_evolution = atomic_evolution
@@ -177,6 +188,7 @@ def settings(self) -> dict[str, typing.Any]:
             "insert_barriers": self.insert_barriers,
             "cx_structure": self._cx_structure,
             "wrap": self._wrap,
+            "preserve_order": self.preserve_order,
         }
 
     def _normalize_coefficients(
@@ -239,3 +251,61 @@ def real_or_fail(value, tol=100):
         return np.real(value)
 
     raise ValueError(f"Encountered complex value {value}, but expected real.")
+
+
+def reorder_paulis(
+    paulis: Sequence[SparsePauliLabel],
+    strategy: rx.ColoringStrategy = rx.ColoringStrategy.Saturation,
+) -> list[SparsePauliLabel]:
+    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 terms act on the same qubit (for example, the terms :math:`IXX` and
+    :math:`IYI` would be connected, but not :math:`IXX` and :math:`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 ``paulis``.
+
+    This method is deterministic and invariant under permutation of the Pauli
+    term in ``paulis``.
+
+    Args:
+        paulis: The operator whose terms to reorder.
+        strategy: The coloring heuristic to use, see ``ColoringStrategy`` [#].
+            Default is ``ColoringStrategy.Saturation``.
+
+    .. [#] https://www.rustworkx.org/apiref/rustworkx.ColoringStrategy.html#coloringstrategy
+
+    """
+
+    def _term_sort_key(term: SparsePauliLabel) -> typing.Any:
+        # sort by index, then by pauli
+        return (term[1], term[0])
+
+    # Do nothing in trivial cases
+    if len(paulis) <= 1:
+        return paulis
+
+    terms = sorted(paulis, key=_term_sort_key)
+    graph = rx.PyGraph()
+    graph.add_nodes_from(terms)
+    indexed_nodes = list(enumerate(graph.nodes()))
+    for (idx1, (_, ind1, _)), (idx2, (_, ind2, _)) in combinations(indexed_nodes, 2):
+        # Add an edge between two terms if they touch the same qubit
+        if len(set(ind1).intersection(ind2)) > 0:
+            graph.add_edge(idx1, idx2, None)
+
+    # rx.graph_greedy_color is supposed to be deterministic
+    coloring = rx.graph_greedy_color(graph, strategy=strategy)
+    terms_by_color = defaultdict(list)
+
+    for term_idx, color in sorted(coloring.items()):
+        term = graph.nodes()[term_idx]
+        terms_by_color[color].append(term)
+
+    terms = list(itertools.chain(*terms_by_color.values()))
+    return terms

qiskit/synthesis/evolution/qdrift.py

@@ -25,7 +25,7 @@
 from qiskit.utils.deprecation import deprecate_arg
 from qiskit.exceptions import QiskitError
 
-from .product_formula import ProductFormula
+from .product_formula import ProductFormula, reorder_paulis
 
 if typing.TYPE_CHECKING:
     from qiskit.circuit.library import PauliEvolutionGate
@@ -68,6 +68,7 @@ def __init__(
         ) = None,
         seed: int | None = None,
         wrap: bool = False,
+        preserve_order: bool = True,
     ) -> None:
         r"""
         Args:
@@ -88,8 +89,13 @@ def __init__(
             seed: An optional seed for reproducibility of the random sampling process.
             wrap: Whether to wrap the atomic evolutions into custom gate objects. This only takes
                 effect when ``atomic_evolution is None``.
+            preserve_order: If ``False``, allows 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, preserve_order
+        )
         self.sampled_ops = None
         self.rng = np.random.default_rng(seed)
 
@@ -125,4 +131,8 @@ def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], f
         sampled_paulis = [
             (pauli[0], pauli[1], np.real(np.sign(pauli[2])) * rescaled_time) for pauli in sampled
         ]
+
+        if not self.preserve_order:
+            sampled_paulis = reorder_paulis(sampled_paulis)
+
         return sampled_paulis

qiskit/synthesis/evolution/suzuki_trotter.py

@@ -23,7 +23,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
 
 if typing.TYPE_CHECKING:
     from qiskit.circuit.quantumcircuit import ParameterValueType
@@ -85,6 +85,7 @@ def __init__(
             | None
         ) = None,
         wrap: bool = False,
+        preserve_order: bool = False,
     ) -> None:
         """
         Args:
@@ -104,6 +105,9 @@ def __init__(
                 built.
             wrap: Whether to wrap the atomic evolutions into custom gate objects. This only takes
                 effect when ``atomic_evolution is None``.
+            preserve_order: If ``False``, allows 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
         """
@@ -113,7 +117,15 @@ def __init__(
                 "Suzuki product formulae are symmetric and therefore only defined "
                 f"for when the order is 1 or even, not {order}."
             )
-        super().__init__(order, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
+        super().__init__(
+            order,
+            reps,
+            insert_barriers,
+            cx_structure,
+            atomic_evolution,
+            wrap,
+            preserve_order=preserve_order,
+        )
 
     def expand(
         self, evolution: PauliEvolutionGate
@@ -144,15 +156,20 @@ def expand(
         operators = evolution.operator  # type: SparsePauliOp | list[SparsePauliOp]
         time = evolution.time
 
+        def to_sparse_list(operator):
+            paulis = (time * (2 / self.reps) * operator).to_sparse_list()
+            if not self.preserve_order:
+                return reorder_paulis(paulis)
+
+            return paulis
+
         # construct the evolution circuit
         if isinstance(operators, list):  # already sorted into commuting bits
-            non_commuting = [
-                (2 / self.reps * time * operator).to_sparse_list() for operator in operators
-            ]
+            non_commuting = [to_sparse_list(operator) for operator in operators]
         else:
             # Assume no commutativity here. If we were to group commuting Paulis,
             # here would be the location to do so.
-            non_commuting = [[op] for op in (2 / self.reps * time * operators).to_sparse_list()]
+            non_commuting = [[op] for op in to_sparse_list(operators)]
 
         # normalize coefficients, i.e. ensure they are float or ParameterExpression
         non_commuting = self._normalize_coefficients(non_commuting)

qiskit/transpiler/passes/synthesis/hls_plugins.py

@@ -1463,11 +1463,14 @@ def run(self, high_level_object, coupling_map=None, target=None, qubits=None, **
 class PauliEvolutionSynthesisDefault(HighLevelSynthesisPlugin):
     """Synthesize a :class:`.PauliEvolutionGate` using the default synthesis algorithm.
 
-    This plugin name is :``PauliEvolution.default`` which can be used as the key on
+    This plugin name is:``PauliEvolution.default`` which can be used as the key on
     an :class:`~.HLSConfig` object to use this method with :class:`~.HighLevelSynthesis`.
 
-    The default synthesis simply calls the synthesis algorithm attached to a
-    PauliEvolutionGate.
+    The following plugin option can be set:
+
+    * preserve_order: If ``False``, allow re-ordering the Pauli terms in the Hamiltonian to
+        reduce the circuit depth of the decomposition.
+
     """
 
     def run(self, high_level_object, coupling_map=None, target=None, qubits=None, **options):
@@ -1477,6 +1480,10 @@ def run(self, high_level_object, coupling_map=None, target=None, qubits=None, **
             return None
 
         algo = high_level_object.synthesis
+
+        if "preserve_order" in options and isinstance(algo, ProductFormula):
+            algo.preserve_order = options["preserve_order"]
+
         return algo.synthesize(high_level_object)
 
 
@@ -1528,6 +1535,9 @@ def run(self, high_level_object, coupling_map=None, target=None, qubits=None, **
             )
             return None
 
+        if "preserve_order" in options and isinstance(algo, ProductFormula):
+            algo.preserve_order = options["preserve_order"]
+
         num_qubits = high_level_object.num_qubits
         pauli_network = algo.expand(high_level_object)
 

releasenotes/notes/reorder-trotter-terms-c8a6eb3cdb831f77.yaml

@@ -0,0 +1,40 @@
+---
+features_synthesis:
+  - |
+    Added a new argument ``preserve_order`` to :class:`.ProductFormula`, which allows
+    re-ordering the Pauli terms in the Hamiltonian before the product formula expansion,
+    to compress the final circuit depth. By setting this to ``False``, a term of form
+
+    .. math::
+
+      Z_0 Z_1 + X_1 X_2 + Y_2 Y_3
+
+    will be re-ordered to
+
+    .. math::
+
+      Z_0 Z_1 + Y_2 Y_3 + X_1 X_2
+
+    which will lead to the ``RZZ`` and ``RYY`` rotations being applied in parallel, instead
+    of three sequential rotations in the first part.
+
+    This option can be set via the plugin interface::
+
+        from qiskit import QuantumCircuit, transpile
+        from qiskit.circuit.library import PauliEvolutionGate
+        from qiskit.quantum_info import SparsePauliOp
+        from qiskit.synthesis.evolution import SuzukiTrotter
+        from qiskit.transpiler.passes import HLSConfig
+
+        op = SparsePauliOp(["XXII", "IYYI", "IIZZ"])
+        time, reps = 0.1, 1
+
+        synthesis = SuzukiTrotter(order=2, reps=reps)
+        hls_config = HLSConfig(PauliEvolution=[("default", {"preserve_order": False})])
+
+        circuit = QuantumCircuit(op.num_qubits)
+        circuit.append(PauliEvolutionGate(op, time), circuit.qubits)
+
+        tqc = transpile(circuit, basis_gates=["u", "cx"], hls_config=hls_config)
+        print(tqc.draw())
+