Stabilizer States - Calculating measurement probabilities with the rank of the stabilizer table's X-block

Question

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)
    qc.h(0)
    for k in range(1, num_qubits, 1):
         qc.cx(0,k)
    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:
        prob.append(prob_dict[key])

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 ?

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_)

But...

assert np.linalg.matrix_rank(m) == 2
# fails because np.linalg.matrix_rank(m) returns 3