Implementing efficient tensor contraction algorithm

Hi Miles, Matt and other admins,

I recent wrote a paper about calculating the entanglement asymmetry in a generic MPS scenario (without integrability, interacting spin-1/2 Ising Hamiltonian, in this paper ( https://arxiv.org/pdf/2312.08601.pdf )).

I calculated the entanglement asymmetry in Eqn 7 (from what I call symmetry projected reduced density matrix, Eqn 6) using the ‘brute-force’ way of constructing the reduced density matrix MPO first, and then perform MPO-MPO contraction and tracing in O(D^8) in the worst case (in Appendix D of my paper). In my opinion, in particular, forming an MPO reduced density matrix will make the bond dimension D^2 and I thought that scales up complexity quickly.

I outlined an alternative way of performing tensor contraction by dealing with individual tensor components (see Fig 15). I was wondering if there is a way to implement Fig 15 in Julia, and if so, what the pseudo-code would look like?

Best,
Brian.

Hi Brian,
From the diagram in the paper (copied below for others to see), it looks like the main network is something like an MPS but which has a dangling virtual bond at the right-hand edge. Am I correct in thinking that the four horizontal MPS-like networks are all the same as each other?

As far as pseudo-code, you can just do a pattern like this (say if you want to contract two different MPS psi and phi):

L = ITensor(1.)
for j=1:N-1
L *= psi[j]
L *= dag(phi[j])
end

which would gradually contract and build up the large “block” tensor you show in the lower part of the diagram below. Lastly you can contract with the matrices shown in the last steps of your diagram. The code above might have to be adjusted if psi and phi are just copies of the same MPS to avoid accidentally contracting their shared virtual indices.

Does that help you get started? It’s hard to give a more specific answer without knowing more details about e.g. the assumptions on the individual tensors (are they all different in each row in your diagram, or can some be assumed to be the same, etc).

Hi Miles,

Thanks for the clarification. One more question, for identical MPS to have different index ID for the same leg, do I use deepcopy(psi) to create a copy that has completely different index ID for the same leg?

Best,
Brian.

deepcopy is a perfect copy of the MPS, so it has the same indexes. The method you are searching for is sim[!], I can’t find it on ITensors Documentation, but if you type from Julia REPL ?sim in the output you would find:

  sim[!](linkinds, M::MPS, args...; kwargs...)
  sim[!](linkinds, M::MPO, args...; kwargs...)

  Apply sim to all link indices of an MPS/MPO, returning a new MPS/MPO.

  The ITensors of the MPS/MPO will be a view of the storage of the original ITensors.

  ───────────────────────────────────────────────────────────────────────────────────────────────────────

  sim[!](siteinds, M::MPS, args...; kwargs...)
  sim[!](siteinds, M::MPO, args...; kwargs...)

  Apply sim to all site indices of an MPS/MPO, returning a new MPS/MPO.

  The ITensors of the MPS/MPO will be a view of the storage of the original ITensors.

  ───────────────────────────────────────────────────────────────────────────────────────────────────────

  sim[!](siteinds, commoninds, M1::MPO, M2::MPS, args...; kwargs...)
  sim[!](siteinds, commoninds, M1::MPO, M2::MPO, args...; kwargs...)

  Apply sim to the site indices that are shared by M1 and M2.

  Returns new MPSs/MPOs. The ITensors of the MPSs/MPOs will be a view of the storage of the original
  ITensors.

  ───────────────────────────────────────────────────────────────────────────────────────────────────────

  sim[!](siteinds, uniqueinds, M1::MPO, M2::MPS, args...; kwargs...)

  Apply sim to the site indices of M1 that are not shared with M2. Returns new MPSs/MPOs.

  The ITensors of the MPSs/MPOs will be a view of the storage of the original ITensors.

This creates a view of the original tensor. If what you need is a deepcopy of an ITensor A with different indices you can just use a combination of deepcopy and sim. So if you only need the 2 MPSs to have different site indices id you can do

A_copy = sim(siteinds, deepcopy(A))