[ITensors] [BUG] An unexpected excited state energy for a 4 site S=1/2 Heisenberg model from DMRG

[hczhai]

Description of bug When I try the following script for getting the low energy spectrum of a 4-site S=1/2 Heisenberg model with custom topology, the DMRG with projection generates an (excited) energy level which is not part of the exact diagonalization spectrum.

How to reproduce The following script does ED and then DMRG for this model:

    using ITensors

    N = 4
    sites = siteinds("S=1/2",N; conserve_qns=true)

    fe4 = [
        (1, 2, 1.17),
        (3, 4, 1.17),
        (1, 4, 1.00),
        (2, 3, 1.00),
        (1, 3, 1.00),
        (2, 4, 1.00),
    ]

    global ampo = OpSum()
    for j=1:6
      global ampo += 1.0 * fe4[j][3],"Sz",fe4[j][1],"Sz",fe4[j][2]
      global ampo += 1/2 * fe4[j][3],"S+",fe4[j][1],"S-",fe4[j][2]
      global ampo += 1/2 * fe4[j][3],"S-",fe4[j][1],"S+",fe4[j][2]
    end

    H = MPO(ampo,sites)

    # exact diag

    hfull = prod(H)
    d, u = eigen(hfull, ishermitian=true)

    @show d

    # dmrg

    state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
    psi0 = randomMPS(sites, state, linkdims=5)

    sweeps = Sweeps(5)
    setmaxdim!(sweeps, 200,500,2000)
    setcutoff!(sweeps, 0.0)
    setnoise!(sweeps, 1E-2, 1E-3, 1E-4, 0)

    energy, psi = dmrg(H,psi0, sweeps)
    energy, psi2 = dmrg(H,[psi],psi0, sweeps)
    energy, psi3 = dmrg(H,[psi,psi2],psi0, sweeps)
    energy, psi4 = dmrg(H,[psi,psi2,psi3],psi0, sweeps)
    energy, psi5 = dmrg(H,[psi,psi2,psi3,psi4],psi0, sweeps)

Expected behavior

The DMRG should reproduce some of the ED energy levels in the singlet (SZ = 0) sector (-1.755, -1.415, -0.585, -0.415, ...)

Actual behavior

The DMRG produces one additional energy level -0.755 which is not part of the ED spectrum.

Version information

