Prepare PauliEvolutionGate for Rustiq & port it to Rust

pyproject.toml

@@ -100,6 +100,7 @@ sk = "qiskit.transpiler.passes.synthesis.solovay_kitaev_synthesis:SolovayKitaevS
 "qft.line" = "qiskit.transpiler.passes.synthesis.hls_plugins:QFTSynthesisLine"
 "qft.default" = "qiskit.transpiler.passes.synthesis.hls_plugins:QFTSynthesisFull"
 "permutation.token_swapper" = "qiskit.transpiler.passes.synthesis.hls_plugins:TokenSwapperSynthesisPermutation"
+"evolution.default" = "qiskit.transpiler.passes.synthesis.hls_plugins:PauliEvolutionSynthesisDefault"
 
 [project.entry-points."qiskit.transpiler.init"]
 default = "qiskit.transpiler.preset_passmanagers.builtin_plugins:DefaultInitPassManager"

qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py

@@ -935,6 +935,14 @@ def to_list(self, array: bool = False):
             return labels
         return labels.tolist()
 
+    def to_sparse_list(self):
+        """Convert to a sparse Pauli list format with elements (pauli, qubits, coefficient)."""
+        pauli_labels = self.paulis.to_labels()
+        sparse_list = [
+            (*sparsify_label(label), coeff) for label, coeff in zip(pauli_labels, self.coeffs)
+        ]
+        return sparse_list
+
     def to_matrix(self, sparse: bool = False, force_serial: bool = False) -> np.ndarray:
         """Convert to a dense or sparse matrix.
 
@@ -1176,5 +1184,13 @@ def apply_layout(
         return new_op.compose(self, qargs=layout)
 
 
+def sparsify_label(pauli_string):
+    """Return a sparse format of a Pauli string, e.g. "XIIIZ" -> ("XZ", [0, 4])."""
+    qubits = [i for i, label in enumerate(reversed(pauli_string)) if label != "I"]
+    sparse_label = "".join(pauli_string[~i] for i in qubits)
+
+    return sparse_label, qubits
+
+
 # Update docstrings for API docs
 generate_apidocs(SparsePauliOp)

qiskit/synthesis/evolution/lie_trotter.py

@@ -14,18 +14,15 @@
 
 from __future__ import annotations
 
-import inspect
 from collections.abc import Callable
 from typing import Any
-import numpy as np
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
-from qiskit.utils.deprecation import deprecate_arg
 
-from .product_formula import ProductFormula
+from .suzuki_trotter import SuzukiTrotter
 
 
-class LieTrotter(ProductFormula):
+class LieTrotter(SuzukiTrotter):
     r"""The Lie-Trotter product formula.
 
     The Lie-Trotter formula approximates the exponential of two non-commuting operators
@@ -40,7 +37,7 @@ class LieTrotter(ProductFormula):
 
     .. math::
 
-        e^{-it(XX + ZZ)} = e^{-it XX}e^{-it ZZ} + \mathcal{O}(t^2).
+        e^{-it(XI + ZZ)} = e^{-it XI}e^{-it ZZ} + \mathcal{O}(t^2).
 
     References:
 
