Implement Qdrift as ProductFormula (#7281)

* cast qdrift as a product formula

* release note and minor style changes

* add test and modify error message

* change evo to evolution_circuit and delete suzuki check

* fix lint

* fixes from PR comments

* remove sampled_ops property

* seed the random generator

* Update qiskit/synthesis/evolution/qdrift.py

Co-authored-by: Julien Gacon <gaconju@gmail.com>

* fix constant coefficients and cyclic import

* fix small mistake in the construction of pauli_list

* remove old leftover comment

* num_gates should be int & add functional test

* update test and evolution time in sampled_ops

Co-authored-by: Julien Gacon <gaconju@gmail.com>
Co-authored-by: mergify[bot] <37929162+mergify[bot]@users.noreply.github.com>
This commit is contained in:
Almudena Carrera Vazquez 2021-11-29 18:30:05 +01:00 committed by GitHub
parent a0e2a29a3a
commit 4a94b770e3
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 187 additions and 20 deletions

View File

@ -37,4 +37,5 @@ from .evolution import (
LieTrotter,
SuzukiTrotter,
MatrixExponential,
QDrift,
)

View File

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

View File

@ -71,7 +71,7 @@ class LieTrotter(ProductFormula):
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 @@ class LieTrotter(ProductFormula):
# 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

View File

@ -33,7 +33,7 @@ class MatrixExponential(EvolutionSynthesis):
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 @@ class MatrixExponential(EvolutionSynthesis):
exp = expm(-1j * time * matrix)
evo.unitary(exp, evo.qubits)
evolution_circuit.unitary(exp, evolution_circuit.qubits)
return evo
return evolution_circuit

View File

@ -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

View File

@ -0,0 +1,100 @@
# 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 Union, Optional, Callable
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
from .lie_trotter import LieTrotter
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
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), 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
weights = np.abs(coeffs)
lambd = np.sum(weights)
num_gates = int(np.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 = algorithm_globals.random.choice(
np.array(pauli_list, dtype=object),
size=(num_gates,),
p=weights / lambd,
)
# Update the coefficients of sampled_ops
self.sampled_ops = [(op, evolution_time) for op, coeff in self.sampled_ops]
# pylint: disable=cyclic-import
from qiskit.circuit.library.pauli_evolution import PauliEvolutionGate
from qiskit.opflow import PauliOp
# 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
)
evolution_circuit = PauliEvolutionGate(
sum(PauliOp(op) for op, coeff in self.sampled_ops),
time=evolution_time,
synthesis=lie_trotter,
).definition
return evolution_circuit

View File

@ -84,7 +84,7 @@ class SuzukiTrotter(ProductFormula):
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 @@ class SuzukiTrotter(ProductFormula):
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

View File

@ -0,0 +1,28 @@
---
features:
- |
Reformulate the former :class:`~qiskit.opflow.evolutions.trotterizations.QDrift`
as a synthesis method instead of a :obj:`~qiskit.opflow.evolutions.TrotterizationBase`. The
other available synthesis methods are
* :class:`~qiskit.synthesis.LieTrotter` - first order Trotterization
* :class:`~qiskit.synthesis.SuzukiTrotter` - higher order Trotterization
* :class:`~qiskit.synthesis.MatrixExponentiation` - exact, matrix-based evolution
.. code-block:: python
from qiskit.circuit import QuantumCircuit
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import QDrift
qdrift = QDrift(reps=2)
operator = (X ^ 3) + (Y ^ 3) + (Z ^ 3)
time = 2.345 # evolution time
evolution_gate = PauliEvolutionGate(operator, time, synthesis=qdrift)
circuit = QuantumCircuit(3)
circuit.append(evolution_gate, range(3))
print(circuit.draw())

View File

@ -12,22 +12,29 @@
"""Test the evolution gate."""
import numpy as np
import scipy
from ddt import ddt, data
from ddt import ddt, data, 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
from qiskit.quantum_info import Operator, SparsePauliOp, Pauli
from qiskit.quantum_info import Operator, SparsePauliOp, Pauli, Statevector
from qiskit.utils import algorithm_globals
@ddt
class TestEvolutionGate(QiskitTestCase):
"""Test the evolution gate."""
def setUp(self):
super().setUp()
# fix random seed for reproducibility (used in QDrift)
algorithm_globals.random_seed = 2
def test_matrix_decomposition(self):
"""Test the default decomposition."""
op = (X ^ 3) + (Y ^ 3) + (Z ^ 3)
@ -113,6 +120,39 @@ class TestEvolutionGate(QiskitTestCase):
self.assertEqual(evo_gate.definition.decompose(), expected)
@data(
(X + Y, 0.5, 1, [(Pauli("X"), 0.5), (Pauli("X"), 0.5)]),
(X, 0.238, 2, [(Pauli("X"), 0.238)]),
)
@unpack
def test_qdrift_manual(self, op, time, reps, sampled_ops):
"""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 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.assertTrue(Operator(evo_gate.definition).equiv(expected))
def test_qdrift_evolution(self):
"""Test QDrift on an example."""
op = 0.1 * (Z ^ Z) + (X ^ I) + (I ^ X) + 0.2 * (X ^ X)
reps = 20
qdrift = PauliEvolutionGate(op, time=0.5 / reps, synthesis=QDrift(reps=reps)).definition
exact = scipy.linalg.expm(-0.5j * op.to_matrix()).dot(np.eye(4)[0, :])
def energy(evo):
return Statevector(evo).expectation_value(op.to_matrix())
self.assertAlmostEqual(energy(exact), energy(qdrift), places=2)
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)]