I have the state vector $|p\rangle$ made up of 4 qubits. Say system A is made up of the first and second qubits while system B is made up of qubits 3 and 4. I want to find the reduced density matrix of system A.

I know I could separately extract qubits 1,2 and 3,4 into their own state vectors then find their density matrices and compute the reduced density matrix for system A.

I want to figure out how to do this without having to extract and separate the systems. First I would find the density matrix of $|p\rangle$ and then do a partial trace with respect to system B. I am not sure how to do the partial trace of system B since the system contains 2 qubits.

Can anyone help me figure this out? I am using Python and NumPy for reference.

  • $\begingroup$ No other imports? Ok to import something that already has partial traces? $\endgroup$
    – AHusain
    Aug 26, 2019 at 23:04
  • $\begingroup$ regarding the implementation with python, note that qutip already handles this with ptrace, and if you only want to use numpy the easiest (and general) way is probably via einsum (if you're ok reshaping the array into a higher-dimensional tensor) $\endgroup$
    – glS
    May 17, 2021 at 10:25

2 Answers 2


Wikipedia entry has a nice description of the partial trace and how to compute it.

In your case the partial trace of $M=|p\rangle\langle p|$ over $B$ can be computed as $$ \text{Tr}_B(M) = \text{Tr}_B\bigg(\sum_{k,l} M_{k,l} \otimes |k\rangle\langle l| \bigg) = \sum_{k} M_{k,k} $$
where $k,l \in \{00,01,10,11\}$.

Those $M_{k,l}$ are $4 \times 4$ submatrices of the $16 \times 16$ matrix $M$. What kind of submatrices depends on the actual encoding. In the big-endian encoding – where the matrix row order corresponds to $|0000\rangle, |0001\rangle, ... , |1111\rangle$ – matrix $M_{k,l}$ has element in the position $(i,j)$ that equals to an element in the position $(k,l)$ of the $(i,j)$-th block of the big matrix $M$. So the resulting sum is actually equivalent to a matrix whose elements are traces of $4 \times 4$ blocks of the matrix $M$.

Partial trace over the first subsystem A in the big-endian encoding is easier – it's just the sum of 4 blocks from the diagonal.


Computationally, the easiest way to do this is probably as follows:

Let your state be $$ |\psi\rangle=\sum_{i,j,k,l}c_{ijkl}|ij\rangle_A|kl\rangle_B $$ Rewrite this as a matrix $$ C=\sum_{i,j,k,l}c_{ijkl}|ij\rangle\langle kl| $$ Effectively, you just have to reshape your numpy array.

Then, you can calculate $$ \rho_A=CC^\dagger $$ or $$ \rho_B=C^\dagger C. $$


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.