Calculating fidelity with iMPS

Question

Edit: I have made some progress with this but I'm on holiday at the moment. I will add something in the next couple of weeks.
Edit 2: I came up with the following code to implement the transfer matrix, which can then be used with arnoldi. However it runs quite slowly for larger unit cells and bond lengths. Please note that I've only tested this code briefly and it should be used with caution

#ifndef __TRANSFER_H_
#define __TRANSFER_H_

#include "itensor/all.h"

using namespace itensor;
using namespace std;

struct TMatrix {
  ITensor left, right;
  ITensor T;
  int N;
  Index link1, link2;
  Index left1, left2, right1, right2;
  TMatrix(MPS psi1, MPS psi2) {
    N = length(psi1);

    // Use primes to distinguish some indices
    psi1.position(1);
    psi2.position(1);
    link1 = commonIndex(psi1.A(1), psi1.A(N), {"Link"});
    link2 = commonIndex(psi2.A(1), psi2.A(N), {"Link"});
    left1 = link1;
    left1.prime(1);
    left2 = link2;
    left2.prime(2);
    right1 = link1;
    right1.prime(3);
    right2 = link2;
    right2.prime(4);
    T = dag(psi1.A(1)) * delta(link1, left1);
    T *= psi2.A(1) * delta(link2, left2);
    auto s1 = findIndex(psi1.A(1), "Site");
    auto s2 = findIndex(psi2.A(1), "Site");
    if (s1 != s2)
      T *= delta(findIndex(psi1.A(1), "Site"), findIndex(psi2.A(1), "Site"));
    for (int i = 2; i <= N; i++) {
      ITensor z = dag(psi1.A(i));
      if (s1 != s2)
        z *= delta(findIndex(psi1.A(i), "Site"), findIndex(psi2.A(i), "Site"));
      if (i == N)
        z *= delta(link1, right1);
      T *= z;

      z = psi2.A(i);
      if (i == N)
        z *= delta(link2, right2);
      T *= z;
    }
  }

  ITensor make_guess() { return randomITensor(right1, right2); }
  size_t size() { return dim(right1) * dim(right2); }

  void product(ITensor &other, ITensor &result) {
    result = other * T;
    result *= delta(right1, left1);
    result *= delta(right2, left2);
  }
};

#endif // __TRANSFER_H_

Dear Miles/Matt/ITensor,

I've been using your iDMRG example code (idmrg.h etc) and I've adapted it to some models I'm interested in. My goal is to use the fidelity and fidelity per site (and their susceptibilities) to examine the phase diagram. For a finite system, the fidelity is simply the overlap of the ground states at two different values of some parameter of the Hamiltonian. I'm able to calculate the "overlap" of two iMPS's with some code adapted from your examples - this required some modifications to make the site indices match if they were created from different SiteSet instances. I believe this is appropriate for the fidelity but I'm not sure, and it offers no way to calculate the fidelity per site. I've added this code below.

In McCulloch's, "Infinite size density matrix renormalization group, revisited" (section D. Fidelity), he suggests a method using the eigenvalues of a transfer matrix Eq. (39). I'm having difficulty understanding this section, in particular what 'E' denotes in Eq. (39). I've gone through the references but I haven't been able to decipher it. However, McCulloch seems to imply the procedure is fairly simple.

So, I was wondering if you have any suggestions on how to calculate the fidelity (per site) from iMPS's.

Cheers,
Alex Henry

Complex ioverlap(MPS &psi1, MPS &psi2, bool same) {
  if (same) {
    // if psi1 and psi2 already have the same site indices
    auto O = dag(psi1(1)) * psi2(1);
    for (auto n : range1(2, length(psi2)))
      O = O * dag(psi1(n)) * psi2(n);
    auto ovrlap = real(eltC(O));
    return ovrlap;
  }
  auto O = dag(psi1(1)) * psi2(1);
  auto l1 = findIndex(dag(psi1(1)), "Site");
  auto l2 = findIndex(psi2(1), "Site");
  O = O * delta(l1, l2);
  for (auto n : range1(2, length(psi2))) {
    auto a = dag(psi1(n));
    auto b = psi2(n);
    l1 = findIndex(a, "Site");
    l2 = findIndex(b, "Site");
    auto link = delta(l1, l2);
    a = a * link;
    O = O * a * b;
  }
  auto overlap = real(eltC(O));
  return overlap;
}

Answer 1

The short answer is that in that equation, E is the dominant eigenvector of the transfer matrix. When you write down the overlap of an infinite MPS in canonical form with itself, it reduces to repeated unit cells of the transfer matrix. A helpful place to look to understand this is Section 5 of: https://arxiv.org/pdf/1306.2164.pdf. So you would have to transform the MPS that is output from iDMRG into canonical form, and then use the canonical form MPS to compute the dominant eigenvector of the transfer matrix. Infinite MPS techniques are generally more technical and subtle than finite MPS techniques. Have you tried using finite MPS techniques first, and what limitations did you see?

We are working on an infinite MPS package in the Julia version which would make performing calculations like this easier (https://github.com/mtfishman/ITensorInfiniteMPS.jl) but it is not quite ready for general usage yet, I'm hoping to have some basic functionality like what you are asking for ready in the next month or so.

Comments

Thanks for your response, Matt.

I was able to figure some of that out, and I've added some code I wrote to the original post above. It appears to give appropriate results, but it runs fairly slowly for larger unit cells (N~6 seems to be necessary for the systems I'm working with). There is probably a more efficient way to do the contraction, but since the MPS is periodic it's not the same as the basic examples. I may be able to figure this out myself. I will have a look at that reference and at your Julia code.
Edit: to clarify, the construction of the tensor T in the above code is what runs slowly, taking 10 or more times longer than the arnoldi routine.

I have used finite DMRG already quite extensively. The reason I'm using iDMRG is that the systems I'm working with (Zn-symmetric Potts-type models) are computationally intensive and don't converge to the thermodynamic limit quickly, seemingly requiring some hundreds of sites. iDMRG already allowed me to compute the (L->infinity) ground state energy much more efficiently, but getting the fidelity turned out to be harder than I expected. The main goal is to use the fidelity per site to examine the phase diagrams of these systems.