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. $<b|U\rho U^\dagger|b>$. 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)
    for k in range(1, num_qubits, 1):
    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 $<b|U\rho U^\dagger|b>$, 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

  • $\begingroup$ 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. $\endgroup$ Feb 7, 2023 at 20:55
  • $\begingroup$ 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† :) $\endgroup$
    – Coryn7
    Feb 7, 2023 at 22:51
  • $\begingroup$ Scott Aaronson's CHP performs / shows how to perform measurements on stabilizer states in $O(n^2)$ time. $\endgroup$
    – ChrisD
    Feb 8, 2023 at 3:42
  • 1
    $\begingroup$ 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? $\endgroup$ Feb 8, 2023 at 21:41
  • $\begingroup$ you are right, i should have specified - rho is a pure, n-qubit ghz state, calculations are planned for n in [3, 10] qubits $\endgroup$
    – Coryn7
    Feb 8, 2023 at 22:59

1 Answer 1


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$.

  • $\begingroup$ thank you very very much, this has been extremely helpful (and also embarrassingly something i should have thought of myself haha) - nevertheless, this speeds the procedure up by a lot, thanks ! :) $\endgroup$
    – Coryn7
    Feb 9, 2023 at 16:55
  • $\begingroup$ No problem. Since you are an undergrad, my advice to you would be to always work out a small instance of the problem on paper before trying to write code. Much time will be saved this way. :-) $\endgroup$ Feb 9, 2023 at 20:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.