I am trying to compute the trace distance of two general $4 \times 4$ density matrices as such:

$D=\frac{1}{2} \, \mathrm{tr} \, |\Delta\rho|_1$ where $\Delta\rho$ is the difference between two density matrices $\rho_1, \rho_2$ and $|A|_1=(A^\dagger A)^{1/2}$. Since density matrices are Hermitian one may write $|\Delta\rho|_1=(\Delta\rho^2)^{1/2}$ hence one ends up with $D=\frac{1}{2} \, \mathrm{tr} \,|\Delta\rho|_1=\frac{1}{2}\sum_i |\lambda_i|$ where the $|\lambda_i|$ are the eigenvalues of the density matrix $\Delta\rho$.

Is there an efficient manner in which one may calculate this, especially for higher dimensional density matrices?

  • $\begingroup$ Can you show an example of rho1, rho2 and give us an idea of what "higher dimensions" might be? $\endgroup$
    – MarcoB
    Commented Dec 29, 2020 at 16:44
  • $\begingroup$ @MarcoB: The rho's are any Hermitian matrices which have trace 1. By "higher dimensions"I mean larger density matrices for example of a multilevel or a composite system (they tend to get very messy to calculate). $\endgroup$
    – dan
    Commented Dec 29, 2020 at 18:31

1 Answer 1

randomDensityMatrix[n_] := ConjugateTranspose[#2].(#1 #2) &[
  Normalize[RandomReal[{1, 10}, n], Total], 

n = 100;
ρ1 = randomDensityMatrix[n];
ρ2 = randomDensityMatrix[n];

Both the following may work:

0.5 Re[Tr[MatrixPower[ConjugateTranspose[#].# &[ρ1 - ρ2], 1/2]]]
0.5 Total[Sqrt[Eigenvalues[ConjugateTranspose[#].# &[ρ1 - ρ2]]]]

But the second version seems to be ten times faster for n = 100.


J.M. pointed out in a comment that Total[Abs[Eigenvalues[ρ1 - ρ2]]]/2 might be a better idea because the condition number of ConjugateTranspose[ρ1 - ρ2].(ρ1 - ρ2) might be quite high. I objected that it actually runs slower and guessed that would happen because ρ1 - ρ2 is indefinite. I have to correct this guess: Now I believe that the speed difference might be caused by Mathematica not correctly detectling that ρ1 - ρ2 is Hermitian. The runtimes below somewhat indicate that Apparently, operations like ConjugateTranspose[#].#& and ConjugateTranspose[#] + #&set some internal flag that makes Eigenvalue branch to a special algorithm for Hermitian matrices:

n = 100;
ρ1 = randomDensityMatrix[n];
ρ2 = randomDensityMatrix[n];

0.5 Re[Tr[MatrixPower[ConjugateTranspose[#].# &[ρ1 - ρ2], 1/2]]] // RepeatedTiming
0.5 Total[Sqrt[Eigenvalues[ConjugateTranspose[#].# &[ρ1 - ρ2]]]] // RepeatedTiming
Total[Abs[Eigenvalues[ρ1 - ρ2]]]/2 // RepeatedTiming
Total[Abs[SingularValueList[ρ1 - ρ2]]]/2 // RepeatedTiming

B = ρ1 - ρ2;
B = 1/2 (ConjugateTranspose[B] + B);

Total[Abs[Eigenvalues[B]]]/2 // RepeatedTiming
Total[Abs[SingularValueList[B]]]/2 // RepeatedTiming

{0.013, 0.289995}

{0.00075, 0.289995}

{0.0054, 0.289995}

{0.0014, 0.289995}

{0.00064, 0.289995}

{0.00129, 0.289995}

As you can see, Eigenvalues and SingularValueList operate much faster on the matrix B as on ρ1 - ρ2. So

Total[Abs[Eigenvalues[ConjugateTranspose[#] + # &[ρ1 - ρ2]]]]/4

seems to be the better choice, no matter what the condition number of ρ1 - ρ2 is.

  • 1
    $\begingroup$ Since ρ1 - ρ2 is Hermitian, Total[Abs[Eigenvalues[ρ1 - ρ2]]]/2 (equivalently, Norm[Eigenvalues[ρ1 - ρ2], 1]/2) should also work, which also avoids the formation of ConjugateTranspose[matrix].matrix. (Otherwise, if ρ1 - ρ2 was not Hermitian, SVD would be necessary.) $\endgroup$ Commented Dec 29, 2020 at 17:17
  • $\begingroup$ @J.M. Yesss, good point... ConjugateTranspose[matrix].matrix may have a much worse condition number than matrix itself. But Total[Abs[Eigenvalues[ρ1 - ρ2]]]/2 seems to branch to a less efficient algorithm for computing the eigenvalues because ρ1 - ρ2 is indefinite... $\endgroup$ Commented Dec 29, 2020 at 17:20
  • 1
    $\begingroup$ In that case, the choice boils down to whether ρ1 - ρ2 is well-conditioned (i.e. the ratio of the largest to smallest eigenvalue is not overly large) or not. If it is not badly conditioned, you could get away with forming ConjugateTranspose[matrix].matrix. $\endgroup$ Commented Dec 29, 2020 at 17:23
  • $\begingroup$ Another things that's odd: Total[Abs[SingularValueList[ρ1 - ρ2]]]/2 is actually faster than Total[Abs[Eigenvalues[ρ1 - ρ2]]]/2... oO $\endgroup$ Commented Dec 29, 2020 at 17:25
  • $\begingroup$ @J.M. I think I figured out the reason for the performance degradation and edited my post. $\endgroup$ Commented Dec 29, 2020 at 18:10

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.