# How to compute the trace distance of a density matrix

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?

• Can you show an example of rho1, rho2 and give us an idea of what "higher dimensions" might be? – MarcoB Dec 29 '20 at 16:44
• @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). – dan Dec 29 '20 at 18:31

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.

• 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.) – J. M.'s ennui Dec 29 '20 at 17:17
• @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... – Henrik Schumacher Dec 29 '20 at 17:20
• 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. – J. M.'s ennui Dec 29 '20 at 17:23
• Another things that's odd: Total[Abs[SingularValueList[ρ1 - ρ2]]]/2 is actually faster than Total[Abs[Eigenvalues[ρ1 - ρ2]]]/2... oO – Henrik Schumacher Dec 29 '20 at 17:25
• @J.M. I think I figured out the reason for the performance degradation and edited my post. – Henrik Schumacher Dec 29 '20 at 18:10