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
    Dec 29 '20 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
    Dec 29 '20 at 18:31
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$
    – J. M.'s torpor
    Dec 29 '20 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$ Dec 29 '20 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$
    – J. M.'s torpor
    Dec 29 '20 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$ Dec 29 '20 at 17:25
  • $\begingroup$ @J.M. I think I figured out the reason for the performance degradation and edited my post. $\endgroup$ Dec 29 '20 at 18:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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