Ground state energy sign is non-deterministically flipping

[MasonProtter]

I've encountered what I believe to be a bug. I have the following system taken from https://doi.org/10.1103/PhysRevLett.124.120503 which is basically a chain where the nodes are Fermions and the links hold spins. The Fermion hopping is modulated by the Sz operator, there's an electric field term, and then an extra term I've added whose role is to energetically enforce a gauge constraint on the Fermions.

using ITensors

"""
     H = ∑ᵢ( -t c(i)' σᶻ(i,i+1) c(i+1) - h σˣ(i, i+1) + ζ (𝟙 - G(i))    

where 

     G(i) = σˣ(i-1,i) (-1)^nf(i) σˣ(i, 1)
"""
function constrained_Hamiltonian_MPO(;L, t, h, ζ=100.0, s)
    ampo = OpSum()
    for n ∈ 1:L
        i  = 2n
        ij = i + 1
        j  = i + 2
        if n < L
            ampo += -2t, "Cdag", i, "Sz", ij, "C", j 
            ampo += -2t, "Cdag", j, "Sz", ij, "C", i
            ampo += -2h, "Sx", ij 
        end
        ampo +=   ζ, "Id", i
        ampo += -4ζ, "Sx", ij-2, "Sx", ij
        ampo +=  4ζ, "Sx", ij-2, "Cdag", i, "C", i, "Sx", ij

    end
    MPO(ampo, s)
end

function constrained_ground_via_dmrg(;L, t=1.0, h, ζ=100.0, outputlevel=0)
    s = siteinds(n -> isodd(n) ? "S=1/2" : "Fermion", 2L+1)
    H = constrained_Hamiltonian_MPO(;L, t, h, ζ, s)

    sweeps = Sweeps(10)
    maxdim!(sweeps,10,10,20,40,80,100,140,180,200)
    cutoff!(sweeps,1e-8)

    ψ0 = randomMPS(s, linkdims=4)

    λ, ϕ = dmrg(H, ψ0, sweeps; outputlevel)
end

I've studied this system (for small sizes) using exact diagonalization, but comparing to the DMRG, I'm finding that ITensors.jl is non-deterministically flipping the sign of the returned energy! Here's an example:

julia> map(1:5) do _
           constrained_ground_via_dmrg(L=8, h=1.0)[1]
       end
5-element Vector{Float64}:
 -6.99999999999987
 -6.9999999999998375
  7.0000000000002185
 -6.99999999999985
 -6.999999999999846

and here is is again for another choice of parameters

julia> map(1:5) do _
           constrained_ground_via_dmrg(L=5, h=-5.2)[1]
       end
5-element Vector{Float64}:
 -20.799999999999926
  20.800000000000065
 -20.799999999999887
 -20.79999999999991
 -20.799999999999915

The correct answer according to my exact diagonalizations in both cases is the negative one. Any ideas on how to fix this? I'm on ITensors.jl v0.2.6 and I've tried Julia v1.7.0-beta4 and v1.6.1 on separate computers.

[mtfishman]

There's something a bit funny about the definition of that Hamiltonian. I'm not sure how those final terms in the OpSum proportional to ζ correspond to the term you wrote in the docstring of the function. For example, when you write ampo += ζ, "Id", i, that will add a new identity term proportional to ζ for every site in the lattice, so it seems like you would be adding on a term L ζ 𝟙 to your overall Hamiltonian. Also, when I run the code and look at the bond dimensions, they are all 1. It makes me think that with ζ = 100, somehow you are just finding product eigenestates of the identity operator or something a bit contrived like that. When I make ζ smaller I don't notice spurious positive energies.

[MasonProtter]

The L ζ 𝟙 is just so that doing - ζ ∑ᵢ G(i) doesn't offset the energy. (I want to be in the G(i) = +1 eigenstates), however this made me realize that there is a different error, I wrote

        ampo += -4ζ, "Sx", ij-2, "Sx", ij
        ampo +=  4ζ, "Sx", ij-2, "Cdag", i, "C", i, "Sx", ij

but that should actually be

        ampo += -4ζ, "Sx", ij-2, "C", i, "Cdag", i, "Sx", ij
        ampo +=  4ζ, "Sx", ij-2, "Cdag", i, "C", i, "Sx", ij

in order to correctly get σˣ(i-1,i) (-1)^nf(i) σˣ(i, 1). I made this error on the ED side as well, things seem to be better now.

However, I still find this non-deterministic sign flipping quite suspicious! It seems to indicate some sort of problem if the algorithm is sometimes giving an energy that's positive and sometimes negative, right?

[emstoudenmire]

Hi Mason, based on Matt’s analysis it doesn’t look like this currently points to a problem with our DMRG code. The reason is that DMRG is (1) never really guaranteed to find the ground state, and in fact can many times be sensitive to choice of initial state and (2) especially if the ground state happens to have a bond dimension of 1, that could make DMRG struggle or totally fail.

The reason for (2) is that DMRG is doing two optimizations in tandem: one to find a better basis in which to represent the ground state and another to minimize the energy in a given basis. If the bond dimension is 1, DMRG is basically stuck with a very small basis. If such a small basis is wrong in any way, then DMRG will not be able to improve the energy, nor improve the basis (basically because there are no other basis states to ‘discover’ or mix together to make improvements, and for other subtler reasons).

So it seems likely that your Hamiltonian may be a kind of “worst case” for a DMRG approach to minimization and/or may still have an issue with its definition.

It would be good to analyze the properties of the ground state and Hamiltonian (i.e. trying other values of zeta as Matt mentioned) some more in a way that doesn’t use DMRG to see if all the definitions (e.g. OpSum inputs) are what they should be and the ground states have ‘typical’ enough properties for DMRG to help. One way would be to literally merge all your MPO tensors together on a small system, then use eigen to diagonalize the resulting H tensor and study the results that way.

Hope that’s helpful.

Miles

[emstoudenmire]

Going to close this for now. Please reopen if you determine that there is a bug versus DMRG getting 'stuck' in a local minimum and/or an issue with the Hamiltonian. Thanks for reporting this -