Eigenvalue decomposition of a density matrix not reproducing original density matrix

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.

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