Problem using TEBD algorithm to represent Gibbs states

General Discussion
,
1

Hi,

I am trying to get thermal states \rho_\beta = e^{-\beta H} mean values of local observables for long-range interacting system. More concretely the Long Range Spin chain:

H_{\text{LR, Ising}}= J \sum_{i \neq j}^N \frac{\sigma_i^x \sigma_{j}^x}{|i-j|^{\alpha}} + h \sum_{i=1}^N \sigma_i^z

Notice that the smaller the \alpha, the “least” local the model is, and thereby MPs and MPOs representations are less efficient as the entanglement of the system is greater.

We construct the time evolution operator of the LR-Ising model explicitly in the short-range case(Radware Bot Manager Captcha) and then applying a super-extensive amount of swap gates to exactly simulate the power-law interactions in a finite system size. This yields a Trotterized representation of the time evolution operator U(\delta t) = \mathrm{e}^{-iH \delta t} in form of an MPO. We concatenate the Trotterized time evoltion to yield the desired final imaginary time or \beta U(t) \simeq \prod_{i=1}^{n_t} U(\delta t), with t = n_t \cdot \delta t.

When I perform the simulations for \alpha \leq 1.5, despite the bond dimension not to collapse the evolution stops and we have the following error:

ERROR: LoadError: BoundsError: attempt to access 0-element Vector{Pair{QN, Int64}} at index [1]
Stacktrace:
[1] throw_boundserror(A::Vector{Pair{QN, Int64}}, I::Tuple{Int64})
@ Base ./essentials.jl:14
[2] getindex
@ ./essentials.jl:916 [inlined]
[3] combineblocks(qns::Vector{Pair{QN, Int64}})
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/qn/qnindex.jl:413
[4] combineblocks
@ ~/.julia/packages/ITensors/4M2ep/src/qn/qnindex.jl:475 [inlined]
[5] combiner(inds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}; dir::Nothing, tags::String)
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/qn/qnitensor.jl:444
[6] combiner(::Index{Vector{Pair{QN, Int64}}}, ::Vararg{Index{Vector{Pair{QN, Int64}}}}; kwargs::@Kwargs{dir::Nothing})
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/tensor_operations/itensor_combiner.jl:7
[7] svd(A::ITensor, Linds::Vector{Index{Vector{Pair{QN, Int64}}}}; leftdir::Nothing, rightdir::Nothing, lefttags::TagSet, righttags::TagSet, mindim::Int64, maxdim::Int64, cutoff::Float64, alg::Nothing, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing, min_blockdim::Nothing, utags::TagSet, vtags::TagSet)
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/tensor_operations/matrix_decomposition.jl:153
[8] svd
@ ~/.julia/packages/ITensors/4M2ep/src/tensor_operations/matrix_decomposition.jl:110 [inlined]
[9] factorize_svd(A::ITensor, Linds::Vector{Index{Vector{Pair{QN, Int64}}}}; singular_values!::Nothing, ortho::String, alg::Nothing, dir::Nothing, mindim::Int64, maxdim::Int64, cutoff::Float64, tags::TagSet, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing, min_blockdim::Nothing)
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/tensor_operations/matrix_decomposition.jl:614
[10] factorize(A::ITensor, Linds::Vector{Index{Vector{Pair{QN, Int64}}}}; mindim::Int64, maxdim::Int64, cutoff::Float64, ortho::String, tags::TagSet, plev::Nothing, which_decomp::Nothing, eigen_perturbation::Nothing, svd_alg::Nothing, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing, min_blockdim::Nothing, singular_values!::Nothing, dir::Nothing)
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/tensor_operations/matrix_decomposition.jl:808
[11] contract(::NDTensors.BackendSelection.Algorithm{:zipup, @NamedTuple{}}, A::MPO, B::MPO; cutoff::Float64, maxdim::Int64, mindim::Int64, kwargs::@Kwargs{})
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/4M2ep/src/ITensorMPS/mpo.jl:845
[12] contract
@ ~/.julia/packages/ITensors/4M2ep/src/ITensorMPS/mpo.jl:815 [inlined]
[13] contract(A::MPO, B::MPO; alg::String, kwargs::@Kwargs{cutoff::Float64, maxdim::Int64})
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/4M2ep/src/ITensorMPS/mpo.jl:808
[14] contract
@ ~/.julia/packages/ITensors/4M2ep/src/ITensorMPS/mpo.jl:807 [inlined]
[15] product(A::MPO, B::MPO; kwargs::@Kwargs{cutoff::Float64, maxdim::Int64})
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/4M2ep/src/ITensorMPS/mpo.jl:885
[16] top-level scope
@ ~/jssegovia/Codes/CorrelationOfDecayLR/CorrelationOfDecayLR/GibbsExpectationValue_tosend.jl:87
in expression starting at /lustre/home/ift/jssegovia/jssegovia/Codes/CorrelationOfDecayLR/CorrelationOfDecayLR/GibbsExpectationValue_tosend.jl:83

