Trace distance of two MPS


I want to get Trace distance D which is defined as

where p and q are two quantum states(MPS).
Now I use


to get two MPO and make


How can I finish the second line in picture with ITensor? The X is a MPO and dagger(X)*X is also a MPO. I have to do a sqrt and trace to a it. I have tried some ways but failed.

Unfortunately I don’t know if there’s an efficient way to get \text{Tr}\sqrt{X^\dagger X} for an MPO. There is always the inefficient way of multiplying all the MPO tensors together and diagonalizing the resulting matrix to get the eigenvalues but that will scale exponentially with number of sites.

Here is a thought though: X^\dagger X will be an operator whose spectrum is the magnitude squared of all the eigenvalues of X. Then taking the square root will take the absolute value. So if X has a positive, real spectrum then you can leave out the squaring and square root steps and just take the trace of X. But I don’t think that will hold for X=p-q in general because q enters with a minus sign.

Lastly, would you consider using the 2-norm instead of 1-norm for your problem? The 2-norm is much more natural to compute and could be computed efficiently.

I try to diagonalize the resulting matrix to get the eigenvalues and sometimes it works. Fortunately, my number of sites is small so the cost is acceptable. However, when I transform the MPO to an ITensor, the number of “Index” of resulted ITensor has to be less than 14. So the sites of MPO must be less than 7.

When you say it has to be less than 14, do you mean the code runs too slowly otherwise or that you get an error message from ITensor?

I think I got an error message from ITensor. When I set 7 sites MPO(it will be made to ITensor with 14 Index), I got:

