Implement Qdrift as ProductFormula

qiskit/synthesis/__init__.py

@@ -37,4 +37,5 @@
     LieTrotter,
     SuzukiTrotter,
     MatrixExponential,
+    QDrift,
 )

qiskit/synthesis/evolution/__init__.py

@@ -17,3 +17,4 @@
 from .product_formula import ProductFormula
 from .lie_trotter import LieTrotter
 from .suzuki_trotter import SuzukiTrotter
+from .qdrift import QDrift

qiskit/synthesis/evolution/lie_trotter.py

@@ -71,7 +71,7 @@ def synthesize(self, evolution):
         time = evolution.time
 
         # construct the evolution circuit
-        evo = QuantumCircuit(operators[0].num_qubits)
+        evolution_circuit = QuantumCircuit(operators[0].num_qubits)
         first_barrier = False
 
         if not isinstance(operators, list):
@@ -87,12 +87,12 @@ def synthesize(self, evolution):
                 # add barriers
                 if first_barrier:
                     if self.insert_barriers:
-                        evo.barrier()
+                        evolution_circuit.barrier()
                 else:
                     first_barrier = True
 
-                evo.compose(
+                evolution_circuit.compose(
                     self.atomic_evolution(op, coeff * time / self.reps), wrap=wrap, inplace=True
                 )
 
-        return evo
+        return evolution_circuit

qiskit/synthesis/evolution/matrix_synthesis.py

@@ -33,7 +33,7 @@ def synthesize(self, evolution):
         time = evolution.time
 
         # construct the evolution circuit
-        evo = QuantumCircuit(operators[0].num_qubits)
+        evolution_circuit = QuantumCircuit(operators[0].num_qubits)
 
         if not isinstance(operators, list):
             matrix = operators.to_matrix()
@@ -42,6 +42,6 @@ def synthesize(self, evolution):
 
         exp = expm(-1j * time * matrix)
 
-        evo.unitary(exp, evo.qubits)
+        evolution_circuit.unitary(exp, evolution_circuit.qubits)
 
-        return evo
+        return evolution_circuit

qiskit/synthesis/evolution/product_formula.py

@@ -283,13 +283,13 @@ def cnot_fountain(pauli: Pauli) -> QuantumCircuit:
 def _default_atomic_evolution(operator, time, cx_structure):
     if isinstance(operator, Pauli):
         # single Pauli operator: just exponentiate it
-        evo = evolve_pauli(operator, time, cx_structure)
+        evolution_circuit = evolve_pauli(operator, time, cx_structure)
     else:
         # sum of Pauli operators: exponentiate each term (this assumes they commute)
         pauli_list = [(Pauli(op), np.real(coeff)) for op, coeff in operator.to_list()]
         name = f"exp(it {[pauli.to_label() for pauli, _ in pauli_list]})"
-        evo = QuantumCircuit(operator.num_qubits, name=name)
+        evolution_circuit = QuantumCircuit(operator.num_qubits, name=name)
         for pauli, coeff in pauli_list:
-            evo.compose(evolve_pauli(pauli, coeff * time, cx_structure), inplace=True)
+            evolution_circuit.compose(evolve_pauli(pauli, coeff * time, cx_structure), inplace=True)
 
-    return evo
+    return evolution_circuit

qiskit/synthesis/evolution/qdrift.py

@@ -0,0 +1,110 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# 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
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""QDrift Class"""
+
+from typing import List, Union, Optional, Callable, Tuple
+import numpy as np
+from qiskit.circuit.quantumcircuit import QuantumCircuit
+from qiskit.quantum_info.operators import SparsePauliOp, Pauli
+from qiskit.utils import algorithm_globals
+
+from .product_formula import ProductFormula
+
+# pylint: disable=invalid-name
+
+
+class QDrift(ProductFormula):
+    r"""The QDrift Trotterization method, which selects each each term in the
+    Trotterization randomly, with a probability proportional to its weight. Based on the work
+    of Earl Campbell in Ref. [1].
+
+    References:
+        [1]: E. Campbell, "A random compiler for fast Hamiltonian simulation" (2018).
+         `arXiv:quant-ph/1811.08017 <https://arxiv.org/abs/1811.08017>`_
+    """
+
+    def __init__(
+        self,
+        reps: int = 1,
+        insert_barriers: bool = False,
+        cx_structure: str = "chain",
+        atomic_evolution: Optional[
+            Callable[[Union[Pauli, SparsePauliOp], float], QuantumCircuit]
+        ] = None,
+    ) -> None:
+        r"""
+        Args:
+            reps: The number of times to repeat the Trotterization circuit.
+            insert_barriers: Whether to insert barriers between the atomic evolutions.
+            cx_structure: How to arrange the CX gates for the Pauli evolutions, can be
+                "chain", where next neighbor connections are used, or "fountain", where all
+                qubits are connected to one.
+            atomic_evolution: A function to construct the circuit for the evolution of single
+                Pauli string. Per default, a single Pauli evolution is decomopsed in a CX chain
+                and a single qubit Z rotation.
+        """
+        super().__init__(1, reps, insert_barriers, cx_structure, atomic_evolution)
+        self.sampled_ops = None
+
+    @property
+    def sampled_ops(self) -> List[Tuple[Pauli, float]]:
+        """returns the list of sampled Pauli ops and their coefficients"""
+        return self._sampled_ops
+
+    @sampled_ops.setter
+    def sampled_ops(self, sampled_ops: List[Tuple[Pauli, float]]) -> None:
+        """sets the list of sampled Pauli ops and their coefficients"""
+        self._sampled_ops = sampled_ops
+
+    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()]
+            coeffs = [np.real(coeff) for op, coeff in operators.to_list()]
+        else:
+            pauli_list = [(op, 1) for op in operators]
+            coeffs = [1 for op in operators]
+
+        # We artificially make the weights positive, TODO check approximation performance
+        weights = np.abs(coeffs)
+        lambd = np.sum(weights)
+
+        N = 2 * (lambd ** 2) * (time ** 2)
+        factor = lambd * time / N * self.reps
+        # The protocol calls for the removal of the individual coefficients,
+        # and multiplication by a constant factor.
+        scaled_ops = [(op, factor / coeff) for op, coeff in pauli_list]
+        self.sampled_ops = algorithm_globals.random.choice(
+            np.array(scaled_ops, dtype=object),
+            size=(int(np.ceil(N * self.reps)),),
+            p=weights / lambd,
+        )
+
+        # construct the evolution circuit
+        evolution_circuit = QuantumCircuit(operators[0].num_qubits)
+        first_barrier = False
+
+        for op, coeff in self.sampled_ops:
+            # add barriers
+            if first_barrier:
+                if self.insert_barriers:
+                    evolution_circuit.barrier()
+            else:
+                first_barrier = True
+
+            evolution_circuit.compose(self.atomic_evolution(op, coeff), wrap=True, inplace=True)
+
+        return evolution_circuit

qiskit/synthesis/evolution/suzuki_trotter.py

@@ -84,7 +84,7 @@ def synthesize(self, evolution):
 
             single_rep.compose(self.atomic_evolution(op, coeff), wrap=True, inplace=True)
 
-        evo = QuantumCircuit(operators[0].num_qubits)
+        evolution_circuit = QuantumCircuit(operators[0].num_qubits)
         first_barrier = False
 
         for _ in range(self.reps):
@@ -95,15 +95,12 @@ def synthesize(self, evolution):
             else:
                 first_barrier = True
 
-            evo.compose(single_rep, inplace=True)
+            evolution_circuit.compose(single_rep, inplace=True)
 
-        return evo
+        return evolution_circuit
 
     @staticmethod
     def _recurse(order, time, pauli_list):
-        if order < 1:
-            raise ValueError("This bitch empty -- yeet!")
-
         if order == 1:
             return pauli_list
 

releasenotes/notes/qdrift-as-product-formula-044a37a106a45a47.yaml

@@ -0,0 +1,6 @@
+---
+prelude: >
+    Reformulate the former QDrift class in qiskit.opflow.evolutions.trotterizations as a ProductFormula instead of a TrotterizationBase.
+other:
+  - |
+    Reformulate the former QDrift class in qiskit.opflow.evolutions.trotterizations as a ProductFormula instead of a TrotterizationBase.

test/python/circuit/library/test_evolution_gate.py

@@ -13,11 +13,11 @@
 """Test the evolution gate."""
 
 import scipy
-from ddt import ddt, data
+from ddt import ddt, data, idata, unpack
 
 from qiskit.circuit import QuantumCircuit, Parameter
 from qiskit.circuit.library import PauliEvolutionGate
-from qiskit.synthesis import LieTrotter, SuzukiTrotter, MatrixExponential
+from qiskit.synthesis import LieTrotter, SuzukiTrotter, MatrixExponential, QDrift
 from qiskit.converters import circuit_to_dag
 from qiskit.test import QiskitTestCase
 from qiskit.opflow import I, X, Y, Z, PauliSumOp
@@ -113,6 +113,29 @@ def test_suzuki_trotter_manual(self):
 
         self.assertEqual(evo_gate.definition.decompose(), expected)
 
+    @idata(
+        [
+            [X + Y, 0.5, 1],
+            [X, 0.238, 2],
+        ]
+    )
+    @unpack
+    def test_qdrift_manual(self, op, time, reps):
+        """Test the evolution circuit of Suzuki Trotter against a manually constructed circuit."""
+        qdrift = QDrift(reps=reps)
+        evo_gate = PauliEvolutionGate(op, time, synthesis=qdrift)
+        evo_gate.definition.decompose()
+
+        # manually construct expected evolution
+        expected = QuantumCircuit(1)
+        for pauli in qdrift.sampled_ops:
+            if pauli[0].to_label() == "X":
+                expected.rx(2 * pauli[1], 0)
+            elif pauli[0].to_label() == "Y":
+                expected.ry(2 * pauli[1], 0)
+
+        self.assertEqual(evo_gate.definition.decompose(), expected)
+
     def test_passing_grouped_paulis(self):
         """Test passing a list of already grouped Paulis."""
         grouped_ops = [(X ^ Y) + (Y ^ X), (Z ^ I) + (Z ^ Z) + (I ^ Z), (X ^ X)]