Implementing local basis optimization

Hi there! I am trying to implement local basis optimization after this paper. I’ve done my best to follow their algorithm (figure 2 is very helpful in this respect), but I am definitely doing something wrong – my GS energy for an 8-site Hubbard Holstein system is on the order of +25, while I know that it should be closer to -2. Here is the gist of what I am doing:

function dmrg(...)
# Randomly initialize a rotation matrix (R) on every site
# All other pre-DMRG steps, like calling orthogonalize! and position! on site 1
for sw in 1:nsweep(sweeps)
 for (b, ha) in sweepnext(N) 
  if do_LBO
   # Perform LBO (see below)
  end
  # Perform the remainder of 2-site DMRG using ITensors code
 end
end

Now, I will get into more detail and include snippets of actual code. Here is how I am randomly initializing my rotation matrices:

 if LBO
    Rs = []
    for i in 1:N
      # Start off with having the same dimension 
      S_ind = getfirst(x -> hastags(x, "Site"), inds(psi[i]))
      Trunc_ind = settags(copy(S_ind),"Trunc")
      R = randomITensor(prime(S_ind), dag(Trunc_ind))
      push!(Rs,R)
    end
  end

And here is my LBO implementation:

## HELPER: Get the left indices of the tensor we want to SVD##
function get_Linds(O::ITensor, tag::String)
    Rind = getfirst(x -> hastags(x, tag), inds(O))
    return noncommoninds([Rind], inds(O))
end

# Put into canonical form
orthogonalize!(psi,b)
M = psi[b] 
R = Rs[b]
H = PH.H[b]

# Perform the first SVD
X, Λ, Y = svd(M, get_Linds(M, "Site"), lefttags="svd1,u", righttags="svd1,v")
u = getfirst(x -> hastags(x, "svd1,u"), inds(Λ))
s = getfirst(x -> hastags(x, "Site"), inds(Y))
t = getfirst(x -> hastags(x, "Trunc"), inds(R))

# Make and optimize R̃
R̃ = R * Y * Λ * dag(delta(t,s))
X̃, Λ̃, R = svd(R̃, u, lefttags="svd2,u", righttags="Trunc")

# Perform updates
t = getfirst(x -> hastags(x, "Trunc"), inds(R))
if dim(t) >= 4 # Keep at least the minimum Hubbard dimension
    # Make the new M̃ and H̃
    M̃ = Λ̃ * X̃ * X
    H̃ = dag(prime(R,t)) * H * noprime(R)

    # Rearrange/relabel some indices to get back to starting arrangement
    replacetags!(M̃, "Trunc", tags(s))
    replacetags!(H̃, "Trunc", tags(s))
    s_old = getfirst(x -> hastags(x, "Site"), inds(R))
    s_new = getfirst(x -> hastags(x, "Site"), inds(H̃))
    R = R * delta(dag(s_old), s_new)

    # Update
    psi[b] = M̃
    Rs[b] = R
    PH.H[b] = H̃
end

I know that my relabelling/ rearranging of indices might be a bit inelegant, but I am not sure if this is behind my problem. Would love to hear any thoughts you might have on this issue! Thanks so much, and have a great week.

I don’t recall a lot about this method, but is the idea to apply unitary transformations to the MPS during DMRG sweeping to lower the entanglement?

If so, then I think the place I would start to debug such a code would be to check the norm and energy of the MPS immediately before and after applying each unitary. If that does change, then it can’t be the basis transformations applied to the state that is changing the energy. Perhaps you can see in this way that it’s at a particular step that the energy starts to change and go up instead of going down at every step as it should?

1 Like

Hi Miles,

Thanks so much for the reply! I really appreciate it. I think the idea behind this particular implementation is to express the local density matrix on site l as \rho_l^{\sigma_l,\sigma_l'} = U_l \Lambda_l U_l^\dagger, and then to truncate U_l \in \mathbb{R}^{d \times d} by getting rid of the d-d_0 rows of the matrix that correspond to the smallest eigenvalues of \rho_l. This gives us a rotation matrix R^{\sigma_l,\tilde{\sigma_l}} \in \mathbb{R}^{d \times d_0}.

Because I am truncating, I believe the norm is changing(?) Maybe it would be a good idea to renormalize after each truncation step. I can also try applying the rotation without a truncation, which should just be a unitary transform and thus a good way to check that the norm and energy of the MPS do not change with each application.

Thanks again for all the help! Have a great day.

Thanks for explaining. Yes, it does sound very important to account for the norm changing, whether by restoring it back to 1, or by computing the energy in a way that the norm is accounted for. I hope the idea of checking the energy sooner and at each step helps!

Got it to work! The key here was to initialize R properly. See below for functioning code, which can be inserted into the dmrg(...)function right before each two-site optimization step:

if LBO && sw>1

    # Get the left indices of a tensor 
    function get_Linds(O::ITensor, tag::String)
      Rind = getfirst(x -> hastags(x, tag), inds(O))
      return noncommoninds([Rind], inds(O))
    end

    # Put into canonical form
    orthogonalize!(psi,b)
    M = psi[b] 
    H = PH_original.H[b] # we start fresh each time

    # If we are on the first half-sweep, initialize the vector of Rs
    if !hastags(M,"Trunc")#sw==1 && ha==1
      M,R = factorize(M, ortho="right", get_Linds(M, "Site"), tags="Trunc")
      push!(Rs, R)
    else 
      R = Rs[b]
    end

    Linds = get_Linds(M, "Trunc")

    # Perform the first SVD
    X, Λ, Y = svd(M, Linds, lefttags="svd1,u", righttags="svd1,v")
    u = getfirst(x -> hastags(x, "svd1,u"), inds(Λ))
    s = getfirst(x -> hastags(x, "Site"), inds(R))
    t = getfirst(x -> hastags(x, "Trunc"), inds(R))

    # Make and optimize R̃
    R̃ = R * Y * Λ 
    X̃, Λ̃, R = svd(R̃, u, lefttags="svd2,u", righttags="Trunc", maxdim=max_LBO_dim)

    t = getfirst(x -> hastags(x, "Trunc"), inds(R))
    if dim(t) >= min_LBO_dim # Keep at least the minimum Hubbard dimension
      # Make the new M̃ and H̃
      M̃ = Λ̃ * X̃ * X
      H̃ = dag(prime(R)) * H * R

      # Update
      psi[b] = M̃
      Rs[b] = R
      PH.H[b] = H̃
    end
  end
  normalize!(psi)

At the very end, before returning psi, it is critical that you rotate it back into the original basis like so:

if LBO
   for b in 1:length(Rs)
     M = psi[b]
     R = Rs[b]
     psi[b] = M*R # Transform back into original basis 
   end
 end

Glad you figured it out! And thanks for posting an update, which will help this topic to become more helpful to others reading it.