Repository: qiskit-community/povm-toolbox

Issue:
[feature] Optionally twirl measurements
It would be nice if the `RandomizedProjectiveMeasurement` (and by extension its subclasses) had built-in support for measurement twirling.

Discussion:

Patch:
diff --git a/povm_toolbox/library/classical_shadows.py b/povm_toolbox/library/classical_shadows.py
index 3bcb1d9..97f3f81 100644
--- a/povm_toolbox/library/classical_shadows.py
+++ b/povm_toolbox/library/classical_shadows.py
@@ -24,6 +24,7 @@ class ClassicalShadows(LocallyBiasedClassicalShadows):
     def __init__(
         self,
         n_qubit: int,
+        measurement_twirl: bool = False,
         measurement_layout: list[int] | None = None,  # TODO: add | Layout
         shot_batch_size: int = 1,
         seed_rng: int | Generator | None = None,
@@ -36,6 +37,10 @@ class ClassicalShadows(LocallyBiasedClassicalShadows):
 
         Args:
             n_qubits: the number of qubits.
+            measurement_twirl : option to randomly twirl the measurements. For each single-qubit
+                projective measurement, random twirling is equivalent to randomly flipping the
+                measurement. This is equivalent to randomly taking the opposite Bloch vector in
+                the Bloch sphere representation.
             measurement_layout: list of indices specifying on which qubits the POVM
                 acts. If None, two cases can be distinguished: 1) if a circuit supplied
                 to the :meth:`.compose_circuits` has been transpiled, its final
@@ -51,6 +56,7 @@ class ClassicalShadows(LocallyBiasedClassicalShadows):
         super().__init__(
             n_qubit=n_qubit,
             bias=bias,
+            measurement_twirl=measurement_twirl,
             measurement_layout=measurement_layout,
             shot_batch_size=shot_batch_size,
             seed_rng=seed_rng,
diff --git a/povm_toolbox/library/locally_biased_classical_shadows.py b/povm_toolbox/library/locally_biased_classical_shadows.py
index 97db5a9..d869bb1 100644
--- a/povm_toolbox/library/locally_biased_classical_shadows.py
+++ b/povm_toolbox/library/locally_biased_classical_shadows.py
@@ -25,6 +25,7 @@ class LocallyBiasedClassicalShadows(RandomizedProjectiveMeasurements):
         self,
         n_qubit: int,
         bias: np.ndarray,
+        measurement_twirl: bool = False,
         measurement_layout: list[int] | None = None,  # TODO: add | Layout
         shot_batch_size: int = 1,
         seed_rng: int | Generator | None = None,
@@ -40,6 +41,10 @@ class LocallyBiasedClassicalShadows(RandomizedProjectiveMeasurements):
                 for measuring in each of the PVMs. I.e., its length equals the number of PVMs (3).
                 These floats should sum to 1. If 2D, it will have a new set of biases for each
                 qubit.
+            measurement_twirl : option to randomly twirl the measurements. For each single-qubit
+                projective measurement, random twirling is equivalent to randomly flipping the
+                measurement. This is equivalent to randomly taking the opposite Bloch vector in
+                the Bloch sphere representation.
             measurement_layout: list of indices specifying on which qubits the POVM
                 acts. If None, two cases can be distinguished: 1) if a circuit supplied
                 to the :meth:`.compose_circuits` has been transpiled, its final
@@ -58,6 +63,7 @@ class LocallyBiasedClassicalShadows(RandomizedProjectiveMeasurements):
             n_qubit=n_qubit,
             bias=bias,
             angles=angles,
+            measurement_twirl=measurement_twirl,
             measurement_layout=measurement_layout,
             shot_batch_size=shot_batch_size,
             seed_rng=seed_rng,
diff --git a/povm_toolbox/library/randomized_projective_measurements.py b/povm_toolbox/library/randomized_projective_measurements.py
index b650942..46fc0e5 100644
--- a/povm_toolbox/library/randomized_projective_measurements.py
+++ b/povm_toolbox/library/randomized_projective_measurements.py
@@ -38,6 +38,7 @@ class RandomizedProjectiveMeasurements(POVMImplementation[RPMMetadata]):
         n_qubit: int,
         bias: np.ndarray,
         angles: np.ndarray,
+        measurement_twirl: bool = False,
         measurement_layout: list[int] | None = None,  # TODO: add | Layout
         shot_batch_size: int = 1,
         seed_rng: int | Generator | None = None,
