Need some help with Purity calculation

Question

maybe you could help me a little about my calculation of a quantum pure state with purification. I have this density matrix:

\begin{equation} \rho= \begin{pmatrix} 0.4489 & 0.2304 & 0.2162 & 0\\ 0.2304 & 0.2518 & 0.2399 & 0\\ 0.2162 & 0.2399 & 0.2993 & 0\\ 0 & 0 & 0 & 0 \end{pmatrix} \end{equation}

Then I calculated the Eigenvectors and Eigenvalues of this matrix:

Eigenvalues: \begin{equation} e=[0.79987375 \space\space 0.16872495 \space\space 0.03140131\space\space 0 ] \end{equation} Eigenvectors: \begin{equation} eigenvectors= \begin{pmatrix} 0.66857314 & 0.73407551 & 0.11892473 & 0\\ 0.51561267 & -0.34235582 & -0.78545278 & 0\\ 0.53586708 & -0.58645173 & 0.60738854 & 0\\ 0 & 0 & 0 & 1 \end{pmatrix} \end{equation}

After that I used the formula for purification with the computational bases: $|\Psi\rangle = \sum \sqrt{p_i} |\phi_i\rangle \otimes \lvert\psi_i\rangle$.

computational bases: \begin{equation} computational bases= \begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{pmatrix} \end{equation}

If I put this all in the equation I got this big $\Psi$

\begin{equation} \Psi= \begin{pmatrix} 0.59794281 \\ 0.21179379 \\ 0.09495786 \\ 0 \\ 0.65652529 \\ -0.14062656 \\ -0.10392167\\ 0\\ 0.10636112\\ -0.3226337\\ 0.10763176\\ 0\\ 0\\ 0\\ 0\\ 0 \end{pmatrix} \end{equation}

Then I checked if the Purifiaction is correct, so I calculated the density matrix of $\Psi$ and then I traced out B and thought the result is the density matrix $\rho$, but this is my result: \begin{equation} \rho= \begin{pmatrix} 0.41140921 & 0.35291256 & 0.00548653 & 0\\ 0.35291256 & 0.461601 & 0.10401436 & 0\\ 0.00548653 & 0.10401436 & 0.12698979 & 0\\ 0 & 0 & 0 & 1 \end{pmatrix} \end{equation}

But this isn´t equal to the density matrix $\rho$. For the calculation of the partial trace i used the qiskit lib and Python:

import qiskit.quantum_info
qubits=[0,1]
rho= qiskit.quantum_info.partial_trace(psi_density,qubits)

My question is why I didn´t get the origin rho matrix back? Did I make a mistake with the Purity calculation? I would be glad to get an explanation for this. edit: For Purifiaction is used this piece of code:

w, v = LA.eig(rho)
null_base = np.array([1,0,0,0])
one_base = np.array([0,1,0,0])
two_base = np.array([0,0,1,0])
three_base = np.array([0,0,0,1])
v_1=v[0]
v_2=v[1]
v_3=v[2]
v_4=v[3]
v_1_trans = v_1.reshape(-1,1)
v_2_trans = v_2.reshape(-1,1)
v_3_trans = v_3.reshape(-1,1)
v_4_trans = v_4.reshape(-1,1)
null_base_trans = null_base.reshape(-1,1)
one_base_trans = one_base.reshape(-1,1)
two_base_trans = two_base.reshape(-1,1)
three_base_trans = three_base.reshape(-1,1)
sum_1 = np.tensordot(v_1_trans, null_base_trans, 0) * np.sqrt(w[0])
sum_2 = np.tensordot(v_2_trans, one_base_trans, 0) * np.sqrt(w[1])
sum_3 = np.tensordot(v_3_trans, two_base_trans, 0) * np.sqrt(w[2])
sum_4 = np.tensordot(v_4_trans, three_base_trans, 0) * np.sqrt(w[3])
psi=sum_1+sum_2+sum_3+sum_4
psi=psi.reshape(16,1)

so w are the eigenvalues and v are the eigenvectors. The trans variables just only make column vectors for tensordot.

Answer 1

The fundamental problem is that you've transposed the matrix of eigenvectors compared to what you need. Yes, as you've written it, the eigenvectors are the columns of the matrix. But I think you must be using the rows in your actual code.

One of the basic checks that I did (and this is crucial whenever you code something) is to be able to do simple test cases where you can check the output. This means knowing how everything works mathematically rather than hoping that the computer will magically take care of it for you. So, what I did when checking your code was to think about the term $$ \sqrt{p_1}|\phi_1\rangle|00\rangle. $$ This is going to yield amplitudes on the components 0000,0100 and 1000 (which are binary values 0, 4, 8), which therefore correspond to the elements of the vector 1, 5 and 9. If I can't even be bothered to calculate their actual values, I at least know that, relatively, they must have the same relation. So, if the first element of $|\phi_1\rangle$ is bigger than the second, the first element of the output must be bigger than the fifth.

Eventually, I threw this into Mathematica (there are a couple of irrelevant phases that appear in this calculation compared to yours,

rho = ({{0.4489, 0.2304, 0.2162, 0},{0.2304, 0.2518, 0.2399, 0},{0.2162, 0.2399, 0.2993, 0},{0, 0, 0, 0}});
Eigenvalues[rho]
{vals, vecs} = Eigensystem[rho]
basis = IdentityMatrix[4];
Sum[KroneckerProduct[{vecs[[i]]}, {basis[[i]]}]*Sqrt[vals[[i]]], {i, 1, 4}]

This gave me the answer

{-0.597943, 0.30153, 0.021074, 0., -0.461142, -0.140627, -0.139185, 0., -0.479256, -0.240892, 0.107632, 0., 0., 0., 0., 0.}

You'll notice that there are a few elements that match your answers, but most of them don't. I found the pattern of ones that matched quite telling (basically corresponding to the diagonal elements of the eigenvectors matrix) and, indeed, if I transpose my eigenvectors matrix, I get the answer you got (up to those irrelevant phase factors).