# Measuring entanglement entropy using a stabilizer circuit simulator

I'm trying to simulate stabilizer circuits using the Clifford tableau formalism that lets you scale up to hundreds of qubits. What I want to do is find the entanglement entropy on by splitting my quantum state (defined as a line of $$N$$ qubits) at the end of my quantum circuit after applying gates and some measurements. From this paper, the authors calculate the entanglement entropy for a state defined by its stabilizers. In particular, I'm looking at Equations 10 (this is their "clipped gauge"), A16 and A17. If $$\mathcal{S}$$ is the set of stabilizers for the state, then the entropy is given by (Equation A16):

$$S = |A| - \log_2 |\mathcal{S}_A|,$$

where $$|A|$$ is the size of the bipartition of the quantum state and $$\mathcal{S}_A$$ is the part of $$\mathcal{S}$$ which acts on $$\bar{A}$$ with the identity.

I want to simulate my quantum circuit and calculate the entanglement entropy like they do in their paper, but I'm not sure what's the easiest way to do so. A lot of the tools for simulating stabilizer circuits aren't the most transparent to use. In particular, I'm trying to understand how to go from the tableau representation a lot of simulators output and the set of generators I need to calculate the entropy.

Is there a simple procedure to go from the tableau representation to the entropy? I'm trying to think of how to implement this in code.

For the actual simulator, I see there are a few options. I need measurements, so while Qiskit does offer Clifford simulation, I can't seem to do measurements with it. The others that offer a Python interface are:

If anyone has experience with these and can explain how to go from the tableau representation to the calculation of the entropy, that would be great, since these simulators usually seem to be for giving shots in the computational basis.

I think stim is the right tool for the job here, because it gives you access to the stabilizer generators and also it defines stim.PauliString which you can use to represent the stabilizers and more easily implement the manipulations being described in the paper.

Ultimately, it will be up to you to translate what the paper wants into a specific combination of methods that implements the logic they describe. Your question isn't really specific enough about which method you want to use from the paper (there's a lot of details hiding behind that equation). Feel free to ask other more specific questions about how to achieve each task in stim. In this answer I'll just be showing how you get the generators.

In stim, when using stim.TableauSimulator, you can get a current list of stabilizer generators like this:

tableau = simulator.current_inverse_tableau() ** -1
return [tableau.z_output(k) for k in range(len(tableau))]


So, basically, you'll have some method sampling stabilizer generators by running a circuit and extracting the stabilizer generators similar to this:

from typing import List, Collection
import stim
import numpy as np

def sample_stabilizers() -> List[stim.PauliString]:
# Simulate some operations. Whatever you want.
s = stim.TableauSimulator()
s.h(0)
s.cnot(0, 1)
s.measure(0)
s.do(stim.Circuit("""
H 3
CZ 1 3
"""))
...

# Get stabilizers.
tableau: stim.Tableau = s.current_inverse_tableau() ** -1
n = len(tableau)
return [tableau.z_output(k) for k in range(n)]


Once you have the stabilizer generators, you can start manipulating them. For example, here's a method that tells you the index of the leftmost non-identity term that appears in a stabilizer:

def left(stabilizer: stim.PauliString):
return min((i for i, p in enumerate(stabilizer) if p), default=0)

• Thanks for the ideas. One thing I don't understand about what you wrote is that you only extract the stabilizers from the zs part of the tableau, right? Would you need to do a similar thing for the xs as well? After some communication with the authors, what I think I really need is the binary matrix corresponding to the stabilizers of the state. I can then manipulate this directly to find the entropy (I'll post in a separate answer once I get that). If I want to extract this from Stim, is there an easy way to convert the Pauli strings into their matrix form to do linear algebra on them? – Germ Mar 27 at 16:50
• @Germ There is an x_output method for getting the xs. simulator.current_inverse_tableau() is a stabilizer tableau containing X and Z parts, and you can do algebra on it. E.g. you may have noticed I raised it to -1 to compute its inverse. The tableau class is currently quite rigid, in order to guarantee its commutativity invariants can't be broken (e.g. it would never allow you to just zero out terms). So I suspect you will need to pull out the PauliString columns or rows, as I did, in order to do some of the things described by the paper. – Craig Gidney Mar 27 at 17:05
• Okay, this makes sense. I'm slowly understanding the framework you have set up. If I may ask one more question: For the zs, if I convert a PauliString into a list, the list will contain 0s and 3s, since the tableau representation has a $Z$ being $10$ in binary. Am I making a mistake if I convert this to a $1$ when trying to build the binary matrix? I just want to make sure I'm not leaving something out. – Germ Mar 27 at 17:24
• @Germ Actually, the integers correspond to Id=0=0b00, X=1=0b01, Y=2=0b10, Z=3=0b11. Also, it's not a python list it's a stim.PauliString which you can access in many ways as if it was a list but it also supports algebra like multiplying two Pauli strings. If you really need X=1 Z=2 Y=3 instead you can workaround with e.g. array = numpy.array(pauli_string, dtype=numpy.uint8); array ^= array >> 1. – Craig Gidney Mar 27 at 17:27
• Thanks for the explanation, I see this now. What I was trying to do was convert a stim.PauliString to a list object and then use this for to create the tableau matrix. I must not understand something about how the xs and zs are stored though, since a stim.PauliString in xs has a $Z$ term, which seems like it shouldn't happen. Is there somewhere I should look to understand this better? – Germ Mar 27 at 22:29

With the help of Craig Gidney and some others, I was able to pin down the procedure to calculate the entropy. Here are the steps.

# Create your circuit with a stabilizer simulator