The code is the following, and I’ve used the parameters: \delta t=0.01, SVDcutoff=1e-18, and the evolution stopped at \beta=2.85 when I was aiming \beta=3. The problem gets worse when decreasing \alpha:

The code is something like this:

for (i, t) in enumerate(β_values)
    contraction_U[:] = apply(U_δt, contraction_U; cutoff=SVDcutoff, maxdim=SVDmaxdim) #here U_δt is the infinitesimal imaginary time evolution MPO, and  contraction_U is initialized as the Identity MPO

    trace_values[i] = tr(contraction_U)
    A_values[i]= tr(apply( contraction_U, A; cutoff=SVDcutoff, maxdim=SVDmaxdim))*(1/trace_values[i])

    maxdim = maxlinkdim(contraction_U)

    open(csv_filename, "a") do file
        writedlm(file, [t trace_values[i]  A_values[i]  maxdim], ',')
    end

    @printf("After t=%g\tmaxlinkdim=%d\ttimestep %d out of %d", t, maxdim, i, ntimes)
    println(" runtime: $(Dates.canonicalize(Dates.now() - t1))")
end

I don’t have a good intuition of what is going on, If someone has it would be really helpful!

2

It would be difficult to debug this without a fully runnable code.

Can you share an example of MPOs you are producting together that leads to that error (say by saving them with HDF5 or JLD2) along with the values of cutoff and maxdim you passed to apply?

3

Sure! The functions I use for the MPO representation are :

function build_expH_ising_murg(sites, JXX::Real, hz::Real, dt::Number)

  

  # For real dt this does REAL time evolution
  # I should have already taken into account both the - sign in exp(-iHt)
  # and the overall minus in Ising H= -(JXX+hZ)

  # dt = JXX * dt 
  # JXX= JXX * dt                                                             ## why?

  cosg = cos(hz * dt * 0.5)
  sing = sin(hz * dt * 0.5)

  N = length(sites)
  U_t = MPO(N)

  link_dimension = 2

  # linkindices = hasqns(sites) ? [Index(link_dimension, "Link,l=$(n-1)") for n in 1:N+1] : [Index(link_dimension, "Link,l=$(n-1)") for n in 1:N+1]
  linkindices = [Index(link_dimension, "Link,l=$(n-1)") for n in 1:N+1]

  for n in 1:N
    # siteindex s
    s = sites[n]

    # left link index ll with daggered QN conserving direction (if applicable)
    ll = dag(linkindices[n])
    # right link index rl
    rl = linkindices[n+1]

    combz = (1 - 2 * sing^2) * op(sites, "Id", n) + im * 2 * sing * cosg * op(sites, "Z", n)
    X = op(sites, "X", n)

    if n == 1
      #U_t[n] = ITensor(ComplexF64, dag(s), s', dag(rl))
      U_t[n] = onehot(rl => 1) * sqrt(cos(JXX*dt)) * combz
      U_t[n] += onehot(rl => 2) * sqrt(im * sin(JXX*dt)) * X
    elseif n == N
      #U_t[n] = ITensor(ComplexF64, ll, dag(s), s')
      U_t[n] = onehot(ll => 1) * sqrt(cos(JXX*dt)) * combz
      U_t[n] += onehot(ll => 2) * sqrt(im * sin(JXX*dt)) * X

    else
      #U_t[n] = ITensor(ComplexF64, ll, dag(s), s', dag(rl))

      U_t[n] = onehot(ll => 1, rl => 1) * cos(JXX*dt) * combz
      U_t[n] += onehot(ll => 1, rl => 2) * sqrt(im * sin(JXX*dt)) * sqrt(cos(JXX*dt)) * X
      U_t[n] += onehot(ll => 2, rl => 1) * sqrt(im * sin(JXX*dt)) * sqrt(cos(JXX*dt)) * X
      U_t[n] += onehot(ll => 2, rl => 2) * im * sin(JXX*dt) * combz
    end
  end
  return U_t
end



function XX_NN_layer(sites, J::Float64)
  N = length(sites)
  return [op("Rxx", sites[j], sites[j+1], ϕ=-J) for j in 1:N-1 ]
end

function XX_powerlaw_layer(sites, α::Float64, J::Number)
  N = length(sites)
  layer = Vector{ITensor}(undef,0)
  # number_of_terms = 3 * N * (N - 1)/2 + 3*(N-1)
  for i in 1:N-1
    for j in i+1:N-1
      # verify factor!
      push!(layer, op("Rxx", sites[j - 1], sites[j], ϕ=-J * (1.0 / ((j - i)^α))) )
      push!(layer, op("Swap", sites[j - 1], sites[j]))
    end
    push!(layer, op("Rxx", sites[N - 1], sites[N], ϕ=-J * (1.0 / ((N - i)^α))))
    push!(layer, [op("Swap", sites[j - 1], sites[j]) for j in N-1:-1:i+1]...)
  end
  return layer
