Outer function errors, Mutual Information

[PZuiuc]

In the following code, after obtaining the ground state psi of the Hamiltonian H using DMRG, I verify the normalization by checking inner(psi, psi) and trace(rho), where rho = outer(dag(psi), psi). I find that inner(psi, psi) equals 1 up to an error on the order of 10⁻¹⁴, while trace(rho) deviates from 1 by about 10⁻⁶. I would like to understand the source of this discrepancy. Specifically, does the outer function automatically introduce truncation or compression that could explain the larger error in trace(rho)? Or might there be a mistake in how I’m computing the trace of rho? Thanks!

Nuc=60
sites = siteinds("S=1/2",Nuc)

J=1.0
os = OpSum()
for j=1:Nuc-1
      if j==Nuc
          jp=1
      else
          jp=j+1
      end
    os += -4J,"Sz",j,"Sz",jp
    os += -2,"Sx",j
 end
 os += -2,"Sx",Nuc

 H = MPO(os,sites)

 psi0 = productMPS(sites, n -> isodd(n) ? "Up" : "Dn")

 nsweeps = 15
 maxdim = [10,20,100,100,200,200,300,300,400,500]
 noise=[1E-3,1E-5,1E-7,1E-7,1E-9,1E-9,1E-11,1E-11,1E-13,1E-13,1E-14,1E-15]
 cutoff = [1E-12]

 energy0, psi = dmrg(H,psi0; nsweeps, maxdim, noise, cutoff)
 orthogonalize!(psi, 1)
 rho0=outer(psi',psi)

 println(inner(psi,psi))

 aaa=ITensor(1.0)
 for i=1:Nuc
     rho0[i]=rho0[i]*delta(sites[i],sites[i]')
     aaa*=rho0[i]
 end
 println(aaa)

[miles]

I’m not 100% sure, but your guess about the role of outer here is likely what’s going on. When I ran your code I found that outer gives an MPO with a bond dimension that is less than the square of the MPS bond dimension than is inputted into it. Perhaps this is ok if there are a lot of small or zero singular values of rho0 but it could be that these are being overly truncated.

In general I would not recommend the use of outer except for basic testing on small systems. A much more accurate and efficient way to work with pure density matrices (outer product of an MPS with itself) is to think of forming such a density matrix implicitly, then computing properties like the trace using an algorithm similar to the algorithm for computing an inner product of two MPS.

If you first right-orthogonalize the MPS, as you have done in your code, then actually all the tensors except the first one cancel out if you take a partial trace of \psi \otimes \psi^\dagger on all the sites except site 1. So then the trace of the entire density matrix reduces to just the inner product of the first MPS tensor with itself.

[PZuiuc]

Hi Miles,
Thank you very much for your swift response. Since my goal is to compute the mutual information for a model more complex than the Ising model, I think constructing the full density matrix is indeed necessary. In this case, do you know of any reliable method within ITensor to build the density matrix accurately?

[miles]

There isn’t a function that you can simply call to obtain the density matrix like this. One reason is that it just doesn’t have an efficient form as an MPS or an MPO. It’s only efficient form is as two MPS in a (“lazy” or formal or implicit) outer product with each other.

What sort of regions do you want to compute the mutual information between? There are efficient ways to compute that if the two regions are sufficiently small or have other structure like if they reach to the edges of the system.

There are also some previous discussions on here asking about mutual information though the answers may not be totally satisfactory.

Finally there is an experimental package ITensorEntropyTools which may help you.
< GitHub - ryanlevy/ITensorEntropyTools.jl: Tools for calculating Entanglement Entropy w/ ITensors.jl>

[PZuiuc]

Thanks a lot for sharing this resource. The mutual information I want to calculate is between two spins inside a chain, which are far away from each other. I will try to look at the ITensorEntropyTools to see if I can find solution.

[miles]

Yes, the mutual information between two individual sites is something that is efficiently computable with MPS. But the best approach is a more “custom” or “targeted” method for computing the reduced density matrices involved, rather than ever forming the entire (pure) density matrix of the whole system as a single tensor network like an MPO.

Here is what I mean. The diagrams shown below depict how the three reduced density matrices involved in the mutual information can be written as efficient tensor network contractions (efficient because you can contract all of these from the ‘side’ i.e. from left to right, alternating between the top and bottom layers).

As a possible optimization, if the MPS is “gauged” so that its orthogonality center is at i for the first diagram, j for the second diagram, and somewhere in the range [i,j] for the third diagram, then the outer MPS tensors ‘cancel’ and can be omitted from the computations. Then for \rho_i for example, only the i th MPS tensor is needed for the computation.

So some pseudocode (i.e. to be tested and checked) for the first RDM is:

psi = orthogonalize(psi,i)
rhoi = prime(psi[i],"Site") * dag(psi[i])

Finally, once each RDM is computed, they can all be diagonalized using eigen to find their spectra.