VPE circuit generation (#688)

* Added circuit generation functions, made some small shifts/edits to estimator functions

* Formatting

* Linting

Co-authored-by: Nicholas Rubin <rubinnc0@gmail.com>

src/openfermion/__init__.py

@@ -258,6 +258,9 @@
     error_bound,
     error_operator,
     trotter_steps_required,
+    vpe_single_circuit,
+    vpe_circuits_single_timestep,
+    standard_vpe_rotation_set,
 )
 
 from openfermion.testing import (

src/openfermion/circuits/__init__.py

@@ -101,3 +101,9 @@
     error_operator,
     trotter_steps_required,
 )
+
+from .vpe_circuits import (
+    vpe_single_circuit,
+    vpe_circuits_single_timestep,
+    standard_vpe_rotation_set,
+)

src/openfermion/circuits/vpe_circuits.py

@@ -0,0 +1,101 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""Circuit generation functions for verified phase estimation (2010.02538)"""
+
+from typing import Sequence, Optional
+import numpy
+import cirq
+
+
+def vpe_single_circuit(qubits: Sequence[cirq.Qid], prep: cirq.Circuit,
+                       evolve: cirq.Circuit, initial_rotation: cirq.Gate,
+                       final_rotation: cirq.Gate) -> cirq.Circuit:
+    """
+    Combines the different parts that make up a VPE circuit
+
+    The protocol for VPE requires combining preparation, evolution, and
+    measurement circuits for different values of time in order to estimate
+    the phase function. This function takes these parts and combines them.
+
+    Note that we need not specify the time of evolution as this is contained
+    already within evolve.
+
+    Arguments:
+        prep [cirq.Circuit] -- The circuit to prepare the initial state
+            (|psi_s>+|psi_r>) from |0>+|1>
+        evolve [cirq.Circuit] -- The circuit to evolve for time t
+        initial_rotation [cirq.Gate] -- The initial rotation on the target qubit
+            (Note that the gate should already be targeting the qubit)
+        final_rotation [cirq.Gate] -- The final rotation on the target qubit
+            (Note that the gate should already be targeting the qubit)
+    """
+    circuit = cirq.Circuit()
+    circuit.append(initial_rotation)
+    circuit.append(prep)
+    circuit.append(evolve)
+    circuit.append(cirq.inverse(prep))
+    circuit.append(final_rotation)
+    circuit.append(cirq.measure(*qubits, key='msmt'))
+    return circuit
+
+
+# Turning off yapf here as its formatting suggestion is bad.
+# yapf: disable
+standard_vpe_rotation_set = [
+    [0.25, cirq.ry(numpy.pi / 2), cirq.ry(-numpy.pi / 2)],
+    [-0.25, cirq.ry(numpy.pi / 2), cirq.ry(numpy.pi / 2)],
+    [-0.25j, cirq.ry(numpy.pi / 2), cirq.rx(-numpy.pi / 2)],
+    [0.25j, cirq.ry(numpy.pi / 2), cirq.rx(numpy.pi / 2)],
+    [0.25, cirq.rx(numpy.pi / 2), cirq.rx(-numpy.pi / 2)],
+    [-0.25, cirq.rx(numpy.pi / 2), cirq.rx(numpy.pi / 2)],
+    [0.25j, cirq.rx(numpy.pi / 2), cirq.ry(-numpy.pi / 2)],
+    [-0.25j, cirq.rx(numpy.pi / 2), cirq.ry(numpy.pi / 2)],
+]
+# yapf: enable
+
+
+def vpe_circuits_single_timestep(qubits: Sequence[cirq.Qid],
+                                 prep: cirq.Circuit,
+                                 evolve: cirq.Circuit,
+                                 target_qubit: cirq.Qid,
+                                 rotation_set: Optional[Sequence] = None
+                                ) -> Sequence[cirq.Circuit]:
+    """Prepares the circuits to perform VPE at a fixed time
+
+    Puts together the set of pre- and post-rotations to implement
+    VPE at for a given state preparation and time evolution.
+
+    [description]
+
+    Arguments:
+        prep [cirq.Circuit] -- The circuit to prepare the target state
+            (|psi_s>+|psi_r>) from |0>+|1>
+        evolve [cirq.Circuit] -- The circuit to evolve for time t
+        target_qubit [cirq.Qid] -- The qubit on which the phase
+            function is encoded
+        rotation_set [Sequence] -- A set of initial and final rotations for the
+            target qubit. We average the phase function estimation over multiple
+            such rotations to cancel out readout noise, final T1 decay, etc.
+            The standard rotation set is typically sufficient for these
+            purposes. The first element of each gate is the multiplier to get
+            the phase function; we do not need this for this function.
+
+            If rotation_set is set to None, the 'standard rotation set' of all
+            possible X and Y rotations before and after the circuit is used.
+    """
+    if rotation_set is None:
+        rotation_set = standard_vpe_rotation_set
+    circuits = [
+        vpe_single_circuit(qubits, prep, evolve, rdata[1].on(target_qubit),
+                           rdata[2].on(target_qubit)) for rdata in rotation_set
+    ]
+    return circuits

src/openfermion/circuits/vpe_circuits_test.py

@@ -0,0 +1,62 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""Tests for vpe_circuits.py"""
+
+import numpy
+import cirq
+
+from openfermion.measurements import get_phase_function
+
+from .vpe_circuits import (
+    vpe_single_circuit,
+    vpe_circuits_single_timestep,
+)
+
+
+def test_single_circuit():
+    q0 = cirq.GridQubit(0, 0)
+    q1 = cirq.GridQubit(0, 1)
+    qubits = reversed([q0, q1])
+    prep = cirq.Circuit([cirq.FSimGate(theta=numpy.pi / 4, phi=0).on(q0, q1)])
+    evolve = cirq.Circuit(
+        [cirq.rz(numpy.pi / 2).on(q0),
+         cirq.rz(numpy.pi / 2).on(q1)])
+    initial_rotation = cirq.ry(numpy.pi / 2).on(q0)
+    final_rotation = cirq.rx(-numpy.pi / 2).on(q0)
+    circuit = vpe_single_circuit(qubits, prep, evolve, initial_rotation,
+                                 final_rotation)
+    assert len(circuit) == 6
+    simulator = cirq.Simulator()
+    result = simulator.run(circuit, repetitions=100)
+    data_counts = result.data['msmt'].value_counts()
+    assert data_counts[1] == 100
+
+
+def test_single_timestep():
+    q0 = cirq.GridQubit(0, 0)
+    q1 = cirq.GridQubit(0, 1)
+    qubits = [q0, q1]
+    prep = cirq.Circuit([cirq.FSimGate(theta=numpy.pi / 4, phi=0).on(q0, q1)])
+    evolve = cirq.Circuit(
+        [cirq.rz(numpy.pi / 2).on(q0),
+         cirq.rz(numpy.pi / 2).on(q1)])
+    target_qubit = q0
+    circuits = vpe_circuits_single_timestep(qubits, prep, evolve, target_qubit)
+    results = []
+    simulator = cirq.Simulator()
+    for circuit in circuits:
+        this_res = simulator.run(circuit, repetitions=10000)
+        results.append(this_res)
+    pf_estimation = get_phase_function(results, qubits, 0)
+    # Loose bound - standard deviation on this measurement should be 0.005, so
+    # 95% of times this should pass. If it fails, check the absolute number.
+    assert numpy.abs(pf_estimation - 1j) < 1e-2