@@ -52,21 +49,6 @@ class LieTrotter(ProductFormula):
         `arXiv:math-ph/0506007 <https://arxiv.org/pdf/math-ph/0506007.pdf>`_
     """
 
-    @deprecate_arg(
-        name="atomic_evolution",
-        since="1.2",
-        predicate=lambda callable: callable is not None
-        and len(inspect.signature(callable).parameters) == 2,
-        deprecation_description=(
-            "The 'Callable[[Pauli | SparsePauliOp, float], QuantumCircuit]' signature of the "
-            "'atomic_evolution' argument"
-        ),
-        additional_msg=(
-            "Instead you should update your 'atomic_evolution' function to be of the following "
-            "type: 'Callable[[QuantumCircuit, Pauli | SparsePauliOp, float], None]'."
-        ),
-        pending=True,
-    )
     def __init__(
         self,
         reps: int = 1,
@@ -100,26 +82,6 @@ def __init__(
         """
         super().__init__(1, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
 
-    def synthesize(self, evolution):
-        # get operators and time to evolve
-        operators = evolution.operator
-        time = evolution.time
-
-        # construct the evolution circuit
-        single_rep = QuantumCircuit(operators[0].num_qubits)
-
-        if not isinstance(operators, list):
-            pauli_list = [(Pauli(op), np.real(coeff)) for op, coeff in operators.to_list()]
-        else:
-            pauli_list = [(op, 1) for op in operators]
-
-        for i, (op, coeff) in enumerate(pauli_list):
-            self.atomic_evolution(single_rep, op, coeff * time / self.reps)
-            if self.insert_barriers and i != len(pauli_list) - 1:
-                single_rep.barrier()
-
-        return single_rep.repeat(self.reps, insert_barriers=self.insert_barriers).decompose()
-
     @property
     def settings(self) -> dict[str, Any]:
         """Return the settings in a dictionary, which can be used to reconstruct the object.

qiskit/synthesis/evolution/product_formula.py

@@ -6,7 +6,7 @@
 # obtain a copy of this license in the LICENSE.txt file in the root directory
 # of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
 #
-# Any modifications or derivative works of this code must retain this
+# typing.Any modifications or derivative works of this code must retain this
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
@@ -16,8 +16,8 @@
 
 import inspect
 from collections.abc import Callable
-from typing import Any
 from functools import partial
+import typing
 import numpy as np
 from qiskit.circuit.parameterexpression import ParameterExpression
 from qiskit.circuit.quantumcircuit import QuantumCircuit
@@ -26,6 +26,9 @@
 
 from .evolution_synthesis import EvolutionSynthesis
 
+if typing.TYPE_CHECKING:
+    from qiskit.circuit.library import PauliEvolutionGate
+
 
 class ProductFormula(EvolutionSynthesis):
     """Product formula base class for the decomposition of non-commuting operator exponentials.
@@ -108,8 +111,49 @@ def wrap_atomic_evolution(output, operator, time):
         else:
             self.atomic_evolution = atomic_evolution
 
+    def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], float]]:
+        """Apply the product formula to expand the Hamiltonian in the evolution gate.
+
+        Args:
+            evolution: The :class:`.PauliEvolutionGate`, whose Hamiltonian we expand.
+
+        Returns:
+            A list of Pauli rotations in a sparse format, where each element is
+            ``(paulistring, qubits, coefficient)``. For example, the Lie-Trotter expansion
+            of ``H = XI + ZZ`` would return ``[("X", [1], 1), ("ZZ", [0, 1], 1)]``.
+        """
+        raise NotImplementedError(
+            f"The method ``expand`` is not implemented for {self.__class__}. Implement it to "
+            f"automatically enable the call to {self.__class__}.synthesize."
+        )
+
+    def synthesize(self, evolution: PauliEvolutionGate) -> QuantumCircuit:
+        """Synthesize a :class:`.PauliEvolutionGate`.
+
+        Args:
+            evolution: The evolution gate to synthesize.
+
+        Returns:
+            QuantumCircuit: A circuit implementing the evolution.
+        """
+        pauli_rotations = self.expand(evolution)
+
+        num_qubits = evolution.num_qubits
+
+        circuit = QuantumCircuit(num_qubits)
+
+        # we could cache the circuit decomposition for the circuits we already saw
+        for i, (pauli_string, qubits, coeff) in enumerate(pauli_rotations):
+            op = SparsePauliOp.from_sparse_list([(pauli_string, qubits, 1)], num_qubits).paulis[0]
+
+            self.atomic_evolution(circuit, op, np.real_if_close(coeff))
+            if self.insert_barriers and i != len(pauli_rotations) - 1:
+                circuit.barrier()
+
+        return circuit
+
     @property
