Repository: PennyLaneAI/pennylane
Pull request #4959 (2023-12-18 18:06:20)

Improve entropy estimates from classical shadow density matrix estimators
As [report by a user in the form](https://discuss.pennylane.ai/t/different-outputs-of-two-methods-for-the-same-entropy/3800/3) we get different behaviors for `qml.math.vn_entropy` and `ClassicalShadow.entropy`.

Our previous solution was not optimal. We now project to the closest valid density matrix via the algorithm outlined below eq. (16) in https://arxiv.org/pdf/1106.5458.pdf

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
Qottmann: There is still an underlying issue of when reconstructed states of classical shadows can be negative and for that case what is the _correct_ way to treat it (if there is any)
josh146: Thanks @Qottmann! Don't forget to update the changelog.

I don't think this is the case, but would this be considered a breaking change? E.g., would any of the numerics of the classical shadow entropy change?
Qottmann: > I don't think this is the case, but would this be considered a breaking change? 
No I wouldnt say, though there is now a deprecation (I started a cycle)
> E.g., would any of the numerics of the classical shadow entropy change?
If anything, they should all improve i.e. give better estimates when comparing with exact state vector calculations 🙂
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/4959?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Comparison is base [(`f496a63`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/f496a638f74a87a8b941e4b9cdbfaa55b2e6ea5c?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.50% compared to head [(`8b8f089`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/4959?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.49%.
> Report is 1 commits behind head on master.


<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #4959      +/-   ##
==========================================
- Coverage   99.50%   99.49%   -0.01%     
==========================================
  Files         390      390              
  Lines       35496    35295     -201     
==========================================
- Hits        35320    35118     -202     
- Misses        176      177       +1     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/4959?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

Qottmann: sc-52334-remove-atol-keyword-in-classicalshadow-entropy
Qottmann: > Since this fix came from a user, can we turn that comparison into a test that would have caught the problematic behaviour?

I was thinking about this but left it with just a check that the entropies from classical shadows are computed vor valid density matrices (in `test_closest_density_matrix()`). For the case of using `qml.math.vn_entropy()`, there is already the option to `check_state=True`, which should be done on states reconstructed with classical shadows.
If you have a good idea for an additional test I'm happy to implement it though :)
Qottmann: amazing, thanks for spotting all these little quirks @lillian542 ! 💯 nicer now!

Patch:
diff --git a/doc/development/deprecations.rst b/doc/development/deprecations.rst
index 85931939a..c96e68538 100644
--- a/doc/development/deprecations.rst
+++ b/doc/development/deprecations.rst
@@ -93,6 +93,13 @@ Pending deprecations
         some_qfunc(params)
         return qml.expval(Hamiltonian)
 
+* ``ClassicalShadow.entropy()`` no longer needs an ``atol`` keyword as a better 
+  method to estimate entropies from approximate density matrix reconstructions
+  (with potentially negative eigenvalues) has been implemented.
+
+  - Deprecated in v0.34
+  - Will be removed in v0.35
+
 Completed deprecation cycles
 ----------------------------
 
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index eb143426c..fc8eec90f 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -290,6 +290,11 @@
 
 <h3>Improvements 🛠</h3>
 
+* `ClassicalShadow.entropy` now uses the algorithm outlined in 
+  [1106.5458](https://arxiv.org/abs/1106.5458) to project the approximate density matrix
+  (with potentially negative eigenvalues) onto the closest valid density matrix.
+  [(#4959)](https://github.com/PennyLaneAI/pennylane/pull/4959)
+
 <h4>Community contributions 🥳</h4>
 
 * The `+=` operand can now be used with a `PauliSentence`, which has also provided
@@ -528,6 +533,10 @@
   Users should now validate these properties manually.
   [(#4773)](https://github.com/PennyLaneAI/pennylane/pull/4773)
 
+* With an algorithmic improvement to `ClassicalShadow.entropy`, the keyword `atol`
+  becomes obsolete and will be removed in v0.35.
+  [(#4959)](https://github.com/PennyLaneAI/pennylane/pull/4959)
+
 <h3>Documentation 📝</h3>
 
 * Documentation for unitaries and operations decompositions was moved from `qml.transforms` to `qml.ops.ops_math`.
@@ -661,6 +670,7 @@ Josh Izaac,
 Juan Giraldo,
 Emiliano Godinez Ramirez,
 Ankit Khandelwal,
+Korbinian Kottmann,
 Christina Lee,
 Vincent Michaud-Rioux,
 Anurav Modak,
diff --git a/pennylane/shadows/classical_shadow.py b/pennylane/shadows/classical_shadow.py
index 55c180f17..8807b6bec 100644
--- a/pennylane/shadows/classical_shadow.py
+++ b/pennylane/shadows/classical_shadow.py
@@ -324,7 +324,7 @@ class ClassicalShadow:
 
         return qml.math.squeeze(results)
 
-    def entropy(self, wires, snapshots=None, alpha=2, k=1, base=None, atol=1e-5):
+    def entropy(self, wires, snapshots=None, alpha=2, k=1, base=None, atol=None):
         r"""Compute entropies from classical shadow measurements.
 
         Compute general Renyi entropies of order :math:`\alpha` for a reduced density matrix :math:`\rho` in terms of
@@ -337,7 +337,10 @@ class ClassicalShadow:
 
         In the case of :math:`\alpha = 2`, the Renyi entropy becomes the logarithm of the purity of the reduced state
 
-        .. math:: S_{\alpha=2}(\rho) = - \log\left(\text{tr}(\rho^2) \right)
+        .. math:: S_{\alpha=2}(\rho) = - \log\left(\text{tr}(\rho^2) \right).
+
+        Since density matrices reconstructed from classical shadows can have negative eigenvalues, we use the algorithm described in
+        `1106.5458 <https://arxiv.org/abs/1106.5458>`_ to project the estimator to the closest valid state.
 
         .. warning::
 
@@ -356,7 +359,6 @@ class ClassicalShadow:
             k (int): Allow to split the snapshots into ``k`` equal parts and estimate the snapshots in a median of means fashion. There is no known advantage to do this for entropies.
                 Thus, ``k=1`` is default and advised.
             base (float): Base to the logarithm used for the entropies.
-            atol (float): Absolute tolerance for eigenvalues close to 0 that are taken into account.
 
         Returns:
             float: Entropy of the chosen subsystem.
@@ -381,7 +383,7 @@ class ClassicalShadow:
             bits, recipes = max_entangled_circuit()
             shadow = qml.ClassicalShadow(bits, recipes)
 
-            entropies = [shadow.entropy(wires=[0], alpha=alpha, atol=1e-2) for alpha in [1., 2., 3.]]
+            entropies = [shadow.entropy(wires=[0], alpha=alpha) for alpha in [1., 2., 3.]]
 
         >>> np.isclose(entropies, entropies[0], atol=1e-2)
         [ True,  True,  True]
@@ -404,42 +406,47 @@ class ClassicalShadow:
             bitstrings, recipes = qnode(x)
             shadow = qml.ClassicalShadow(bitstrings, recipes)
 
-        >>> [shadow.entropy(wires=wires, alpha=alpha, atol=1e-10) for alpha in [1., 2., 3.]]
+        >>> [shadow.entropy(wires=wires, alpha=alpha) for alpha in [1., 2., 3.]]
         [1.5419292874423107, 1.1537924276625828, 0.9593638767763727]
 
         """
 
+        if atol is not None:
+            warnings.warn(
+                "Trying to use atol in ClassicalShadow.entropy. This feature is deprecated and has been replaced by a more accurate version that projects the density matrix estimator to the closest valid density matrix in L2 norm. You can simply remove the atol keyword in your code.",
+                qml.PennyLaneDeprecationWarning,
+            )
+
         global_snapshots = self.global_snapshots(wires=wires, snapshots=snapshots)
         rdm = median_of_means(global_snapshots, k, axis=0)
 
         # Allow for different log base
         div = np.log(base) if base else 1
 
-        # This was returning negative values in most cases, so commenting it out
-        # until we figure out the issue.
-        # if alpha == 2:
-        #     # special case of purity
-        #     res = -qml.math.log(qml.math.trace(rdm @ rdm))
-        #     return res / div
-
-        # Else
-        # Compute Eigenvalues and choose only those >>0
-        evs = qml.math.eigvalsh(rdm)
-        mask0 = qml.math.logical_not(qml.math.isclose(evs, 0, atol=atol))
-        mask1 = qml.math.where(evs > 0, True, False)
-        mask = qml.math.logical_and(mask0, mask1)
-        # Renormalize because of cropped evs
-        evs_nonzero = qml.math.gather(evs, mask)
-        evs_nonzero = evs_nonzero / qml.math.sum(evs_nonzero)
-
+        evs_nonzero = _project_density_matrix_spectrum(rdm)
         if alpha == 1:
             # Special case of von Neumann entropy
-            return -qml.math.sum(evs_nonzero * qml.math.log(evs_nonzero)) / div
+            return qml.math.entr(evs_nonzero) / div
 
         # General Renyi-alpha entropy
         return qml.math.log(qml.math.sum(evs_nonzero**alpha)) / (1.0 - alpha) / div
 
 
+def _project_density_matrix_spectrum(rdm):
+    """Project the estimator density matrix rdm with possibly negative eigenvalues onto the closest true density matrix in L2 norm"""
+    # algorithm below eq. (16) in https://arxiv.org/pdf/1106.5458.pdf
+    evs = qml.math.eigvalsh(rdm)[::-1]  # order from largest to smallest
+    d = len(rdm)
+    a = 0.0
+    for i in range(d - 1, -1, -1):
+        if evs[i] + a / (i + 1) > 0:
+            break
+        a += evs[i]
+
+    lambdas = evs[: i + 1] + a / (i + 1)
+    return lambdas[::-1]
+
+
 # Util functions
 def median_of_means(arr, num_batches, axis=0):
     r"""


Test patch:
diff --git a/tests/shadow/test_shadow_entropies.py b/tests/shadow/test_shadow_entropies.py
index 98d611b41..e969d0f4a 100644
--- a/tests/shadow/test_shadow_entropies.py
+++ b/tests/shadow/test_shadow_entropies.py
@@ -20,6 +20,7 @@ import pytest
 import pennylane as qml
 import pennylane.numpy as np
 from pennylane.shadows import ClassicalShadow
+from pennylane.shadows.classical_shadow import _project_density_matrix_spectrum
 
 np.random.seed(777)
 
@@ -53,6 +54,14 @@ def expected_entropy_ising_xx(param, alpha):
     return qml.math.log(qml.math.sum(eigs**alpha)) / (1 - alpha)
 
 
+def test_deprecation_warning_atol():
+    """Test that a depceration warning is raised when trying to use atol"""
+    bits, recipes = max_entangled_circuit(wires=2)()
+    shadow = ClassicalShadow(bits, recipes)
+    with pytest.warns(qml.PennyLaneDeprecationWarning):
+        shadow.entropy([0, 1], atol=1e-3)
+
+
 class TestShadowEntropies:
     """Tests for entropies in ClassicalShadow class"""
 
@@ -64,9 +73,7 @@ class TestShadowEntropies:
         bits, recipes = max_entangled_circuit(wires=n_wires)()
         shadow = ClassicalShadow(bits, recipes)
 
-        entropies = [
-            shadow.entropy(wires=[0], alpha=alpha, atol=1e-2, base=base) for alpha in [1, 2, 3]
-        ]
+        entropies = [shadow.entropy(wires=[0], alpha=alpha, base=base) for alpha in [1, 2, 3]]
         assert np.allclose(entropies, entropies[0], atol=1e-2)
         expected = np.log(2) / np.log(base)
         assert np.allclose(entropies, expected, atol=2e-2)
@@ -125,9 +132,7 @@ class TestShadowEntropies:
             exact = [exact_vn, exact_alpha, exact_2]
 
             # shadow estimate
-            entropies = [
-                shadow.entropy(wires=rdm_wires, alpha=alpha, atol=1e-10) for alpha in [1, 1.5, 2]
-            ]
+            entropies = [shadow.entropy(wires=rdm_wires, alpha=alpha) for alpha in [1, 1.5, 2]]
 
             assert np.allclose(entropies, exact, atol=1e-1)
 
@@ -155,3 +160,40 @@ class TestShadowEntropies:
                     expected_entropy = expected_entropy_ising_xx(param, alpha) / np.log(base)
 
                     assert qml.math.allclose(entropy, expected_entropy, atol=1e-1)
+
+    def test_closest_density_matrix(self):
+        """Test that the closest density matrix from the estimator is valid"""
+        wires = 3
+        dev = qml.device("default.qubit", wires=range(wires))
+
+        # Just a simple circuit
+        @qml.qnode(dev)
+        def qnode(x):
+            for i in range(wires):
+                qml.RY(x[i], wires=i)
+
+            for i in range(wires - 1):
+                qml.CNOT((i, i + 1))
+
+            return qml.state()
+
+        x = np.linspace(0.5, 1.5, num=wires)
+        state = qnode(x)
+        rho = qml.math.dm_from_state_vector(state)
+        lambdas = _project_density_matrix_spectrum(rho)
+        assert np.isclose(np.sum(lambdas), 1.0)
+        assert all(lambdas > 0)
+
+
+rho0 = np.zeros((2**3, 2**3))
+rho0[0, 0] = 1.0
+
+rho1 = np.diag([-0.1, -0.1, -0.1, 1.3])
+rho2 = np.diag([-0.1, -0.1, 0.1, 1.1])
+
+
+@pytest.mark.parametrize("rho", [rho0, rho1, rho2])
+def test_project_density_matrix_spectrum(rho):
+    """Test the function _project_density_matrix_spectrum behaves as expected for trivial case"""
+    new_lambdas = _project_density_matrix_spectrum(rho)
+    assert qml.math.allclose(new_lambdas, [1.0])