Clip probabilities in QuantumState (Qiskit#9762)

* Clip probabilities

* clip implictly

* add reno

* add tests

* add renormalization

* use sum, not norm

* ensure round still works

* normalizing seems to re-introduce errors

* Tighten floating-point tests

In cases of rounding and clipping, it's important that the
floating-point output is bit-for-bit correct, so the fuzzy tests weren't
ideal (some of these were strict, but it was inconsistent).  It's better
to use `assertEqual` rather than `assertTrue` where possible so we get
better errors on failure.

---------

Co-authored-by: Jake Lishman <jake.lishman@ibm.com>

qiskit/quantum_info/states/densitymatrix.py

@@ -480,8 +480,13 @@ def probabilities(self, qargs=None, decimals=None):
         probs = self._subsystem_probabilities(
             np.abs(self.data.diagonal()), self._op_shape.dims_l(), qargs=qargs
         )
+
+        # to account for roundoff errors, we clip
+        probs = np.clip(probs, a_min=0, a_max=1)
+
         if decimals is not None:
             probs = probs.round(decimals=decimals)
+
         return probs
 
     def reset(self, qargs=None):

qiskit/quantum_info/states/quantum_state.py

@@ -234,7 +234,9 @@ def probabilities_dict(self, qargs=None, decimals=None):
             dict: The measurement probabilities in dict (ket) form.
         """
         return self._vector_to_dict(
-            self.probabilities(qargs=qargs, decimals=decimals), self.dims(qargs), string_labels=True
+            self.probabilities(qargs=qargs, decimals=decimals),
+            self.dims(qargs),
+            string_labels=True,
         )
 
     def sample_memory(self, shots, qargs=None):

qiskit/quantum_info/states/statevector.py

@@ -575,8 +575,13 @@ def probabilities(self, qargs=None, decimals=None):
         probs = self._subsystem_probabilities(
             np.abs(self.data) ** 2, self._op_shape.dims_l(), qargs=qargs
         )
+
+        # to account for roundoff errors, we clip
+        probs = np.clip(probs, a_min=0, a_max=1)
+
         if decimals is not None:
             probs = probs.round(decimals=decimals)
+
         return probs
 
     def reset(self, qargs=None):

releasenotes/notes/clip-quantumstate-probabilities-5c9ce05ffa699a63.yaml

@@ -0,0 +1,8 @@
+---
+fixes:
+  - |
+    Clip probabilities in the :meth:`.QuantumState.probabilities` and 
+    :meth:`.QuantumState.probabilities_dict` methods to the interval ``[0, 1]``.
+    This fixes roundoff errors where probabilities could e.g. be larger than 1, leading 
+    to errors in the shot emulation of the :class:`.Sampler`.
+    Fixed `#9761 <https://github.com/Qiskit/qiskit-terra/issues/9761>`__.

test/python/quantum_info/states/test_densitymatrix.py

@@ -1202,6 +1202,27 @@ def test_drawings(self):
             with self.subTest(msg=f"draw('{drawtype}')"):
                 dm.draw(drawtype)
 
+    def test_clip_probabilities(self):
+        """Test probabilities are clipped to [0, 1]."""
+        dm = DensityMatrix([[1.1, 0], [0, 0]])
+
+        self.assertEqual(list(dm.probabilities()), [1.0, 0.0])
+        # The "1" key should be exactly zero and therefore omitted.
+        self.assertEqual(dm.probabilities_dict(), {"0": 1.0})
+
+    def test_round_probabilities(self):
+        """Test probabilities are correctly rounded.
+
+        This is good to test to ensure clipping, renormalizing and rounding work together.
+        """
+        p = np.sqrt(1 / 3)
+        amplitudes = [p, p, p, 0]
+        dm = DensityMatrix(np.outer(amplitudes, amplitudes))
+        expected = [0.33, 0.33, 0.33, 0]
+
+        # Exact floating-point check because fixing the rounding should ensure this is exact.
+        self.assertEqual(list(dm.probabilities(decimals=2)), expected)
+
 
 if __name__ == "__main__":
     unittest.main()

test/python/quantum_info/states/test_statevector.py

@@ -1308,6 +1308,24 @@ def test_statevector_len(self):
         self.assertEqual(len(empty_vector), len(empty_sv))
         self.assertEqual(len(dummy_vector), len(sv))
 
+    def test_clip_probabilities(self):
+        """Test probabilities are clipped to [0, 1]."""
+        sv = Statevector([1.1, 0])
+
+        self.assertEqual(list(sv.probabilities()), [1.0, 0.0])
+        # The "1" key should be zero and therefore omitted.
+        self.assertEqual(sv.probabilities_dict(), {"0": 1.0})
+
+    def test_round_probabilities(self):
+        """Test probabilities are correctly rounded.
+
+        This is good to test to ensure clipping, renormalizing and rounding work together.
+        """
+        p = np.sqrt(1 / 3)
+        sv = Statevector([p, p, p, 0])
+        expected = [0.33, 0.33, 0.33, 0]
+        self.assertEqual(list(sv.probabilities(decimals=2)), expected)
+
 
 if __name__ == "__main__":
     unittest.main()