@@ -57,6 +58,10 @@ class RandomizedProjectiveMeasurements(POVMImplementation[RPMMetadata]):
                 different angles of each effect. I.e. its length equals two times the number of
                 PVMs (since we have 2 angles per PVMs). If 2D, it will have a new set of angles
                 for each qubit.
+            measurement_twirl : option to randomly twirl the measurements. For each single-qubit
+                projective measurement, random twirling is equivalent to randomly flipping the
+                measurement. This is equivalent to randomly taking the opposite Bloch vector in
+                the Bloch sphere representation.
             measurement_layout: list of indices specifying on which qubits the POVM
                 acts. If None, two cases can be distinguished: 1) if a circuit supplied
                 to the :meth:`.compose_circuits` has been transpiled, its final
@@ -110,6 +115,7 @@ class RandomizedProjectiveMeasurements(POVMImplementation[RPMMetadata]):
                 f" ``n_qubit`` ({n_qubit})."
             )
         self.angles = angles.reshape((self.n_qubit, self._n_PVMs, 2))
+        self.measurement_twirl = measurement_twirl
 
         self.msmt_qc = self._build_qc()
 
@@ -251,6 +257,25 @@ class RandomizedProjectiveMeasurements(POVMImplementation[RPMMetadata]):
             ).reshape(  # Reshape to match the shape of ``pvm_idx``.
                 (*circuit_binding.shape, num_batches)
             )
+            # If the twirling option is turned on we double the number of PVMs
+            # because each PVM can be twirled. The encoding works as follows :
+            #   `pvm_idx % self._n_PVMS` is the index of the PVM used and
+            #   `pvm_idx // self._n_PVMS` indicates if the PVM has been twirled.
+            # For the example of :class:`ClassicalShadows`, we have:
+            #   n_PVMs = 3
+            #   pvm_idx == 0 -> untwirled Z-measurement: {|0><0|, |1><1|}
+            #   pvm_idx == 1 -> untwirled X-measurement: {|+><+|, |-><-|}
+            #   pvm_idx == 2 -> untwirled Y-measurement: {|+i><+i|, |-i><-i|}
+            #   pvm_idx == 3 -> twirled Z-measurement: {|1><1|, |0><0|}
+            #   pvm_idx == 4 -> twirled X-measurement: {|-><-|, |+><+|}
+            #   pvm_idx == 5 -> twirled Y-measurement: {|-i><-i|, |+i><+i|}
+            if self.measurement_twirl:
+                pvm_idx[..., i] += self._n_PVMs * self._rng.integers(
+                    2,
+                    size=circuit_binding.size * num_batches,
+                ).reshape(  # Reshape to match the shape of ``pvm_idx``.
+                    (*circuit_binding.shape, num_batches)
+                )
 
         # Transform the PVM indices into actual measurement parameters that are
         # then coerced into a :class:``BindingsArray`` object.
@@ -358,10 +383,13 @@ class RandomizedProjectiveMeasurements(POVMImplementation[RPMMetadata]):
         theta: np.ndarray = np.empty(pvm_idx.shape)
         phi: np.ndarray = np.empty(pvm_idx.shape)
         for multi_index in np.ndindex(pvm_idx.shape):
-            # multi_index.shape = (*pv, num_batches, n_qubit)
+            # multi_index.shape is (*pv, num_batches, n_qubit)
             i_qubit = multi_index[-1]
-            theta[multi_index] = self.angles[i_qubit, pvm_idx[multi_index], 0]
-            phi[multi_index] = self.angles[i_qubit, pvm_idx[multi_index], 1]
+            actual_pvm_idx = pvm_idx[multi_index] % self._n_PVMs
+            twirl = pvm_idx[multi_index] // self._n_PVMs
+
+            theta[multi_index] = self.angles[i_qubit, actual_pvm_idx, 0] + np.pi * (twirl > 0)
+            phi[multi_index] = self.angles[i_qubit, actual_pvm_idx, 1]
 
         return BindingsArray(data={self._qc_theta: theta, self._qc_phi: phi})
 
