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[

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

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

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

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.


1 Answer 1


One can read in the documentation of Eigensystem that

{Λ, U} = Eigensystem[testm];


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


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.