Derandomized + Adaptive Classical Shadows

tangelo/toolboxes/measurements/__init__.py

@@ -14,3 +14,7 @@
 
 from .qubit_terms_grouping import group_qwc, exp_value_from_measurement_bases
 from .estimate_measurements import get_measurement_estimate
+from .classical_shadows.classical_shadows import zero_state, one_state, I, rotations, matrices, traces, ClassicalShadow
+from .classical_shadows.randomized import RandomizedClassicalShadow
+from .classical_shadows.derandomized import DerandomizedClassicalShadow
+from .classical_shadows.adaptive import AdaptiveClassicalShadow

tangelo/toolboxes/measurements/classical_shadows/adaptive.py

@@ -0,0 +1,201 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# 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.
+
+"""This file provides an API enabling the use of adaptive classical shadows.
+This algorithm is described in C. Hadfield, ArXiv:2105.12207 [Quant-Ph] (2021).
+"""
+
+import random
+
+import numpy as np
+
+from tangelo.toolboxes.measurements import ClassicalShadow
+from tangelo.linq.circuit import Circuit
+from tangelo.linq.helpers.circuits.measurement_basis import measurement_basis_gates, pauli_string_to_of
+
+
+class AdaptiveClassicalShadow(ClassicalShadow):
+    """Classical shadows using adaptive single Pauli measurements, as defined
+    in C. Hadfield, ArXiv:2105.12207 [Quant-Ph] (2021).
+    """
+
+    def build(self, n_shots, qu_op):
+        """Adaptive classical shadow building method to define relevant
+        unitaries depending on the qubit operator.
+
+        Args:
+            n_shots (int): The number of desired measurements.
+            qu_op (QubitOperator): The observable that one wishes to measure.
+
+        Returns:
+            list of str: The list of Pauli words that describes the measurement
+                basis to use.
+        """
+
+        measurement_procedure = [self._choose_measurement(qu_op) for _ in range(n_shots)]
+
+        self.unitaries = measurement_procedure
+        return measurement_procedure
+
+    def _choose_measurement(self, qu_op):
+        """Algorithm 1 from the publication.
+
+        Args:
+            qu_op (QubitOperator): The operator that one wishes to maximize the
+                measurement budget over.
+
+        Returns:
+            str: Pauli words for one measurement.
+        """
+
+        i_qubit_random = random.sample(range(self.n_qubits), self.n_qubits)
+        inverse_map = np.argsort(i_qubit_random)
+
+        single_measurement = [None] * self.n_qubits
+
+        for iteration, i_qubit in enumerate(i_qubit_random):
+            cbs = self._generate_cbs(qu_op,
+                                     i_qubit_random[0:iteration],
+                                     single_measurement[0:iteration],
+                                     i_qubit)
+            sum_cbs = sum(cbs)
+
+            # If sum is 0., the distribution is set to be uniform.
+            if sum_cbs < 1.e-7:
+                single_measurement[iteration] = random.choice(["X", "Y", "Z"])
+            else:
+                random_val = random.random()
+
+                # Depending on the cbs, the probabilities of drawing a specific
+                # Pauli gate are shifted.
+                if random_val < cbs[0] / sum_cbs:
+                    single_measurement[iteration] = "X"
+                elif random_val < (cbs[0] + cbs[1]) / sum_cbs:
+                    single_measurement[iteration] = "Y"
+                else:
+                    single_measurement[iteration] = "Z"
+
+        # Reordering according to the qubit indices 0, 1, 2, ... self.n_qubits.
+        reordered_measurement = [single_measurement[inverse_map[j]] for j in range(self.n_qubits)]
+
+        return "".join(reordered_measurement)
+
+    def _generate_cbs(self, qu_op, prev_qubits, prev_paulis, curr_qubit):
+        """Generates the cB values from which the Pauli basis is determined for
+        the current qubit (curr_qubit), as shown in Algorithm 2 from the paper.
+
+        Args:
+            qu_op (QubitOperator) : The operator one wishes to get the
+                expectation value of.
+            prev_qubits (list) : list of previous qubits from which the
+                measurement basis is already determined
+            curr_qubit (int) : The current qubit being examined.
+            bs (list) : the Pauli word for prev_qubits.
+
+        Returns:
+            list of float: cB values for X, Y and Z.
+        """
+
+        cbs = [0.] * 3
+        map_pauli = {"X": 0, "Y": 1, "Z": 2}
+
+        for term, coeff in qu_op.terms.items():
+
+            # Default conditions to compute cbs.
+            same_qubit_flag = False
+            same_pauli_flag = True
+
+            for i_qubit, pauli in term:
+                # Checking if the current qubit index has been detected in the
+                # term.
+                if i_qubit == curr_qubit:
+                    same_qubit_flag = True
+
+                # Checking if the Pauli basis in the term has already been
+                # chosen for this qubit.
+                for prev_qubit, prev_pauli in zip(prev_qubits, prev_paulis):
+                    if i_qubit == prev_qubit and pauli != prev_pauli:
+                        same_pauli_flag = False
+
+                if same_qubit_flag and same_pauli_flag:
+                    cbs[map_pauli[pauli]] += abs(coeff)**2
+
+        return np.sqrt(cbs)
+
+    def get_basis_circuits(self, only_unique=False):
+        """Outputs a list of circuits corresponding to the adaptive single-Pauli
+        unitaries.
+
+        Args:
+            only_unique (bool): Considering only unique unitaries.
+
+        Returns:
+            list of Circuit or tuple: All basis circuits or a tuple of unique
+                circuits (first) with the numbers of occurence (last).
+        """
+
+        if not self.unitaries:
+            raise ValueError(f"A set of unitaries must de defined (can be done with the build method in {self.__class__.__name__}).")
+
+        unitaries_to_convert = self.unique_unitaries if only_unique else self.unitaries
+
+        basis_circuits = list()
+        for pauli_word in unitaries_to_convert:
+            # Transformation of a unitary to quantum gates.
+            pauli_of = pauli_string_to_of(pauli_word)
+            basis_circuits += [Circuit(measurement_basis_gates(pauli_of), self.n_qubits)]
+
+        # Counting each unique circuits (use for reversing to a full shadow from an experiement on hardware).
+        if only_unique:
+            unique_basis_circuits = [(basis_circuits[i], self.unitaries.count(u)) for i, u in enumerate(unitaries_to_convert)]
+            return unique_basis_circuits
+        else:
+            return basis_circuits
+
+    def get_term_observable(self, term, coeff=1.):
+        """Returns the estimated observable for a term and its coefficient.
+
+        Args:
+            term (tuple): Openfermion style of a qubit operator term.
+            coeff (float): Multiplication factor for the term.
+
+        Returns:
+            float: Observable estimated with the shadow.
+        """
+
+        sum_product = 0
+        n_match = 0
+
+        zero_state, one_state = 1, -1
+
+        # For every single_measurement in shadow_size.
+        for snapshot in range(self.size):
+            match = True
+            product = 1
+
+            # Checking if there is a match for all Pauli gate in the term.
+            # Works also with operator not on all qubits (e.g. X1 will hit Z0X1,
+            # Y0X1 and Z0X1).
+            for i_qubit, pauli in term:
+                if pauli != self.unitaries[snapshot][i_qubit]:
+                    match = False
+                    break
+                state = zero_state if self.bitstrings[snapshot][i_qubit] == "0" else one_state
+                product *= state
+
+            # No quantity is considered if there is no match.
+            sum_product += match*product
+            n_match += match
+
+        return sum_product / n_match * coeff if n_match > 0 else 0.