@@ -372,6 +400,11 @@ class RandomizedProjectiveMeasurements(POVMImplementation[RPMMetadata]):
     ) -> tuple[int, ...]:
         """Transform a PVM index and a bitstring outcome to a POVM outcome.
 
+        The method takes into account the possible twirling of the measurements
+        and un-do its effect. For single-qubit projective measurements, it means
+        to perform a bit-flip of the classical outcome for each twirled projective
+        measurement.
+
         Args:
             pvm_idx: qubit-wise index indicating which PVM was used to perform the measurement.
             bitstring_outcomes: the outcome of the measurements performed with the PVM label
@@ -381,7 +414,10 @@ class RandomizedProjectiveMeasurements(POVMImplementation[RPMMetadata]):
             A tuple of indices indicating the POVM outcomes on each qubit. For each qubit,
             the index goes from :math:``0`` to :math:``2 * self.n_PVM - 1``.
         """
-        return tuple(pvm_idx[i] * 2 + int(bit) for i, bit in enumerate(bitstring_outcome[::-1]))
+        return tuple(
+            (pvm_idx[i] % self._n_PVMs) * 2 + (int(bit) + pvm_idx[i] // self._n_PVMs) % 2
+            for i, bit in enumerate(bitstring_outcome[::-1])
+        )
 
     def definition(self) -> ProductPOVM:
         """Return the POVM corresponding to this implementation."""
diff --git a/pyproject.toml b/pyproject.toml
index ec0831d..c3cfb76 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -132,7 +132,7 @@ ignore = [
 explicit-preview-rules = true
 
 [tool.ruff.lint.pylint]
-max-args = 7
+max-args = 10
 
 [tool.ruff.lint.extend-per-file-ignores]
 "docs/**/*" = [


Test patch:
diff --git a/test/library/test_randomized_projective_measurements.py b/test/library/test_randomized_projective_measurements.py
index 6703ab4..3778fe6 100644
--- a/test/library/test_randomized_projective_measurements.py
+++ b/test/library/test_randomized_projective_measurements.py
@@ -13,7 +13,13 @@
 from unittest import TestCase
 
 import numpy as np
-from povm_toolbox.library import RandomizedProjectiveMeasurements
+from numpy.random import default_rng
+from povm_toolbox.library import ClassicalShadows, RandomizedProjectiveMeasurements
+from povm_toolbox.post_processor import POVMPostProcessor
+from povm_toolbox.sampler import POVMSampler
+from qiskit import QuantumCircuit
+from qiskit.primitives import StatevectorSampler
+from qiskit.quantum_info import SparsePauliOp
 
 
 class TestRandomizedPMs(TestCase):
@@ -51,6 +57,39 @@ class TestRandomizedPMs(TestCase):
 
             self.assertEqual(qc.num_qubits, n_qubit)
 
+    def test_twirling(self):
+        """Test if the twirling option works correctly."""
+        rng = default_rng(13)
+
+        qc = QuantumCircuit(2)
+        qc.h(0)
+
+        n_qubit = qc.num_qubits
+        measurement = ClassicalShadows(n_qubit, seed_rng=rng, measurement_twirl=True)
+
+        rng2 = default_rng(26)
+
+        sampler = StatevectorSampler(seed=rng2)
+        povm_sampler = POVMSampler(sampler=sampler)
+
+        job = povm_sampler.run([qc], shots=128, povm=measurement)
+        pub_result = job.result()[0]
+
+        post_processor = POVMPostProcessor(pub_result)
+
+        observable = SparsePauliOp(["ZI"], coeffs=[1.0])
+        exp_value, _ = post_processor.get_single_exp_value_and_std(observable)
+        self.assertAlmostEqual(exp_value, 0.9843749999999998)
+        observable = SparsePauliOp(["IZ"], coeffs=[1.0])
+        exp_value, _ = post_processor.get_single_exp_value_and_std(observable)
+        self.assertAlmostEqual(exp_value, 0.07031249999999983)
+        observable = SparsePauliOp(["ZY"], coeffs=[1.0])
+        exp_value, _ = post_processor.get_single_exp_value_and_std(observable)
+        self.assertAlmostEqual(exp_value, 0.0)
+        observable = SparsePauliOp(["IX"], coeffs=[1.0])
+        exp_value, _ = post_processor.get_single_exp_value_and_std(observable)
+        self.assertAlmostEqual(exp_value, 1.1718749999999998)
+
     # TODO: write a unittest for each public method of RandomizedProjectiveMeasurements
 
     # TODO: write a unittest to assert the correct handling of invalid inputs (i.e. verify that