Reorder Pauli terms before Trotterization

qiskit/circuit/library/pauli_evolution.py

@@ -110,18 +110,24 @@ 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
-        self.synthesis = synthesis
+        self._synthesis = synthesis
+
+    @property
+    def synthesis(self) -> EvolutionSynthesis:
+        """Return the synthesis used."""
+        if self._synthesis is not None:
+            return self._synthesis
+
+        # pylint: disable=cyclic-import
+        from qiskit.synthesis.evolution import LieTrotter
+
+        return LieTrotter()
 
     @property
     def time(self) -> ParameterValueType:
@@ -143,7 +149,14 @@ def time(self, time: ParameterValueType) -> None:
 
     def _define(self):
         """Unroll, where the default synthesis is matrix based."""
-        self.definition = self.synthesis.synthesize(self)
+        if self.synthesis is None:
+            from qiskit.synthesis.evolution import LieTrotter
+
+            synthesis = LieTrotter()
+        else:
+            synthesis = self.synthesis
+
+        self.definition = synthesis.synthesize(self)
 
     def validate_parameter(self, parameter: ParameterValueType) -> ParameterValueType:
         """Gate parameters should be int, float, or ParameterExpression"""

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: Whether to allow preserve_ordering the terms of the operator to
+                potentially yield a shallower evolution circuit. Not relevant
+                when synthesizing operator with a single term. Defaults to ``True``.
         """
-        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

@@ -16,8 +16,11 @@
 
 import inspect
 from collections.abc import Callable
+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 +32,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 +68,7 @@ def __init__(
             | None
         ) = None,
         wrap: bool = False,
+        preserve_order: bool = False,
     ) -> None:
         """
         Args:
@@ -84,11 +90,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: Whether to allow preserve_ordering 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 +187,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 +250,74 @@ 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: SparsePauliLabel, strategy: rx.ColoringStrategy = rx.ColoringStrategy.Saturation
+) -> SparsePauliOp | list[SparsePauliOp]:
+    r"""
+    Creates an equivalent operator by preserve_ordering 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 :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 ``operator``. Finally, a new
+    :class:`~qiskit.quantum_info.SparsePauliOp` is created where terms of the
+    same color are grouped together.
+
+    In trivial cases, i.e. when either
+    - the input is a :class:`~qiskit.quantum_info.SparsePauliOp` with
+      less than two Pauli terms, or
+    - the input is a list containing a single
+      :class:`~qiskit.quantum_info.SparsePauliOp` which has less than two
+      Pauli terms,
+    this method does nothing.
+
+    If ``operators`` is a list of :class:`~qiskit.quantum_info.SparsePauliOp`,
+    then preserve_ordering is applied to every operator independently, and the list of
+    preserve_ordered operators is returned.
+
+    This method is deterministic and invariant under permutation of the Pauli
+    term in ``operators``.
+
+    Args:
+        paulis: The operator whose terms to preserve_order.
+        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 = sum(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,14 @@ 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: Whether to allow reordering the terms of the sampled
+                operator before synthesizing the evolution circuit. Setting this
+                to ``False`` can potentially yield a shallower evolution circuit.
+                Defaults to ``True``.
         """
-        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 +132,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: 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. Defaults to ``True``.
         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

@@ -282,6 +282,7 @@
     synth_clifford_ag,
     synth_clifford_bm,
 )
+from qiskit.synthesis.evolution import ProductFormula
 from qiskit.synthesis.linear import (
     synth_cnot_count_full_pmh,
     synth_cnot_depth_line_kms,
@@ -1047,8 +1048,22 @@ def run(self, high_level_object, coupling_map=None, target=None, qubits=None, **
 
 
 class PauliEvolutionSynthesisDefault(HighLevelSynthesisPlugin):
-    """The default implementation calling the attached synthesis algorithm."""
+    """The default implementation calling the attached synthesis algorithm.
+
+    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 following plugin option can be set:
+
+    * preserve_order: If ``True``, 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):
         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)