Contraction resulted in ITensor with 14 indices, which is greater than or equal to 
the ITensor order warning threshold 14. You can modify the threshold with macros 
like `@set_warn_order N`, `@reset_warn_order`, and `@disable_warn_order` 
or functions like `ITensors.set_warn_order(N::Int)`, `ITensors.reset_warn_order()`, and `ITensors.disable_warn_order()`.

  [1] _contract(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/KifqM/src/itensor.jl:1778
  [2] contract(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/KifqM/src/itensor.jl:1872
  [3] *(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/KifqM/src/itensor.jl:1860
  [4] top-level scope
    @ ~/code/try_trace.jl:33
  [5] eval
    @ ./boot.jl:373 [inlined]
  [6] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1196
  [7] include_string(m::Module, txt::String, fname::String)
    @ Base ./loading.jl:1206
  [8] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base ./essentials.jl:716
  [9] invokelatest(::Any, ::Any, ::Vararg{Any})
    @ Base ./essentials.jl:714
 [10] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/eval.jl:211
 [11] (::VSCodeServer.var"#65#69"{Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/eval.jl:155
 [12] withpath(f::VSCodeServer.var"#65#69"{Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/repl.jl:184
 [13] (::VSCodeServer.var"#64#68"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/eval.jl:153
 [14] hideprompt(f::VSCodeServer.var"#64#68"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/repl.jl:36
 [15] (::VSCodeServer.var"#63#67"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/eval.jl:124
 [16] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging ./logging.jl:511
 [17] with_logger
    @ ./logging.jl:623 [inlined]
 [18] (::VSCodeServer.var"#62#66"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/eval.jl:201
 [19] #invokelatest#2
    @ ./essentials.jl:716 [inlined]
 [20] invokelatest(::Any)
    @ Base ./essentials.jl:714
 [21] macro expansion
    @ ~/.vscode/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/eval.jl:34 [inlined]
 [22] (::VSCodeServer.var"#60#61")()
    @ VSCodeServer ./task.jl:423

But it works for 6 sites MPO.

Thanks. So if you take a look at the error message, what it is saying is that we have a limit of 14 indices that is turned on by default. But you can remove this limit by following the steps explained in the message.

Yeah, 6 sites is sufficient for me. If I need more I will try to remove this limit. Thanks for your help.

I realize this is an old question, but since you are considering state vectors (pure states), this is completely unnecessary. For a pair of pure states the trace distance is just a function of the overlap.
See, e.g., the last property on Wikipedia: Trace distance - Wikipedia

Anyway, I found this question because I’m looking how to compute the trace norm of an MPO efficiently.
However, it seems that @miles has replied to my question and one goes back to exponential scaling.

But since this is a very active field, I guess it doesn’t hurt to ask: was there any development in the past year about this computation?
I’ll mention that the trace norm can be written in a variational form (particularly also as a semidefinite program), but I am not aware of a method to write such an optimization as a tensor network algorithm…

Actually, the above response about trace distance has helped me solve the problem.

If you talk about the norm of density matrix, I try to transform a “ITensor” to a tensor and get a matrix formulation. So I can get the trace. Above process works for systems. (I have tried about 8 sites)

I’ve never seen the trace norm. Is it related to the trace of density matrix?

For sure you can compute the trace distance between two pure states as discussed in the previous replies, but I think it’s unnecessarily complicated, since it is just a simple function of the overlap between the two states.

The trace norm is exactly the norm ||X||_1 in your original question. It’s just a name popular in quantum information for the Schatten-1 matrix norm (also known as nuclear norm, i.e. sum of the singular values of the operator).

Basically I’d like to know if there is a way to compute the trace norm of an arbitrary MPO, beyond a few sites, i.e. with a complexity that is not exponential in the number of sites.

I may not know a simple way to compute the trace norm of any MPO. I just get a big tensor after contraction and use the function “reshape” to get a matrix. It is not a smart way but it is sufficient
to may small system.

If you are interested in above program, here are mine and hope it give you some help.


for i1=1:2
    for j1=1:2
        for k1=1:2
            for m1=1:2
                el = mm[d1=>i1,d2=>j1,d3=>k1,d4=>m1]
                be[i1,k1,j1,m1]= el


Here, “M1” is a very small MPO with four legs. Then the “af” is density matrix we want and you can get the trace.

One can definitely compute the trace of an MPO efficiently. If my understanding is correct, for symmetric, positive semi definite MPO’s (viewed as a large matrix), then the trace will be the same as the trace norm ||X||_1.

For a more general MPO that isn’t symmetric and positive semi-definite, I don’t think there is a general sub-exponential way to compute this norm. But if the spectrum of singular values falls quickly enough, I could envision using DMRG to compute the squares of the dominant singular values by squaring the MPO and inputting it into the DMRG algorithm. Then one could take the square roots of these values and add them up.

Yes, the trace norm is simply the trace for a positive semidefinite matrix.
In most cases of interest in quantum information one wants to compute the trace norm of MPOs that are not positive definite (the most common case being the trace distance in the original question), and in some cases not even Hermitian.

I am not sure if I would expect the singular values to decay quickly, but the suggestion to treat X^\dag X as an Hamiltonian and estimate its spectrum is certainly intriguing.
Thanks a lot for the input!

Glad that’s helpful. In general, tensor networks are only efficient for objects (e.g. operators) where there is some notion of singular values that decay quickly along some kind of “cut” (i.e. decomposition or bipartition) so if you are considering operators that don’t have that property, a tensor network may just not be the right tool.

(Take this with a grain of salt, however, because some kinds of ‘loopy’ tensor networks like PEPS can efficiently represent objects with a huge number of large singular values along a 1D cut, if you think of the cut as dividing a 2D system into two halves, say. Also tensor networks can help in cases where there is low-rank structure between spatial scales, so it’s not always necessary to have low-rank structure across different spatial regions as the only situation in which tensor networks help.)

1 Like

I wanted to provide some quantum information background that might help clarify the difficulty in computing trace distance (TD) versus Hilbert-Schmidt distance (HSD) in general. For clarity, TD is given by \frac{1}{2}||\rho - \sigma||_1 whereas HSD is given by ||\rho - \sigma||_2 for ||\cdot||_p is Schatten p-norm..

In short, computing HSD is always easier than computing TD. For pure states (and hence MPS), they are equivalent distances measures, so it is better to compute HSD. For mixed states (and hence generally MPO density matrices), they can be quite different. Nevertheless, it is known that

||\rho - \sigma||_2 \leq ||\rho - \sigma||_1 \leq \sqrt{\text{rank}(\rho - \sigma)} || \rho - \sigma||_2

from monotonicity of norms and Cauchy inequality. So in worst case all you can say ||\rho - \sigma||_1 \leq \sqrt{d} ||\rho - \sigma||_2 for d \times d density matrices. In case where computing TD is too hard (common), best one can do (without additional info on \rho, \sigma) is compute ||\rho - \sigma||_2 as a loose upper bound on TD.