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.