-    def settings(self) -> dict[str, Any]:
+    def settings(self) -> dict[str, typing.Any]:
         """Return the settings in a dictionary, which can be used to reconstruct the object.
 
         Returns:

qiskit/synthesis/evolution/qdrift.py

@@ -16,14 +16,18 @@
 
 import inspect
 import math
+import typing
+from itertools import chain
 from collections.abc import Callable
 import numpy as np
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
 from qiskit.utils.deprecation import deprecate_arg
 
 from .product_formula import ProductFormula
-from .lie_trotter import LieTrotter
+
+if typing.TYPE_CHECKING:
+    from qiskit.circuit.library import PauliEvolutionGate
 
 
 class QDrift(ProductFormula):
@@ -88,44 +92,33 @@ def __init__(
         self.sampled_ops = None
         self.rng = np.random.default_rng(seed)
 
-    def synthesize(self, evolution):
-        # get operators and time to evolve
+    def expand(self, evolution: PauliEvolutionGate) -> list[tuple[str, tuple[int], float]]:
         operators = evolution.operator
-        time = evolution.time
+        time = evolution.time  # used to determine the number of gates
 
-        if not isinstance(operators, list):
-            pauli_list = [(Pauli(op), coeff) for op, coeff in operators.to_list()]
-            coeffs = [np.real(coeff) for op, coeff in operators.to_list()]
+        # QDrift is based on first-order Lie-Trotter, hence we can just concatenate all
+        # Pauli terms and ignore commutations
+        if isinstance(operators, list):
+            paulis = list(chain.from_iterable([op.to_sparse_list() for op in operators]))
         else:
-            pauli_list = [(op, 1) for op in operators]
-            coeffs = [1 for op in operators]
+            paulis = operators.to_sparse_list()
+
+        coeffs = [np.real(coeff) for _, _, coeff in paulis]
 
         # We artificially make the weights positive
         weights = np.abs(coeffs)
         lambd = np.sum(weights)
 
         num_gates = math.ceil(2 * (lambd**2) * (time**2) * self.reps)
+
         # The protocol calls for the removal of the individual coefficients,
         # and multiplication by a constant evolution time.
-        evolution_time = lambd * time / num_gates
-
-        self.sampled_ops = self.rng.choice(
-            np.array(pauli_list, dtype=object),
-            size=(num_gates,),
-            p=weights / lambd,
-        )
-
-        # pylint: disable=cyclic-import
-        from qiskit.circuit.library.pauli_evolution import PauliEvolutionGate
-
-        # Build the evolution circuit using the LieTrotter synthesis with the sampled operators
-        lie_trotter = LieTrotter(
-            insert_barriers=self.insert_barriers, atomic_evolution=self.atomic_evolution
+        sampled = self.rng.choice(
+            np.array(paulis, dtype=object), size=(num_gates,), p=weights / lambd
         )
-        evolution_circuit = PauliEvolutionGate(
-            sum(SparsePauliOp(np.sign(coeff) * op) for op, coeff in self.sampled_ops),
-            time=evolution_time,
-            synthesis=lie_trotter,
-        ).definition
 
-        return evolution_circuit
+        rescaled_time = lambd * time / num_gates
+        sampled_paulis = [
+            (pauli[0], pauli[1], np.sign(pauli[2]) * rescaled_time) for pauli in sampled
+        ]
+        return sampled_paulis

qiskit/synthesis/evolution/suzuki_trotter.py

@@ -16,8 +16,7 @@
 
 import inspect
 from collections.abc import Callable
-
-import numpy as np
+from itertools import chain
 
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.quantum_info.operators import SparsePauliOp, Pauli
@@ -44,7 +43,7 @@ class SuzukiTrotter(ProductFormula):
 
     .. math::
 
-        e^{-it(XX + ZZ)} = e^{-it/2 ZZ}e^{-it XX}e^{-it/2 ZZ} + \mathcal{O}(t^3).
+        e^{-it(XI + ZZ)} = e^{-it/2 XI}e^{-it ZZ}e^{-it/2 XI} + \mathcal{O}(t^3).
 
     References:
         [1]: D. Berry, G. Ahokas, R. Cleve and B. Sanders,
@@ -105,51 +104,67 @@ def __init__(
             ValueError: If order is not even
         """
 