tangelo/toolboxes/measurements/classical_shadows/classical_shadows.py

@@ -0,0 +1,165 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# 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.
+
+"""This file provides an API enabling the use of classical shadows. The original
+idea is described in H.Y. Huang, R. Kueng, and J. Preskill, Nature Physics 16,
+1050 (2020).
+"""
+
+import abc
+import warnings
+
+import numpy as np
+
+from tangelo.linq.circuit import Circuit
+
+
+# State |0> or |1>.
+zero_state = np.array([1, 0])
+one_state = np.array([0, 1])
+
+# Pauli matrices.
+I = np.array([[1, 0], [0, 1]])
+X = np.array([[0, 1], [1, 0]])
+Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
+Z = np.array([[1, 0], [0, -1]])
+matrices = {"X": X, "Y": Y, "Z": Z}
+
+# Traces of each Pauli matrices.
+traces = {pauli: np.trace(matrix) for pauli, matrix in matrices.items()}
+
+# Reverse channels to undo single Pauli rotations.
+S_T = np.array([[1, 0], [0, -1j]], dtype=complex)
+H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
+I = np.array([[1, 0], [0, 1]])
+rotations = {"X": H, "Y": H @ S_T, "Z": I}
+
+
+class ClassicalShadow(abc.ABC):
+    """Abstract class for the classical shadows implementation. Classical
+    shadows is a mean to characterize a quantum state (within an error treshold)
+    with the fewest measurement possible.
+    """
+
+    def __init__(self, circuit, bitstrings=None, unitaries=None):
+        """Default constructor for the ClassicalShadow object. This class is
+        the parent class for the different classical shadows flavors. The object
+        is defined by the bistrings and unitaries used in the process. Abstract
+        methods are defined to take into account the procedure to inverse the
+        channel.
+
+        Args:
+            bistrings (list of str): Representation of the outcomes for all
+                snapshots. E.g. ["11011", "10000", ...].
+            unitaries (list of str): Representation of the unitary for every
+                snapshot, used to reverse the channel.
+        """
+
+        self.circuit = circuit
+        self.bitstrings = list() if bitstrings is None else bitstrings
+        self.unitaries = list() if unitaries is None else unitaries
+
+        # If the state has been estimated, it is stored into this attribute.
+        self.state_estimate = None
+
+    @property
+    def n_qubits(self):
+        """Returns the number of qubits the shadow represents."""
+        return self.circuit.width
+
+    @property
+    def size(self):
+        """Number of shots used to make the shadow."""
+        return len(self.bitstrings)
+
+    @property
+    def unique_unitaries(self):
+        """Returns the list of unique unitaries."""
+        return list(set(self.unitaries))
+
+    def __len__(self):
+        """Same as the shadow size."""
+        return self.size
+
+    def append(self, bitstring, unitary):
+        """Append method to merge new snapshots to an existing shadow.
+
+        Args:
+            bistring (str or list of str): Representation of outcomes.
+            unitary (str or list of str): Relevant unitary for those outcomes.
+        """
+        if isinstance(bitstring, list) and isinstance(unitary, list):
+            assert len(bitstring) == len(unitary)
+            self.bitstrings += bitstring
+            self.unitaries += unitary
+        elif isinstance(bitstring, str) and isinstance(unitary, str):
+            self.bitstrings.append(bitstring)
+            self.unitaries.append(unitary)
+        else:
+            raise ValueError("bistring and unitary arguments must be consistent strings or list of strings.")
+
+    def get_observable(self, qubit_op, *args, **kwargs):
+        """Getting an estimated observable value for a qubit operator from the
+        classical shadow. This function loops through all terms and calls, for
+        each of them, the get_term_observable method defined in the child class.
+        Other arguments (args, kwargs) can be passed to the method.
+
+        Args:
+            qubit_op (QubitOperator): Operator to estimate.
+        """
+        observable = 0.
+        for term, coeff in qubit_op.terms.items():
+            observable += self.get_term_observable(term, coeff, *args, **kwargs)
+
+        return observable
+
+    def simulate(self, backend, initial_statevector=None):
+        """Simulate, using a predefined backend, a shadow from a circuit or a
+        statevector.
+
+        Args:
+            backend (Simulator): Backend for the simulation of a shadow.
+            initial_statevector(list/array) : A valid statevector in the format
+                supported by the target backend.
+        """
+
+        if not self.unitaries:
+            raise ValueError(f"The build method of {self.__class__.__name__} must be called before simulation.")
+
+        if backend.n_shots != 1:
+            warnings.warn(f"Changing number of shots to 1 for the backend (classical shadows).")
+            backend.n_shots = 1
+
+        # Different behavior if circuit or initial_statevector is defined.
+        one_shot_circuit_template = self.circuit if self.circuit is not None else Circuit(n_qubits=self.n_qubits)
+
+        for basis_circuit in self.get_basis_circuits(only_unique=False):
+            one_shot_circuit = one_shot_circuit_template + basis_circuit if (basis_circuit.size > 0) else one_shot_circuit_template
+
+            # Frequencies returned by simulate are of the form {'0100...01': 1.0}.
+            # We add the bitstring to the shadow.
+            freqs, _ = backend.simulate(one_shot_circuit, initial_statevector=initial_statevector)
+            self.bitstrings += [list(freqs.keys())[0]]
+
+    @abc.abstractmethod
+    def build(self):
+        pass
+
+    @abc.abstractmethod
+    def get_basis_circuits(self, only_unique=False):
+        pass
+
+    @abc.abstractmethod
+    def get_term_observable(self):
+        pass