src/openfermion/measurements/vpe_estimators.py

@@ -18,6 +18,7 @@
 import cirq
 
 from openfermion.linalg import fit_known_frequencies
+from openfermion.circuits import standard_vpe_rotation_set
 
 
 class _VPEEstimator(metaclass=abc.ABCMeta):
@@ -74,7 +75,7 @@ def __init__(self, evals: numpy.ndarray, ref_eval: float = 0):
         self.evals = evals
         self.ref_eval = ref_eval
 
-    def get_simulation_points(self) -> numpy.ndarray:
+    def get_simulation_points(self, safe: bool = True) -> numpy.ndarray:
         """Generates time points for estimation
 
         VPE requires estimating the phase function g(t) at multiple points t,
@@ -83,15 +84,28 @@ def get_simulation_points(self) -> numpy.ndarray:
 
         In this case, we fit len(self.energies) complex amplitudes to a complex
         valued signal, we need precisely this number of points in the signal.
+
+        However, it appears numerically that approximately twice as many points
+        are needed to prevent aliasing, so we double this number here.
+
         Then, to prevent aliasing, we need to make sure that the time step
         dt < 2*pi / (E_max-E_min). Here, we choose dt = pi / (E_max-E_min).
         (Importantly, for Pauli operators this reproduces the H test.)
 
+        Args:
+            safe [bool, default True] -- numerical testing shows that taking
+                approximately twice as many points is better for the stability
+                of the estimator; this
+
         Returns:
             times: a set of times t that g(t) should be estimated at.
         """