-        if order % 2 == 1:
+        if order > 1 and order % 2 == 1:
             raise ValueError(
                 "Suzuki product formulae are symmetric and therefore only defined "
-                "for even orders."
+                "for even orders (or order==1)."
             )
         super().__init__(order, reps, insert_barriers, cx_structure, atomic_evolution, wrap)
 
-    def synthesize(self, evolution):
-        # get operators and time to evolve
-        operators = evolution.operator
-        time = evolution.time
-
-        if not isinstance(operators, list):
-            pauli_list = [(Pauli(op), np.real(coeff)) for op, coeff in operators.to_list()]
-        else:
-            pauli_list = [(op, 1) for op in operators]
+    def expand(self, evolution):
+        """
+        H = ZZ + IX --> ("X", [0], 1/2), ("ZZ", [0, 1], 1), ("X", [0], 1/2)
 
-        ops_to_evolve = self._recurse(self.order, time / self.reps, pauli_list)
+        ("X", [0], 1/2), ("ZZ", [0, 1], 1), ("X", [0], 1), ("ZZ", [0, 1], 1), ("X", [0], 1/2)
+        """
+        operators = evolution.operator  # type: SparsePauliOp | list[SparsePauliOp]
+        time = evolution.time
 
         # construct the evolution circuit
-        single_rep = QuantumCircuit(operators[0].num_qubits)
-
-        for i, (op, coeff) in enumerate(ops_to_evolve):
-            self.atomic_evolution(single_rep, op, coeff)
-            if self.insert_barriers and i != len(ops_to_evolve) - 1:
-                single_rep.barrier()
+        if isinstance(operators, list):  # already sorted into commuting bits
+            non_commuting = [
+                (time / self.reps * operator).to_sparse_list() 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 (time / self.reps * operators).to_sparse_list()]
 
-        return single_rep.repeat(self.reps, insert_barriers=self.insert_barriers).decompose()
+        # we're already done here since Lie Trotter does not do any operator repetition
+        product_formula = self._recurse(self.order, non_commuting)
+        flattened = self.reps * list(chain.from_iterable(product_formula))
+        return flattened
 
     @staticmethod
-    def _recurse(order, time, pauli_list):
+    def _recurse(order, grouped_paulis):
         if order == 1:
-            return pauli_list
+            return grouped_paulis
 
         elif order == 2:
-            halves = [(op, coeff * time / 2) for op, coeff in pauli_list[:-1]]
-            full = [(pauli_list[-1][0], time * pauli_list[-1][1])]
+            halves = [
+                [(label, qubits, coeff / 2) for label, qubits, coeff in paulis]
+                for paulis in grouped_paulis[:-1]
+            ]
+            full = [grouped_paulis[-1]]
             return halves + full + list(reversed(halves))
 
         else:
             reduction = 1 / (4 - 4 ** (1 / (order - 1)))
             outer = 2 * SuzukiTrotter._recurse(
-                order - 2, time=reduction * time, pauli_list=pauli_list
+                order - 2,
+                [
+                    [(label, qubits, coeff * reduction) for label, qubits, coeff in paulis]
+                    for paulis in grouped_paulis
+                ],
             )
             inner = SuzukiTrotter._recurse(
-                order - 2, time=(1 - 4 * reduction) * time, pauli_list=pauli_list
+                order - 2,
+                [
+                    [
+                        (label, qubits, coeff * (1 - 4 * reduction))
+                        for label, qubits, coeff in paulis
+                    ]
+                    for paulis in grouped_paulis
+                ],
             )
             return outer + inner + outer