This can be done with whatever simulator you want. As Craig mentioned in his answer, Stim is a great tool for the job. The rest of the answer in this section will assume you're using Stim, but it's not required.

Your code will look something like this:

import stim

circuit = stim.TableauSimulator()
...

# Create the tableau representation
tableau = circuit.current_inverse_tableau() ** -1
zs = [tableau.z_output(k) for k in range(len(tableau))]
zs = np.array(zs)


What you get with tableau is a set of stim.PauliString objects, which are essentially your "stabilizers" and "destabilizers", to use the language of Aaronson's paper on page 4. For the purposes of the entropy, we only care about the stabilizers, which are given by the zs object that we define here.

Essentially, a quantum circuit starts in the product state $$|0 \rangle^{\otimes N}$$, and then gates transform the state. The idea is that the stabilizer generators for this initial state is:

$$$$g_1 = Z_1 \equiv Z \otimes \mathbb{1} \otimes \mathbb{1} \otimes \ldots \otimes \mathbb{1}, \\ g_2 = Z_2 \equiv \mathbb{1} \otimes Z \otimes \mathbb{1} \otimes \ldots \otimes \mathbb{1},$$$$ and so on until the end. If we have $$N$$ qubits, there will be exactly $$8$$ stabilizers. However, it's super important to note that these are not identified with certain qubits. For example, stabilizer $$g_1$$ is not identified with the first qubit (this tripped me up for a bit, so I wanted to note it).

The way Stim stores the stabilizers is the following: $$0$$ means an identity, $$1$$ means an $$X$$ operator, $$2$$ means a $$Y$$ operator, and $$3$$ means a $$Z$$ operator. This is because of the binary notation.

To get into the actual form of the tableau matrix, we need to make a $$N \times 2N$$ matrix, with the left $$N \times N$$ block being for the $$X$$ operators and the right block for the $$Z$$ operators.

So you can just write a little function like this:

def binaryMatrix(zStabilizers):
"""
- Purpose: Construct the binary matrix representing the stabilizer states.
- Inputs:
- zStabilizers (array): The result of conjugating the Z generators on the initial state.
Outputs:
- binaryMatrix (array of size (N, 2N)): An array that describes the location of the stabilizers in the tableau representation.
"""
N = len(zStabilizers)
binaryMatrix = np.zeros((N,2*N))
r = 0 # Row number
for row in zStabilizers:
c = 0 # Column number
for i in row:
if i == 3: # Pauli Z
binaryMatrix[r,N + c] = 1
if i == 2: # Pauli Y
binaryMatrix[r,N + c] = 1
binaryMatrix[r,c] = 1
if i == 1: # Pauli X
binaryMatrix[r,c] = 1
c += 1
r += 1

return binaryMatrix


Now, we're ready to calculate the entropy.

# Calculating the entropy of a cut

Now that we have the matrix corresponding to the evolved quantum state through the circuit, we want to find the entanglement entropy. In this paper, the key equation is Equation A19, but the real helpful comment I found on this was on Footnote 11 of this paper (page 10), which says how to do this numerically. I also communicated with one of the authors (Xiao Chen) from the other paper here, and his comments were also quite helpful.

Our system involves $$N$$ qubits, and now we want to find the entanglement entropy between two subsystems, which we will label $$A$$ and $$B$$. Equation A19 of the paper I referenced above tells us:

$$$$\label{eq:Entropy} S_A = \text{rank}\left( \text{projection}_A \mathcal{S} \right) - N_A.$$$$

In this equation, $$N_A$$ is the number of qubits in part $$A$$, $$\mathcal{S}$$ is the set of stabilizers (our binary matrix we got above), and the projection operator means we want to "truncate" the matrix so that it's only describing the parts on $$A$$.

To do this, remember that our matrix starts off as $$N \times 2N$$. We now want to truncate it so that we don't care about the qubits in region $$B$$. Mathematically, this is what the projection operator does. It pretends everything in region $$B$$ is the identity, which in our matrix means the entries become zero.

But a simpler way to deal with this is to just delete the columns needed to describe region $$B$$, since they won't play a role anyway. Numerically, the following function does the trick:

def getCutStabilizers(binaryMatrix, cut):
"""
- Purpose: Return only the part of the binary matrix that corresponds to the qubits we want to consider for a bipartition.
- Inputs:
- binaryMatrix (array of size (N, 2N)): The binary matrix for the stabilizer generators.
- cut (integer): Location for the cut.
- Outputs:
- cutMatrix (array of size (N, 2cut)): The binary matrix for the cut on the left.
"""
N = len(binaryMatrix)
cutMatrix = np.zeros((N,2*cut))

cutMatrix[:,:cut] = binaryMatrix[:,:cut]
cutMatrix[:,cut:] = binaryMatrix[:,N:N+cut]

return cutMatrix


This truncates our original $$N \times 2N$$ matrix into a $$N \times 2N_A$$ matrix, with everything else deleted.

Now, the equation for the entropy simply requires us to compute the rank of this matrix, and subtract off the number of qubits in region $$A$$. Numerically, you can do this using Gaussian elimination with Boolean variables (in others worlds, modulo 2 arithmetic), but you can also use just a plain old SVD over real variables to get the rank.

Edit: I made a mistake in saying that a regular SVD calculation will work. For reference, I was using the matrix_rank function from NumPy. After comparing it with another approach, it seems like it mostly works, but is sometimes off. As such, I'd recommend doing Gaussian elimination, with something like this.

After that, you should be good to go. The entropy is essentially a matrix rank computation. Also note that you can use whatever stabilizer circuit simulator you like, as long as in the end you get out the stabilizer generators which you can then build your binary matrix from.