5
$\begingroup$

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?

Improve this question
$\endgroup$
2
  • $\begingroup$ Can you show an example of rho1, rho2 and give us an idea of what "higher dimensions" might be? $\endgroup$ 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$ Commented Dec 29, 2020 at 18:31

1 Answer 1

Reset to default
6
$\begingroup$ 
randomDensityMatrix[n_] := ConjugateTranspose[#2].(#1 #2) &[ Normalize[RandomReal[{1, 10}, n], Total], RandomVariate[CircularUnitaryMatrixDistribution[n]] ] 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.

Edit:

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; RandomSeed[20201229]; ρ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.

Improve this answer
$\endgroup$
6
  • 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
  • 1
    $\begingroup$ Nice detective work! ;) Seeing an SVD proceed faster than an eigendecomposition is surely pretty sus (as the kids might say these days) in itself. $\endgroup$ Commented Dec 30, 2020 at 9:42

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.