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.