end

function XX_NTNT_layer(sites, J::Float64)
  N = length(sites)
  layer = Vector{ITensor}(undef,0)
  # number_of_terms = 3 * N * (N - 1)/2 + 3*(N-1)
  for i in 1:N-2
    push!(layer, op("Rxx", sites[i], sites[i+1], ϕ=-J))
    push!(layer, op("Swap", sites[i], sites[i+1]))
    push!(layer, op("Rxx", sites[i+1], sites[i+2], ϕ=-J))
    push!(layer, op("Swap", sites[i+1], sites[i]))
  end
  push!(layer, op("Rxx", sites[N-1], sites[N], ϕ=-J))
  return layer
end


function Rz_rotation_layer(sites, ϕ::Number)
  N = length(sites)
  # verify the factor!!!
  return MPO([op("Rz", sites[j], θ=-2ϕ) for j in 1:N], 0, N+1)
end


function Rx_rotation_layer(sites, ϕ::Number)
  N = length(sites)
  # verify the factor!!!
  return MPO([op("Rx", sites[j], θ=-2ϕ) for j in 1:N], 0, N+1)
end

function timeEvo_Cirac_LR_TFI(sites, α::Float64, J::Float64,hp::Float64,δtn::Number)
  Uz = Rz_rotation_layer(sites,hp*δtn*0.5)
  Uxx = XX_powerlaw_layer(sites,α,J*δtn)

  Uz_Uxx = mapreduce(
    e -> apply(e, Uz; maxdim=SVDmaxdim,cutoff=SVDcutoff),
    (previousResult, newElement) -> apply(newElement, previousResult; maxdim=SVDmaxdim, cutoff=SVDcutoff),
  Uxx
  )

  return mapreduce(e -> apply(e, Uz_Uxx; maxdim=SVDmaxdim,cutoff=SVDcutoff),
  (previousResult, newElement) -> apply(newElement, previousResult; maxdim=SVDmaxdim,cutoff=SVDcutoff),
  Uz)
end

function timeEvo_Cirac_LR_TFI_corr(sites, α::Float64, J::Float64,hp::Float64,δtn::Number)
  Uz = Rz_rotation_layer(sites,hp*δtn*0.5)
  Uxx = XX_powerlaw_layer(sites,α,J*δtn)
  return mapreduce(e -> apply(Uz,e; maxdim=SVDmaxdim,cutoff=SVDcutoff), e-> apply(Uxx,e; maxdim=SVDmaxdim,cutoff=SVDcutoff), e-> apply(Uz,e; maxdim=SVDmaxdim,cutoff=SVDcutoff))

  # Uz_Uxx = mapreduce(
  #   e -> apply(Uz,e; maxdim=SVDmaxdim,cutoff=SVDcutoff),
  #   (previousResult, newElement) -> apply(newElement, previousResult; maxdim=SVDmaxdim, cutoff=SVDcutoff),
  # Uxx
  # )

  # return mapreduce(e -> apply(Uz_Uxx,e; maxdim=SVDmaxdim,cutoff=SVDcutoff),
  # (previousResult, newElement) -> apply(newElement, previousResult; maxdim=SVDmaxdim,cutoff=SVDcutoff),
  # Uz)
end




function timeEvo_Cirac_LR_TFI_NewVersion(sites, α::Float64, J::Float64,hp::Float64,δtn::Number,SVDcutoff::Float64,SVDmaxdim::Int64)
  OperatorMPO = Rz_rotation_layer(sites,hp*δtn*0.5)
  Uz = Rz_rotation_layer(sites,hp*δtn*0.5)
  Uxx = XX_powerlaw_layer(sites,α,J*δtn)
  
  # OperatorMPO[:] =  apply(Uz, OperatorMPO; cutoff=SVDcutoff, maxdim=SVDmaxdim)
  OperatorMPO[:] =  apply(Uxx, OperatorMPO; cutoff=SVDcutoff, maxdim=SVDmaxdim)
  OperatorMPO[:] =  apply(Uz, OperatorMPO; cutoff=SVDcutoff, maxdim=SVDmaxdim)

  return OperatorMPO
end

and the operator is defined as:

const U_δt = timeEvo_Cirac_LR_TFI_NewVersion(sites, α, J, hp, δβ, SVDcutoff, SVDmaxdimOperator)
const contraction_U = MPO(map(n -> op(sites, "Id", n), 1:length(sites)))

Then, the parameters you asked for: const SVDmaxdimOperator = 30, const SVDcutoff = 1e-18, const SVDmaxdim=1800, const \delta \beta = - 0.01* im, const J=1, const hp=0.5. For \alpha\leq 1.5 system size above 55 I start having this problem, and the evolution does not finish.

pd: I am creating a hdf5 file with the contraction, I did not have it before as I simply loaded the data in a csv to save some space

Thank you very much for your help!