ED vs. DMRG overlap != 1.

[dolefeast]

Hi, I’ve just started using the ITensor library, and the first thing I wanted to do is to get the ground state of the TFIM. This has already been done, I know, but I cannot find my specific error.

What I wanted to do is to solve the problem by “vanilla” methods (just with the linear algebra package) and by using the DMRG method from ITensor, and compare them. My strategy is to construct the Hamiltonian as a matrix (H_ED) and as an MPO (H_MPO), find the eigenstates of both using, respectively, LinearAlgebra.eigen and dmrg.
Finally, compare them through their inner product, by transforming the groundstate of the exact diagonalization into a ITensorMPS.MPS and taking inner(groundstate_ED, groundstate_MPS), which I expect to be (very close to) 1.

The MWE:

using ITensors
using MKL
using ITensorMPS
using Plots
using LinearAlgebra
using LaTeXStrings

function Ising_MPO(L::Int, g; sites)
    # Operator terms
    os = OpSum()
    for i in 1:L-1
        os += -1, "Sz", i, "Sz", i+1
        os += -g, "Sx", i
    end
    os += -g, "Sx", L
    return MPO(os, sites)
end

function Ising_ED(L::Int64, g::Float64; sparsify::Bool=false)

    id = [ 1. 0 ; 0 1.]
    sigma_x = [ 0 0.5 ; 0.5 0 ]
    sigma_z = [ 0.5 0 ; 0 -0.5 ]
    H = zeros(Int, 2^L, 2^L) 

    # The first term of the NN Hamiltonian is σᶻ⊗σᶻ⊗I⊗... 
    # the rest are just cyclic permutations
    H_NN_ops = fill(id, L)
    H_NN_ops[1] = sigma_z
    H_NN_ops[2] = sigma_z

    # The transversal field hamiltonian
    H_trans_ops = fill(id, L)
    H_trans_ops[1] = sigma_x

    
    for i in 1:L-1
        # Tensor product of the operators in H_NN_ops
        H -= foldl(kron, H_NN_ops) # Similar to "functools.reduce"
        H_NN_ops = circshift(H_NN_ops, 1) # cycle the ops
    end

    for i in 1:L
        H -= g * foldl(kron, H_trans_ops)
        H_trans_ops = circshift(H_trans_ops, 1)
    end
    
    H
end

function MPS_ED_overlap(L::Int, g::Float64; sparsify=false)
    sites = siteinds("S=1/2", L)

    H_ED = Ising_ED(L, g, sparsify=sparsify) # The Hamiltonian for exact diag.
    H_MPO = Ising_MPO(L, g, sites=sites) # Hamiltonian for DMRG

    #Exact diagonalization

    energy_ED, groundstate_ED = eigen(H_ED)
    display(energy_ED)
    energy_ED = energy_ED[1]
    groundstate_ED = groundstate_ED[:,1]

    groundstate_ED_MPS = MPS(groundstate_ED, sites)

    # DMRG
    state_i  = [ isodd(n) ? "Up" : "Dn" for n=1:L ]
    psi0_i = MPS(sites, state_i)
    sweeps = Sweeps(10);

    #println("L = $L, g = $g")
    energy_MPS, psi0 = dmrg(H_MPO, psi0_i, sweeps, outputlevel=0);
    
    norm_ED2 = inner(groundstate_ED_MPS, groundstate_ED_MPS)
    norm_MPS = inner(psi0, psi0)
    
    # In case any of the resulting groundstates is not normalized
    overlap = inner(groundstate_ED_MPS, psi0) / sqrt(norm_ED2*norm_MPS)
    return energy_ED, energy_MPS, overlap
end

When studying this for different chain lengths and external fields using the following code

p = plot(layout=2)

Ls = 2:2:10
gs = 10 .^ range(-4, stop=0., length=10)
# Checking the overlap between both eigenstates

for L in Ls
    overlaps = zeros(length(gs))
    energies = zeros(length(gs))
    for (i, g) in enumerate(gs)
        energy_ED, energy_MPS, overlap = MPS_ED_overlap(L, g)
        overlaps[i] = abs(overlap)^2
        energies[i] = (energy_ED-energy_MPS)/energy_ED
    end
    plot!(p, gs, overlaps, xscale=:log10, marker=:circle,
        label="L=$L", xlab="g", ylab=L"$\vert\langle \psi_{ED} \vert \psi_{DMRG}\rangle\vert^2$",
        subplot=1)
    plot!(p, gs, energies, xscale=:log10, marker=:circle,
    label="L=$L", xlab="g", ylab=L"$\frac{E_{ED} - E_{MPS}}{E_{ED}}$",
    subplot=2)
  
  
end

plot!()

the following plot is returned:

Left panel is the overlap of both groundstates as a function of the external field g and right panel is the relative energy difference also as a function of the external transverse field.

[miles]

The TFIM has a ground state degeneracy (related to spontaneous symmetry breaking) for smaller values of g, up to some critical value. So most likely DMRG is just finding one of the two ground states at random while ED is likely making a superposition of both ground states. If I’m right about that, neither answer is really right or wrong , but just different ways of conceptualizing the physics.

[dolefeast]

Thanks for the speedy reply!

It does seem to be due to the degeneracy of the ground state.

[dolefeast]

The fix was to calculate also the projection of the DMRG ground state with the second eigenvector given by the eigen decomposition, to account for the degeneracy of the Hamltonian for small `g.