# How to find the reduced density matrix of a four-qubit system?

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.

• No other imports? Ok to import something that already has partial traces? Aug 26 '19 at 23:04
• 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)
– glS
May 17 at 10:25

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

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.