Consider a $n$-Qubit stabilizer state $\newcommand{\ket}[1]{\vert#1\rangle}\newcommand{\bra}[1]{\langle#1\vert}\rho = \ket{\psi}\bra{\psi}$ and its $n \times 2n$ boolean stabilizer tableau.

Any Stabilizer State can be expressed as an equally weighted superposition $$ \ket{\psi} = \frac{1}{\sqrt{|A|}}\sum\limits_{x\in A}i^{f(x)}\ket{x} $$ with $f(x) \in$ {-1, 0, 1, 2}. Such that the probability of measuring any one of the basis states $\ket{x}$ is given by $p_x = \frac{1}{|A|}$. In these lecture notes (page 6) by Scott Aaronson, it is said, that the number of basis states $\ket{x}$ with nonzero amplitude is given by $$ |A| = 2^k $$ where k is the rank of the stabilizer table's X-block.

I am trying to replicate this using Qiskits StabilizerState simulator. For a given stabilizer state, I calculate its measurement probabilities via the probabilities_dict() function which seems to yield the correct result. The number of measurement outcomes matches the calculated probability.

Here's the problem: Calculating the measurement probabilities by using the rank of the stabilizer table's X-block does oftentimes not yield the same result as the aforementioned method. A short working example is provided below:

from numpy.linalg import matrix_rank as rnk
from math import log

from qiskit import QuantumCircuit
from qiskit.quantum_info import random_clifford, StabilizerState

def ghz_stab(num_qubits):
    qc = QuantumCircuit(num_qubits)
    for k in range(1, num_qubits, 1):
    return StabilizerState(qc)

qubits = 10      

rnd_clifford = random_clifford(qubits)
random_state = ghz_stab(qubits).evolve(rnd_clifford)
x_rank = rnk(random_state.clifford.stab_x)
rank_prob = 0.5**x_rank

prob_dict = random_state.probabilities_dict()
prob = []
for key in prob_dict:                  
    if prob_dict[key] not in prob:

print(f'dict prob: {prob[0]}')
print(f'log(number of non-zero amplitudes): {log(len(prob_dict), 2)}') 
print(f'X-block rank: {x_rank}')      
print(f'rank prob: {rank_prob}')
print(f'match: {rank_prob == prob[0]}')

We begin with some known state (e.g. the $n$-Qubit GHZ-state) which undergoes some random Clifford transformation. We now would like to know the basis-measurement probabilities of the new state, the actual outcomes themselves are not of interest. Instead of calculating the entire probability dictionary it seems more efficient to use the rank $k$ of the evolved stabilizer table's X-block as $p=2^{-k}$. But as previously stated, this does not seem to work reliably.

Maybe I misunderstood some part of the theory or made a mistake in the code ?


1 Answer 1


The issue is that numpy.linalg.matrix_rank is assuming you want the rank over real numbers, when actually you want the rank over integers modulo 2.

For example, here's a matrix with rank 2 over Z2:

import numpy as np

m = np.array([
    [1, 1, 0, 0],
    [1, 0, 1, 0],
    [0, 1, 1, 0],
], dtype=np.bool_)


assert np.linalg.matrix_rank(m) == 2
# fails because np.linalg.matrix_rank(m) returns 3
  • $\begingroup$ Is there already a function for calculating the rank of an array with elements from Z2 or will i just have to write some Gauss-Jordan algorithm for that myself ? $\endgroup$
    – Coryn7
    Commented May 10, 2023 at 6:35
  • $\begingroup$ @Coryn7 you probably have to write it $\endgroup$ Commented May 10, 2023 at 16:01
  • $\begingroup$ okay, thats's fine, at least i now know what to do - you have been a great help, thanks ! :) $\endgroup$
    – Coryn7
    Commented May 11, 2023 at 7:19

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.