How to compute the trace distance of a density matrix

Question

I am trying to compute the trace distance of two general $4 times 4$ density matrices as such:
$D=frac{1}{2} , mathrm{tr} , |Deltarho|_1$ where $Deltarho$ 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
$|Deltarho|_1=(Deltarho^2)^{1/2}$ hence one ends up with $D=frac{1}{2} , mathrm{tr} ,|Deltarho|_1=frac{1}{2}sum_i |lambda_i|$ where the $|lambda_i|$ are the eigenvalues of the density matrix $Deltarho$.
Is there an efficient manner in which one may calculate this, especially for higher dimensional density matrices?

Henrik Schumacher · Accepted Answer

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.

How to compute the trace distance of a density matrix

One Answer

Add your own answers!

Ask a Question