# How to get the density matrix from a stabilizer table in qiskit

I am new to qiskit and quantum computing in general, so bear with me please. For my bachelor's thesis, I am programming qiskit to first generate a random Clifford circuit (qc) and to then measure the output state once.

For that I am using

$$\texttt{measurement = StabilizerState(qc).measure()}$$

which returns a tuple, containing the measurement outcome and its stabilizer table, For example

$$\texttt{(011, StabilizerTable: ['ZII', '-IZI', '-IIZ']).}$$

Now I would like to apply the inverse Clifford circuit to the outcome state (thereby transforming the stabilizer table to [U†$$\cdot$$ZII$$\cdot$$U, -U†$$\cdot$$IZI$$\cdot$$U, -U†$$\cdot$$IIZ$$\cdot$$U], where U† is the global unitary representing the inverse Clifford circuit.) For that I am using

$$\texttt{measurement.evolve(StabilizerState(qc.inverse())}$$

which returns a new stabilizer table describing the transformed stabilizer state. I would now like to know some sort of representation of the new state, like its density matrix.

My question is: is there some way to get the underlying state vector or density matrix from a stabilizer table? My first idea, reading through the documentation was

$$\texttt{StabilizerState(input).to_operator}$$

which did not work since it only returns the input in form of an operator (for example qc).

I guess I could always just convert the first measurement result (011) into an array and calculate U† |011> <011| U the old fashioned way, which doesn't seem very efficient to me though.

EDIT: you can build the desired DensityMatrix starting from an empty $$n$$-qubits circuit (in the state $$| \psi \rangle = | 0 \rangle^{\otimes n}$$) and then use the DensityMatrix.evolve method passing your StabilizerState object:

from qiskit import QuantumCircuit
from qiskit.quantum_info import DensityMatrix, StabilizerState
from qiskit.visualization import array_to_latex

qc = QuantumCircuit(2)
rho = DensityMatrix(qc)
qc.h(0)
qc.x(1)
stabstate = StabilizerState(qc)
rho = rho.evolve(stabstate)

array_to_latex(rho)


$$\begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & \tfrac{1}{2} & \tfrac{1}{2} \\ 0 & 0 & \tfrac{1}{2} & \tfrac{1}{2} \\ \end{bmatrix}$$

• first of all, thank you for your response ! :) but sadly, this still doesn't work (maybe I'm doing it wrong). Suppose i have a quantum circuit qc on 2 qubits. on q0 it applys an H gate and on qb1 an X Gate. The output state would be 1/sqrt 2 10 + 11. The function StabilizerState(qc) confirms this with its output: StabilizerTable: ['+IX', '-ZI']. But if i then go on with your suggestion: rho = DensityMatrix(StabilizerState(qc)), i dont get the desired result of XH |00><00|HX but instead it gives me the matrix representation of HX Jan 30 at 17:48
• edit: DensityMatrix(qc) does in fact yield the density matrix of XH |00>. So how can i make it work for the StabilizerState too ? :) Jan 30 at 19:01
• Ok now I got what you mean. I edited my answer to show the correct solution applied to the specific example you did in the comments Jan 31 at 10:14
• thank you very much, that was helpful ! :) Jan 31 at 20:52
• i do have one more question, hopefully this thread is not too old already: i would like to calculate <b| U rho U † |b>, where b is some standard basis state, U is a clifford operation and rho is a stabilizer state. Is the method, you described above the most efficient one (in the sense of Gottesmann and Knills theorem) to get the evolved density matrix, or can it be improved ? Do i have to calculate <b| U rho U † |b> 'by hand', i.e. standard matrix multiplication (with numpy etc), or is there some built in qiskit function for that kind of scenario ? Feb 7 at 13:33

If your input state is $$|000\rangle$$, then you can easily get the output as a Statevector as follows:

sv = Statevector.from_label('000').evolve(stab)
sv.draw('latex')


Similarly, you can get the output as a DensityMatrix:

dm = DensityMatrix.from_label('000').evolve(stab)
dm.draw('latex')

• thanks to you too, i appreciate your help ! :) Jan 31 at 20:52