-        numsteps = len(self.evals)
-        step_size = numpy.pi / (max(self.evals) - min(self.evals))
+        if safe:
+            numsteps = len(self.evals) * 2
+            step_size = numpy.pi / (max(self.evals) - min(self.evals))
+        else:
+            numsteps = len(self.evals)
+            step_size = numpy.pi / (max(self.evals) - min(self.evals))
         maxtime = step_size * (numsteps - 1)
         times = numpy.linspace(0, maxtime, numsteps)
         return times
@@ -131,21 +145,6 @@ def get_expectation_value(self,
         return expectation_value
 
 
-# disabling yapf here as its proposed formatting decreases readability
-# yapf: disable
-standard_rotation_set = [
-    [0.25, cirq.ry(numpy.pi / 2), cirq.ry(-numpy.pi / 2)],
-    [-0.25, cirq.ry(numpy.pi / 2), cirq.ry(numpy.pi / 2)],
-    [-0.25j, cirq.ry(numpy.pi / 2), cirq.rx(-numpy.pi / 2)],
-    [0.25j, cirq.ry(numpy.pi / 2), cirq.rx(numpy.pi / 2)],
-    [0.25, cirq.rx(numpy.pi / 2), cirq.rx(-numpy.pi / 2)],
-    [-0.25, cirq.rx(numpy.pi / 2), cirq.rx(numpy.pi / 2)],
-    [0.25j, cirq.rx(numpy.pi / 2), cirq.ry(-numpy.pi / 2)],
-    [-0.25j, cirq.rx(numpy.pi / 2), cirq.ry(numpy.pi / 2)],
-]
-# yapf: enable
-
-
 def get_phase_function(results: Sequence[cirq.TrialResult],
                        qubits: Sequence[cirq.Qid],
                        target_qid: int,
@@ -185,7 +184,7 @@ def get_phase_function(results: Sequence[cirq.TrialResult],
     """
     hs_index = 2**(len(qubits) - target_qid - 1)
     if rotation_set is None:
-        rotation_set = standard_rotation_set
+        rotation_set = standard_vpe_rotation_set
     phase_function = 0
     if len(results) != len(rotation_set):
         raise ValueError("I have an incorrect number of TrialResults "

src/openfermion/measurements/vpe_estimators_test.py

@@ -26,7 +26,7 @@
 
 def test_requests_simulation_at_pi_for_pauli():
     estimator = PhaseFitEstimator(evals=[-1, +1])
-    sim_points = estimator.get_simulation_points()
+    sim_points = estimator.get_simulation_points(safe=False)
     assert len(sim_points) == 2
     assert numpy.isclose(sim_points[0], 0)
     assert numpy.isclose(sim_points[1], numpy.pi / 2)