Hi everyone,
I struggle to understand some differences when using MPS in iTensor and generic iTN in iTensorNetwork. I am running Julia Version 1.9.4, ITensors v0.3.54, and ITensorNetworks v0.3.10.
I have implemented the dynamics of a generic lattice spin model using MPS and iTN and got different results. The ITensorNetworks implementation is based on the heavy_hex_ising_real_tebd.jl example that was included on Github in ITensorNetworks v0.3.10. The iTensorNetwork is regauged to the Vidal Gauge at every time step.
In the figure on the left, you can see the magnetization in time produced by the MPS and iTN implementations, and the overlap between the wavefunctions of the MPS and iTN. In the right figure, you can see the distance to the Vidal Gauge for the iTN.
The MPS correctly captures the dynamics, while the iTN does not. If I increase the cutoff to e.g. 1e-2 the behaviour of the MPS approaches the one of the iTN. However, given a maxdim=32, cutoff=1e-14 I don’t understand why the iTN does not approach the correct dynamics of the MPS.
Do you have any hint on why this is the case? Does it have something to do with how the regauging to Vidal gauge is implemented below?
using LinearAlgebra
using ITensors
using ITensorNetworks
using ITensorNetworks:
initialize_bond_tensors,
vidal_itn_canonicalness,
vidal_to_symmetric_gauge,
vidal_gauge,
approx_network_region,
optimal_contraction_sequence,
contract,
symmetric_to_vidal_gauge,
norm_network
using ITensorUnicodePlots
using UnicodePlots
using Graphs
using NamedGraphs
using Plots
######### Util functions
"""Take the expectation value of o on an ITN using belief propagation"""
function expect_state_SBP(o::ITensor, ψ::AbstractITensorNetwork, ψψ::AbstractITensorNetwork, mts::DataGraph)
Oψ = apply(o, ψ; cutoff=1e-16)
ψ = copy(ψ)
s = siteinds(ψ)
vs = vertices(s)[findall(i -> (length(commoninds(s[i], inds(o))) != 0), vertices(s))]
vs_braket = [(v, 1) for v in vs]
numerator_network = approx_network_region(
ψψ, mts, vs_braket; verts_tn=ITensorNetwork(ITensor[Oψ[v] for v in vs])
)
denominator_network = approx_network_region(ψψ, mts, vs_braket)
num_seq = contraction_sequence(numerator_network; alg="optimal")
den_seq = contraction_sequence(numerator_network; alg="optimal")
return contract(numerator_network; sequence=num_seq)[] /
contract(denominator_network; sequence=den_seq)[]
end
function tensornetwork_to_array(psi_tn::ITensorNetwork)
"""
converts a tensor network state with n vertices (external indices) to a array of dimension (1,2^n)
"""
if typeof(psi_tn)==ITensorNetwork{Int64}
phi = psi_tn[1]
for j in 2:length(vertices(psi_tn))
phi = phi * psi_tn[j]
end
state_arr = Array{ComplexF64}(phi,reverse([indx for indx in inds(phi)]))
state_vector = reshape(state_arr,(1,prod(size(state_arr))))
end
if typeof(psi_tn)==ITensorNetwork{Tuple{Int64, Int64}}
key_vals = [(i,j) for (i,j) in keys(vertex_data(psi_tn))]
phi = psi_tn[key_vals[1]]
for j in 2:length(key_vals)
phi = phi * psi_tn[key_vals[j]]
end
state_vector = Array{ComplexF64}(phi)
end
return state_vector
end
function mps_to_array(psi::MPS)
"""
converts a matrix product state MPS with n indices to a array of dimension (1,2^n)
"""
phi = psi[1]
for j in 2:length(psi)
phi = phi * psi[j]
end
state_arr = Array{ComplexF64}(phi,reverse([indx for indx in inds(phi)]))
state_vector = reshape(state_arr,(1,prod(size(state_arr))))
return state_vector
end
######### Parameters
n = 3
Jc = -1.0
Gamma = 0.5
ttot = 3.0
dt = 0.1
ttimes = 0:dt:ttot
g = Graphs.complete_graph(n)
g = NamedGraph(g)
g_edges = edges(g)
g_vertices = vertices(g)
J = Jc*ones(ne(g))
######### Initialistation
sites = siteinds("S=1/2",g)
sites_ren=[only(sites[vidx]) for vidx in vertices(sites)]
states = ["Up", "Up", "Up"]
psi_init_mps = MPS(sites_ren, states)
psi_init_tn = ITensorNetwork(sites, v -> states[v])
bond_tensors = initialize_bond_tensors(psi_init_tn)
######### Gate Definition
gates_X = ITensor[]
for v in g_vertices
sj = sites[v]
Xj = op("X", sj)
Gj = exp(-im * dt / 2 * Gamma * Xj)
push!(gates_X, Gj)
#println("Gate X: v ", sj)
end
gates_ZZ = ITensor[]
# do 2nd order trotter H_ZZ with 1/2 factor
for (j,e) in enumerate(g_edges)
vi = e.src
vj = e.dst
si = sites[vi]
sj = sites[vj]
ZZj = J[j] * op("Z", si) * op("Z", sj)
Gj = exp(-im * dt / 2 * (1-Gamma) * ZZj)
push!(gates_ZZ, Gj)
#println("Gate ZZ: Src ", si, " Dst ", sj)
end
# order of gates e^H_ZZ/2 * e^H_X/2 * e^H_X/2 * e^H_ZZ/2 |psi>
gates = vcat(gates_ZZ,gates_X)
append!(gates, reverse(gates))
######### Observables
M_z_vals_mps = []
M_z_vals_tn = []
#MPS
psi_mps = psi_init_mps
Z_val_mps = ITensors.expect(psi_mps, "Z")
push!(M_z_vals_mps,sum(Z_val_mps)/n)
#TNS
psi_tn = psi_init_tn
psi_tn, mts = vidal_to_symmetric_gauge(psi_tn, bond_tensors)
psipsi_tn = norm_network(psi_tn)
Z_val_tn = [expect_state_SBP(op("Z", sites[v]), psi_tn, psipsi_tn, mts) for v in vertices(psi_tn)]
push!(M_z_vals_tn,sum(Z_val_tn)/n) # Magnetization
# Comparison Values
psi_arr_vals_mps = []
psi_arr_vals_tn = []
psi_arr_val_mps = mps_to_array(psi_init_mps)
push!(psi_arr_vals_mps,psi_arr_val_mps)
psi_arr_val_tn = tensornetwork_to_array(psi_init_tn)
push!(psi_arr_vals_tn,psi_arr_val_tn)
# Distance To Vidal Gauge
dist_vidal_vals = []
cur_C = vidal_itn_canonicalness(psi_tn, bond_tensors)
push!(dist_vidal_vals,cur_C)
######### Time Evolutions
apply_kwargs = (; maxdim=32, cutoff=1e-14) # (;cutoff=1e-5) #(; maxdim=8) (; maxdim=8, cutoff=1e-12)
for t in ttimes[2:end]
for gate in gates
psi_mps = apply(gate, psi_mps; apply_kwargs...)
psi_tn, bond_tensors = apply(gate, psi_tn, bond_tensors; normalize=true, apply_kwargs...) #; apply_kwargs...)
end
cur_C = vidal_itn_canonicalness(psi_tn, bond_tensors)
push!(dist_vidal_vals,cur_C)
psi_tn, mts = vidal_to_symmetric_gauge(psi_tn, bond_tensors);
bond_tensors = initialize_bond_tensors(psi_tn);
psi_tn, bond_tensors = vidal_gauge(psi_tn, mts, bond_tensors; eigen_message_tensor_cutoff=1e-14, regularization=1e-14)
# MPS Magnetization
Z_val_mps = ITensors.expect(psi_mps, "Z")
push!(M_z_vals_mps,sum(Z_val_mps)/n)
# Normalization
ITensors.normalize!(psi_mps)
# TN Magnetization
psi_tn, mts = vidal_to_symmetric_gauge(psi_tn, bond_tensors)
psipsi_tn = norm_network(psi_tn)
Z_val_tn = [expect_state_SBP(op("Z", sites[v]), psi_tn, psipsi_tn, mts) for v in vertices(psi_tn)]
push!(M_z_vals_tn,sum(Z_val_tn)/n)
# Comparison Values
psi_arr_val_mps = mps_to_array(psi_mps)
push!(psi_arr_vals_mps,psi_arr_val_mps)
psi_arr_val_tn = tensornetwork_to_array(psi_tn)
push!(psi_arr_vals_tn,psi_arr_val_tn)
end
println("Comparison Evolution Done")
overlap_tn_mps = real([sum(conj(psi_arr_vals_mps[idx]) .* psi_arr_vals_tn[idx]) for idx in 1:length(psi_arr_vals_mps)])
p1 = Plots.plot(ttimes, real(M_z_vals_mps), label="⟨Mz⟩ MPS", xlabel="t")
Plots.scatter!(p1, ttimes, real(M_z_vals_tn), label="⟨Mz⟩ TN")
Plots.plot!(p1, ttimes, overlap_tn_mps, label="Overlap TN & MPS")
p2 = Plots.scatter(ttimes,dist_vidal_vals,label="Distance to Vidal Gauge for TN")
plot_combined = Plots.plot(p1, p2, layout=(1, 2), legend=true, size=(1000,400))