efficient calculation of MPO-MPO multiplication diagonal elements/trace

Hi!

I need to calculate the trace of the multiplication of two MPOs, namely tr(A::MPO, B::MPO). The naive procedure, e.g.,

using ITensors

sites = siteinds("SpinHalf", 3)
A = randomMPO(sites)
B = randomMPO(sites)
println("tr(A*B) = ", tr((A*prime(B)); plev=Pair(0,2)))

works well, but is highly inefficient. It first calculates all of the elements of A*B, while only the diagonals are needed.
Is there a way to contract only the relevant matrices and get the diagonal elements only?

I am using Julia v1.7.0 and ITensors v0.3.18

Many thanks in advance

-Lidia

Good question. The easiest thing to do right now is:

inner(swapprime(dag(A), 0 => 1), B)

That is slightly more inefficient than it needs to be for complex tensors, since inner calls dag internally on the first input so we have to undo it by passing dag(A), but the scaling will be better than tr(apply(A, B)) (equivalent to the code you wrote above) either way.

I’m considering adding an interface like: tr(@lazy apply(A, B)) which would perform the trace efficiently without performing the intermediate MPO-MPO contraction apply(A, B).

2 Likes

thank you very much!