• output from versioninfo() surrounded by backticks (`) Julia Version 1.7.2 Commit bf53498635 (2022-02-06 15:21 UTC) Platform Info: OS: macOS (x86_64-apple-darwin19.5.0) CPU: Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz WORD_SIZE: 64 LIBM: libopenlibm LLVM: libLLVM-12.0.1 (ORCJIT, icelake-client)`
• output from ] status ITensors surrounded by backticks (`)Status '~/.julia/environments/v1.7/Project.toml\' [9136182c] ITensors v0.2.15`

[emstoudenmire]

Thanks for the report. I found the issue, however. The excited state DMRG method is not guaranteed to return excited states, but rather uses a Langrage multiplier or "weight" (whose size can be adjusted) to "penalize" overlap of the new states being computed with the previous ones. Sometimes this strategy can fail, whether because of the initial state used or because of the value of the penalty "weight" value (e.g. if it is too small or too large relative to the energy gaps).

Here's what I find when printing out the overlaps of the states obtained from the above code:

inner(psi, psi2) = -1.5097548993758442e-15 inner(psi, psi3) = -1.0 inner(psi, psi4) = 1.942890293094024e-15 inner(psi, psi5) = 1.0190786078592668e-15 inner(psi2, psi3) = -4.135492081000051e-15 inner(psi2, psi4) = 3.885780586188048e-16 inner(psi2, psi5) = -0.6650294984203701 inner(psi3, psi4) = 1.0547118733938987e-15 inner(psi3, psi5) = 2.7165764727593182e-15 inner(psi4, psi5) = 1.6514567491299204e-15

So you can see that psi3 is actually a second copy of psi, and that psi2 and psi5 also have significant overlap.

[hczhai]

Thanks for the quick reply! So one can check the overlaps between states to figure out which of them are the true ground/excited states then.

[emstoudenmire]

Yes, you have the right idea, and then one can sort them by energy and discard any that aren't orthogonal. Another thing you can do is to make a "mini matrix" out of the matrix elements of the psi's with H like M[i,j] = inner(psi[i],H,psi[j]) and diagonalize this matrix to get even more precision. It's like the idea of a Krylov space.

The main risk of these methods, though, is that it could miss an excited state, like skipping over one. So it's important to be careful about that.

[hczhai]

Thanks! I agree that missing states can be more risky.

Indeed, by setting a larger "weight", I can now successfully get the full correct spectrum using DMRG. Also I guess using different random initial guesses every time may help reducing the chance of missing states.

So now the script looks like this:

    psi0 = randomMPS(sites, state, linkdims=5)
    energy, psi = dmrg(H,psi0, sweeps, weight=4)
    psi0 = randomMPS(sites, state, linkdims=5)
    energy2, psi2 = dmrg(H,[psi],psi0, sweeps, weight=4)
    psi0 = randomMPS(sites, state, linkdims=5)
    energy3, psi3 = dmrg(H,[psi,psi2],psi0, sweeps, weight=4)
    psi0 = randomMPS(sites, state, linkdims=5)
    energy4, psi4 = dmrg(H,[psi,psi2,psi3],psi0, sweeps, weight=4)
    psi0 = randomMPS(sites, state, linkdims=5)
    energy5, psi5 = dmrg(H,[psi,psi2,psi3,psi4],psi0, sweeps, weight=4)
    psi0 = randomMPS(sites, state, linkdims=5)
    energy6, psi6 = dmrg(H,[psi,psi2,psi3,psi4,psi5],psi0, sweeps, weight=4)

And this generates 6 correct states without any repetition/missing states, in the desired order, which is great. Thanks again for your help!

[emstoudenmire]

Those are great ideas to improve the code! Glad it worked much better with those.

[emstoudenmire]

I'll update our excited state code example to use this idea.

[mtfishman]

@hczhai glad this got settled, thanks @emstoudenmire for the help. Good suggestion to form the matrices h[i,j] = inner(psi[i],H,psi[j]) and n[i,j] = inner(psi[i],psi[j]). We should mention tricks like that in the docs section on computing excited states (https://itensor.github.io/ITensors.jl/stable/examples/DMRG.html#Compute-excited-states-with-DMRG), along with setting a higher weight, using random starting states, sweeping back through the excited states, etc.

[emstoudenmire]

Good idea to add those. I actually just did a quick PR to switch to using randomMPS. We could definitely add the other ideas too.

[mtfishman]

Great, thanks!

In the past I've wondered how well it works to diagonalize the h[i,j] = inner(psi[i],H,psi[j]) and n[i,j] = inner(psi[i],psi[j]) matrices and then rotate the MPS states into that basis, have you tried that?

[emstoudenmire]

Just tried it real quickly and seems to work quite well. Julia has a convenient generalized eigen called via eigen(A,B). Here is the code I made for the case of three states:

  psis = [psi0,psi1,psi2]

  h = zeros(3,3)
  n = zeros(3,3)
  for i=1:3,j=1:3
    h[i,j] = inner(psis[i],H,psis[j])
    n[i,j] = inner(psis[i],psis[j])
  end

  display(h)
  display(n)
  vals,vecs = eigen(h,n)
  @show vals

  @show vecs[:,1]
  @show vecs[:,2]
  @show vecs[:,3]

  v0 = vecs[:,1]
  phi0 = sum([v0[1]*psi0,v0[2]*psi1,v0[3]*psi2];cutoff=1E-9)
  normalize!(phi0)

  v1 = vecs[:,2]
  phi1 = sum([v1[1]*psi0,v1[2]*psi1,v1[3]*psi2];cutoff=1E-9)
  normalize!(phi1)

  v2 = vecs[:,3]
  phi2 = sum([v2[1]*psi0,v2[2]*psi1,v2[3]*psi2];cutoff=1E-9)
  normalize!(phi2)

  psis = [phi0,phi1,phi2]

  h = zeros(3,3)
  n = zeros(3,3)
  for i=1:3,j=1:3
    h[i,j] = inner(psis[i],H,psis[j])
    n[i,j] = inner(psis[i],psis[j])
  end

  display(h)
  display(n)

And some sample output (which I've edited a little bit):

Before:
h = 
3×3 Matrix{Float64}:
 -81.1917       -0.00758735   -0.00137357
  -0.00758735  -75.1239       -0.36207
  -0.00137357   -0.36207     -74.9461
n = 
3×3 Matrix{Float64}:
 1.0         9.34533e-5  1.69203e-5
 9.34533e-5  1.0         0.0052778
 1.69203e-5  0.0052778   1.0

vals = [-81.19170871435301, -75.13030162112851, -74.93998903980001]
vecs[:, 1] = [-1.0, 4.560142526693715e-8, 3.424441094794933e-8]
vecs[:, 2] = [9.024964040918464e-5, -1.0, 0.18701242225693476]
vecs[:, 3] = [3.3877974252682016e-5, -0.18191417401455068, -1.0]

After:
h = 
3×3 Matrix{Float64}:
 -81.1917        -2.96081e-5   -0.000444844
  -2.96081e-5   -75.1303        2.86782e-7
  -0.000444844    2.86782e-7  -74.94
n =
3×3 Matrix{Float64}:
 1.0          3.67584e-7   5.4819e-6
 3.67584e-7   1.0         -3.15353e-9
 5.4819e-6   -3.15353e-9   1.0

You can see the first time the h and n matrices are printed out there are significant off-diagonal elements in n and that some of the energies are higher. Then afterward everything is better.

[mtfishman]

Wow, cool! The code is really simple too, we should definitely put something like that in the examples/docs.

In fact if it seems to work well in general we could package that up in a function...