[Unitary Hack]Test and benchmarking the entropy calculation for mixed state

[hongyehu]

In the current implementation, there is an attribute called the rank of the stabilizer state state.r. If state.r=0, then the state is a pure stabilizer state, and all the stabilizers are used to stabilize the state. If state.r>0, then there are state.r null stabilizers, and the state is a mixed state.
In the PyClifford, there is a function called entropy that can calculate the subsystem entropy of the state. When the state is pure, this calculation is tested. The task for this bounty is to benchmark the correctness of the entropy calculation for mixed state.
A successful solution includes 1. testing result and 2. correct entropy implementation if needed

[obliviateandsurrender]

Hey! I spent some time checking the correctness of the entropy calculation in pyclifford, and it does seem to work correctly. Here's the benchmarking code snippet I wrote for this purpose -

from tqdm.auto import tqdm
import pennylane as qml
from time import time
import numpy as np
import pyclifford
import stim

total_qubits = [2, 10, 50, 100, 250, 500, 1000, 2500, 5000, 10000]
times = [[] for _ in range(len(total_qubits))]
prtys = [[] for _ in range(len(total_qubits))]

def pl_entropy(state, mask):
    tableau = stim.Tableau.from_stabilizers([stim.PauliString(repr(stab).strip()) for stab in state.stabilizers])
    z_stabs = qml.math.array([tableau.z_output(wire) for wire in range(state.N)])
    return qml.devices.DefaultClifford()._measure_stabilizer_entropy(z_stabs, mask, log_base=2)

def random_tableau_state(N):
    random_tableau = stim.Tableau.random(N)
    x2x, x2z, z2x, z2z, x_signs, z_signs = random_tableau.to_numpy()
    gs = np.vstack((np.vstack((z2x, z2z)).T.reshape(2 * N, N).T,
                                np.vstack((x2x, x2z)).T.reshape(2 * N, N).T)).astype(int)
    ps = 2 * np.concatenate((z_signs, x_signs)).astype(int) 
    return pyclifford.StabilizerState(gs=gs, ps=ps)
    
for ix, num_qubits in tqdm(enumerate(total_qubits)):

    py_state = random_tableau_state(num_qubits)

    for _ in range(5):
        idx = np.random.choice(num_qubits, size=num_qubits, replace=False)
        for i_, num_traced in tqdm(enumerate(total_qubits[:ix]), leave=False):
            start = time()
            res1 = py_state.entropy(idx[:num_traced])
            end = time()
            res2 = pl_entropy(py_state, idx[:num_traced])
            assert np.allclose(res1, res2)
            
            times[ix].append(end - start)
            prtys[ix].append(2 ** -res1)

and here are some benchmarking results (with average purity encountered in the traced subsystems) -

Some additional points -

1. pyclifford.random_clifford_state doesn't scale beyond 2000 qubits. So I had to rely on using stim to build a tableau and then manually convert the definitions of gs and ps.
2. For the above, I'll admit, it took me some playing around to know about the definition of Tableau being used in the pyclifford. It would be nice to document them more clearly.
3. I did run an initial check with Qutip (up to ~20 qubits). To have a more complete check, I relied on the implementation that, in reality, is based on the results in this paper.

[obliviateandsurrender]

Hey @hongyehu, any updates on this? :)

[hongyehu]

Hi @obliviateandsurrender, sorry for the slow response. It looks great! I will test it tomorrow. If it passes, looks like we have a winner for this bounty :-)

[obliviateandsurrender]

Thanks! Let me know how it goes and if I should be opening any PR for it!