# Stabilizer State - efficient calculation of measurement probabilities - Qiskit

I would like to calculate the probability of measuring some state $$U\rho U^\dagger$$ in the basis state $$b \in (0,1)^{\otimes n}$$, i.e. $$$$. Now, according to Gottesmann and Knill’s theorem, this can be efficiently calculated, if $$\rho$$ is a stabilizer state and U is some Clifford operation.

I made several attempts, figuring this out by myself. In order to achieve at least some result, I resorted to calculate the problem the old fashioned way, matrix multiplication style. But as one expects, the bigger the system size, the longer it takes to calculate. Also i am planning to calculate up to 5000 probabilities, each for different basis states and clifford operations.

In qiskit I tried the following:

from qiskit import QuantumCircuit
from qiskit.quantum_info import DensityMatrix, random_clifford, StabilizerState
from numpy import zeros

qc = ghz(3)
rnd_clifford = random_clifford(3)                    #generates a random clifford operation
rho_ghz = DensityMatrix(qc)
rho_ghz.evolve(rnd_clifford)                         #evolves the GHZ state with random clifford operation
#this should be equivalent to U*rho*U†

prob_b = basis(3, b).transpose()@rho_ghz@basis(3, b)  #calculation of <b|U*rho*U†|b>

def ghz(num_qubits):
qc = QuantumCircuit(num_qubits)
qc.h(0)
for k in range(1, num_qubits, 1):
qc.cx(0,k)
qc.barrier()
return qc

def basis(num_qubits, bitstring):
bas = zeros(2**num_qubits); bas[int(bitstring, 2)] = 1
return bas


I want to iterate this procedure many times and save the results to a list. Especially for many iterations and high dimensionality (exponential in qubits) of the GHZ state, this can take a long while to calculate.

So my question is: Is there a more efficient way to calculate $$$$, given the prerequisites of $$\rho$$ being a stabilizer state and U being some Clifford Operator, i.e. making use of Gottesmann and Knills theorem?

Anoher idea i had, would be to maybe make use of the $$\texttt{qiskit.quantum_info.StabilizerState.probabilities_dict}$$ method, and to then get the dictionary entry for the desired basis b

• Why don't you just use the qasm simulator and take a lot of shots? And use those to estimate the frequencies? If you go down this route, it's better to use Stim, as it will yield even faster simulations. Feb 7 at 20:55
• it's important that this step is in post processing, i.e. i dont want to retrieve the actual probability distribution of the evolved state, but rather calculate a number of single matrix elements <b|UrhoU†|b>, each for different matrices UrhoU† :) Feb 7 at 22:51
• Scott Aaronson's CHP performs / shows how to perform measurements on stabilizer states in $O(n^2)$ time. Feb 8 at 3:42
• You should edit your question to state that you are interested in "probability amplitudes", rather than "measurement probabilities". Some questions (1) $\rho$ drawn from a small set or can be an arbitrary stabilizer state? (2) Is $\rho$ a pure stabilizer state or mixture of them? (3) What is the largest $\rho$ you are going to be working with? Feb 8 at 21:41
• you are right, i should have specified - rho is a pure, n-qubit ghz state, calculations are planned for n in [3, 10] qubits Feb 8 at 22:59

You have indicated that $$\rho$$ is a pure $$n$$-qubit GHZ state, $$|G_n\rangle = \frac{|0\rangle^{\otimes n} + |1\rangle^{\otimes n}}{\sqrt{2}}.$$ Then, $$\langle b| U\rho U^\dagger |b\rangle = \langle b| U |G_n\rangle\langle G_n| U^\dagger |b\rangle = |\langle b| U |G_n\rangle|^2.$$
You have further indicated that $$|b\rangle$$ is the standard computational basis, so there is a 1 in exactly one spot of the vector $$|b\rangle$$. Furthermore, $$|G_n\rangle$$ only has non-zero entries in the first and last spot. So, $$U|G_n\rangle$$ is just the sum of the first and last columns of $$U$$. And $$\langle b| U |G_n\rangle$$ is the $$b$$-th entry of $$U|G_n\rangle$$.
In other words, $$\langle b| U |G_n\rangle$$ is the sum of the matrix elements $$U_{b,0} + U_{b,2^n-1}$$.
So, once you generate your $$U$$, only a trivial computation is needed to compute $$\langle b| U\rho U^\dagger |b\rangle$$.