# Need some help with Purity calculation

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

$$$$\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}$$$$

Then I calculated the Eigenvectors and Eigenvalues of this matrix:

Eigenvalues: $$$$e=[0.79987375 \space\space 0.16872495 \space\space 0.03140131\space\space 0 ]$$$$ Eigenvectors: $$$$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}$$$$

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: $$$$computational bases= \begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{pmatrix}$$$$

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

$$$$\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}$$$$

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: $$$$\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}$$$$

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.

• Are the eigenvectors the rows or the columns of the matrix you gave? Commented Apr 4, 2020 at 18:16
• The eigenvectors are the columns of the matrix. Commented Apr 4, 2020 at 18:30
• In which case your $\Psi$ doesn’t look right. The first and the fifth elements should correspond to the first two elements of the first eigenvector, but the relative sizes are wrong. Commented Apr 5, 2020 at 5:48
• ok, thank you for your answer, but why have the first and the fith element of $\Psi$ has to be equal to the first two elements of the first eigenvector? I edited my post, so you can see, how I calculated $\Psi$. Maybe you could tell me where I did a mistake? Commented Apr 5, 2020 at 17:59

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}]


{-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.}