Eigenvalue decomposition of a density matrix not reproducing original density matrix

Question

A density matrix $\rho$ in quantum mechanics is defined as any self adjoint and positive semidefinite matrix with a trace or 1. It can be expanded into sets of pure states such that $\rho=\sum_{i}p_{i} \vert \phi_{i}\rangle\langle \phi_{i}\vert$ where the $p_{i}$ sum to unity and the $\phi_{i}$ are not unique and not even necessarily orthogonal. However, in the special case that the $\phi_{i}$ are orthogonal they also correspond to the eigenvectors of $\rho$ with the $p_{i}$ the eigenvalues.

When I try to perform this decomposition on random density matrices in Mathematica I recover a valid density matrix however it is not equivalent to the original density matrix. I am curious if there is some subtlety related to the way Mathematica orders or normalizes eigenvectors which is causing this?

Below I give a specific example. I know the code isn't elegant but it should be transparent at least. The initial matrix we wish to decompose is defined as:

testm = {{0.2146034795560191` + 0.` I, -0.06085127765739766` - 
 0.10901433064555709` I, 
0.003998568237967147` + 
 0.1698381572412944` I, -0.19626623770083873` + 
 0.06355374075232331` I}, {-0.06085127765739767` + 
 0.10901433064555709` I, 0.28539652044398106` + 0.` I, 
0.04070225716833632` - 
 0.17141179642617063` I, -0.003998568237967189` - 
 0.16983815724129442` I}, {0.003998568237967161` - 
 0.16983815724129436` I, 
0.040702257168336316` + 0.1714117964261706` I, 
0.2853965204439809` + 0.` I, 
0.06085127765739764` + 
 0.10901433064555702` I}, {-0.1962662377008387` - 
 0.06355374075232331` I, -0.003998568237967209` + 
 0.16983815724129442` I, 
0.06085127765739763` - 0.10901433064555702` I, 
0.21460347955601905` + 0.` I}};

I decompose it using Eigensystem to find the four $\vert \phi_{i}\rangle\langle\phi_{i}\vert$ labelled here as $t_{i}$:

t1 = ConjugateTranspose[{Eigensystem[testm][[2]][[1]]}].{Eigensystem[
   testm][[2]][[1]]};

t2 = ConjugateTranspose[{Eigensystem[testm][[2]][[2]]}].{Eigensystem[
   testm][[2]][[2]]};

t3 = ConjugateTranspose[{Eigensystem[testm][[2]][[3]]}].{Eigensystem[
   testm][[2]][[3]]};

t4 = ConjugateTranspose[{Eigensystem[testm][[2]][[4]]}].{Eigensystem[
   testm][[2]][[4]]};

Which I then sum up using the corresponding eigenvalues as the weights:

testmdecomposition = Eigensystem[testm][[1]][[1]]*t1 + Eigensystem[testm][[1]][[2]]*t2 + Eigensystem[testm][[1]][[3]]*t3 + Eigensystem[testm][[1]][[4]]*t4

Although the final matrix is a valid density matrix it varies substantially in the off diagonals from the original matrix.

Answer 1

One can read in the documentation of Eigensystem that

{Λ, U} = Eigensystem[testm];

satisfy

testm == Transpose[U].DiagonalMatrix[Λ].Conjugate[U]

and not

testm == ConjugateTranspose[U].DiagonalMatrix[Λ].U

With finite precision numbers, we can check that with

Max[Abs[testm - Transpose[U].DiagonalMatrix[Λ].Conjugate[U]]]

8.92542*10^-16

Alternatively, you obtain the projectors onto the eigenspace with

t = Table[ KroneckerProduct[U[[i]], Conjugate[U[[i]]]], {i, 1, Length[testm]}];

and reobtain the matrix testm with

Λ.t