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

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

In your case the partial trace of $M=|pranglelangle p|$ over $B$ can be computed as $$ text{Tr}_B(M) = text{Tr}_Bbigg(sum_{k,l} M_{k,l} otimes |kranglelangle 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 $|0000rangle, |0001rangle, ... , |1111rangle$ – 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 $$ |psirangle=sum_{i,j,k,l}c_{ijkl}|ijrangle_A|klrangle_B $$ Rewrite this as a matrix $$ C=sum_{i,j,k,l}c_{ijkl}|ijranglelangle kl| $$ Effectively, you just have to reshape your numpy array.

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