VUMPS with Infinite MPO (#50)

Author: LHerviou

examples/vumps/development/vumps_mpo.jl

@@ -2,62 +2,43 @@ using ITensors
 using ITensorInfiniteMPS
 using Random
 
-# Helper function to make an MPO
-import ITensors: op
-op(::OpName"Zero", ::SiteType"S=1/2", s::Index) = ITensor(s', dag(s))
-
-function tfi_mpo(s, l, r; J=1.0, h)
-  dₗ = 3 # The link dimension of the TFI
-  Hmat = fill(op("Zero", s), dₗ, dₗ)
-  Hmat[1, 1] = op("Id", s)
-  Hmat[2, 2] = op("Id", s)
-  Hmat[3, 3] = op("Id", s)
-  Hmat[2, 1] = -J * op("X", s)
-  Hmat[3, 2] = -J * op("X", s)
-  Hmat[3, 1] = -h * op("Z", s)
-  H = ITensor()
-  for i in 1:dₗ, j in 1:dₗ
-    H += Hmat[i, j] * setelt(l => i) * setelt(r => j)
-  end
-  return H
+# Parameters
+cutoff = 1e-8
+maxdim = 100
+tol = 1e-8
+maxiter = 20
+outer_iters = 3
+
+N = 2
+model = Model"ising"()
+
+function space_shifted(::Model"ising", q̃sz)
+  return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
 end
 
-Random.seed!(1234)
-
-N = 1
-s = siteinds("S=1/2", N)
-χ = 6
-@assert iseven(χ)
-if any(hasqns, s)
-  space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2)
-else
-  space = χ
-end
-#ψ = InfiniteMPS(ComplexF64, s; space = space)
-ψ = InfiniteMPS(s; space=space)
-randn!.(ψ)
-
-# Use a finite MPO to create the infinite MPO
-H = InfiniteMPO(N)
-
-# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ
-J = 1.0
-h = 2.0
-
-# For the finite unit cell
-s¹ = addtags.(s, "c=1")
-l¹ = [Index(3, "Link,l=$n,c=1") for n in 1:N]
-for n in 1:N
-  # TODO: use a CelledVector here to make the code logic simpler
-  l¹ₙ = n == 1 ? replacetags(dag(l¹[N]), "c=1" => "c=0") : l¹[n - 1]
-  r¹ₙ = l¹[n]
-  H[n] = tfi_mpo(s¹[n], l¹ₙ, r¹ₙ; J=J, h=h)
+space_ = fill(space_shifted(model, 1), N)
+s = infsiteinds("S=1/2", N; space=space_)
+initstate(n) = "↑"
+ψ = InfMPS(s, initstate)
+
+model_params = (J=-1.0, h=0.9)
+vumps_kwargs = (multisite_update_alg="sequential", tol=tol, maxiter=maxiter)
+subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
+
+H = InfiniteITensorSum(model, s; model_params...)
+# Alternate steps of running VUMPS and increasing the bond dimension
+ψ1 = vumps(H, ψ; vumps_kwargs...)
+for _ in 1:outer_iters
+  ψ1 = subspace_expansion(ψ1, H; subspace_expansion_kwargs...)
+  ψ1 = vumps(H, ψ1; vumps_kwargs...)
 end
 
-ψf = vumps(H, ψ; nsweeps=10)
+Hmpo = InfiniteMPOMatrix(model, s; model_params...)
+# Alternate steps of running VUMPS and increasing the bond dimension
+ψ2 = vumps(Hmpo, ψ; vumps_kwargs...)
+for _ in 1:outer_iters
+  ψ2 = subspace_expansion(ψ2, Hmpo; subspace_expansion_kwargs...)
+  ψ2 = vumps(Hmpo, ψ2; vumps_kwargs...)
+end
 
-s⃗ = [uniqueind(ψ[n], ψ[n - 1], ψ[n + 1]) for n in 1:2]
-h = -J * op("X", s⃗, 1) * op("X", s⃗, 2) - h * op("Z", s⃗, 1) * op("Id", s⃗, 2)
-#ψ2 = ψf.AL[1] * ψf.C[1] * ψf.AR[2]
-ψ2 = ψf.AL[1] * ψf.AL[2] * ψf.C[2]
-@show ψ2 * h * prime(dag(ψ2), "Site")
+SzSz = prod(op("Sz", s[1]) * op("Sz", s[2]))

src/ITensorInfiniteMPS.jl

@@ -32,6 +32,7 @@ include("abstractinfinitemps.jl")
 include("infinitemps.jl")
 include("infinitempo.jl")
 include("infinitecanonicalmps.jl")
+include("infinitempomatrix.jl")
 include("transfermatrix.jl")
 include("models/models.jl")
 include("models/ising.jl")
@@ -44,6 +45,7 @@ include("infinitemps_approx.jl")
 include("nullspace.jl")
 include("subspace_expansion.jl")
 include("vumps_localham.jl")
+include("vumps_mpo.jl")
 
 export Cell,
   CelledVector,
@@ -52,6 +54,7 @@ export Cell,
   InfMPS,
   InfiniteSum,
   InfiniteMPO,
+  InfiniteMPOMatrix,
   InfiniteSumLocalOps,
   ITensorMap,
   ITensorNetwork,

src/celledvectors.jl

@@ -56,6 +56,7 @@ end
 
 translatecell(T::ITensor, n::Integer) = ITensors.setinds(T, translatecell(inds(T), n))
 translatecell(T::MPO, n::Integer) = translatecell.(T, n)
+translatecell(T::Matrix{ITensor}, n::Integer) = translatecell.(T, n)
 
 struct CelledVector{T} <: AbstractVector{T}
   data::Vector{T}

src/infinitempo.jl

@@ -4,7 +4,7 @@
 #
 
 mutable struct InfiniteMPO <: AbstractInfiniteMPS
-  data::Vector{ITensor}
+  data::CelledVector{ITensor}
   llim::Int #RealInfinity
   rlim::Int #RealInfinity
   reverse::Bool
@@ -16,3 +16,5 @@ end
 ## InfiniteSumLocalOps(N::Int) = InfiniteSumLocalOps(Vector{ITensor}(undef, N))
 ## InfiniteSumLocalOps(data::Vector{ITensor}) = InfiniteSumLocalOps(CelledVector(data))
 ## getindex(l::InfiniteSumLocalOps, n::Integer) = ITensors.data(l)[n]
+
+# TODO? Instead of having a big quasi empty ITensor, store only the non zero blocks

src/infinitempomatrix.jl

@@ -0,0 +1,40 @@
+
+mutable struct InfiniteMPOMatrix <: AbstractInfiniteMPS
+  data::CelledVector{Matrix{ITensor}}
+  llim::Int #RealInfinity
+  rlim::Int #RealInfinity
+  reverse::Bool
+end
+
+# TODO better printing?
+function Base.show(io::IO, M::InfiniteMPOMatrix)
+  print(io, "$(typeof(M))")
+  (length(M) > 0) && print(io, "\n")
+  for i in eachindex(M)
+    if !isassigned(M, i)
+      println(io, "#undef")
+    else
+      A = M[i]
+      println(io, "Matrix tensor of size $(size(A))")
+      for k in 1:size(A)[1], l in 1:size(A)[2]
+        if !isassigned(A, k + (size(A)[1] - 1) * l)
+          println(io, "[($k, $l)] #undef")
+        elseif isempty(A[k, l])
+          println(io, "[($k, $l)] empty")
+        else
+          println(io, "[($k, $l)] $(inds(A[k, l]))")
+        end
+      end
+    end
+  end
+end
+
+function getindex(ψ::InfiniteMPOMatrix, n::Integer)
+  return ψ.data[n]
+end
+
+function InfiniteMPOMatrix(arrMat::Vector{Matrix{ITensor}})
+  return InfiniteMPOMatrix(arrMat, 0, size(arrMat)[1], false)
+end
+
+#nrange(H::InfiniteMPOMatrix) = [size(H[j])[1] - 1 for j in 1:nsites(H)]

src/infinitemps.jl

@@ -88,10 +88,31 @@ struct InfiniteSum{T}
 end
 InfiniteSum{T}(N::Int) where {T} = InfiniteSum{T}(Vector{T}(undef, N))
 InfiniteSum{T}(data::Vector{T}) where {T} = InfiniteSum{T}(CelledVector(data))
-
 # Automatically converts from ITensor to MPO.
 # XXX: check this conversion is correct.
-InfiniteSum{T}(data::Vector{ITensor}) where {T} = InfiniteSum{T}(CelledVector(T.(data)))
+#InfiniteSum{T}(data::Vector{ITensor}) where {T} = InfiniteSum{T}(CelledVector(T.(data)))
+#This crashed for me when T is ITensor
+function InfiniteSum{MPO}(data::Vector{ITensor})
+  N = length(data)
+  temp_inds = [filterinds(data[n]; plev=0) for n in 1:N]
+  return InfiniteSum{MPO}([
+    MPO(
+      data[n],
+      [(temp_inds[n][j], dag(prime(temp_inds[n][j]))) for j in 1:length(temp_inds[n])],
+    ) for n in 1:N
+  ])
+end
+
+function InfiniteSum{MPO}(infsum::InfiniteSum{ITensor})
+  N = nsites(infsum)
+  temp_inds = [filterinds(infsum[n]; plev=0) for n in 1:N]
+  return InfiniteSum{MPO}([
+    MPO(
+      infsum[n],
+      [(temp_inds[n][j], dag(prime(temp_inds[n][j]))) for j in 1:length(temp_inds[n])],
+    ) for n in 1:N
+  ])
+end
 
 function Base.getindex(l::InfiniteSum, n1n2::Tuple{Int,Int})
   n1, n2 = n1n2

src/models/models.jl

@@ -28,6 +28,12 @@ function InfiniteSum{MPO}(model::Model, s::CelledVector; kwargs...)
   return InfiniteSum{MPO}(mpos)
 end
 
+function InfiniteSum{ITensor}(model::Model, s::CelledVector; kwargs...)
+  N = length(s)
+  itensors = [ITensor(model, s, n; kwargs...) for n in 1:N] #slightly improved version. Note: the current implementation does not really allow for staggered potentials for example
+  return InfiniteSum{ITensor}(itensors)
+end
+
 # MPO building version
 function ITensors.MPO(model::Model, s::CelledVector, n::Int64; kwargs...)
   n1, n2 = 1, 2
@@ -39,3 +45,33 @@ end
 function ITensors.ITensor(model::Model, s::CelledVector, n::Int64; kwargs...)
   return prod(MPO(model, s, n; kwargs...)) #modification to allow for more than two sites per term in the Hamiltonians
 end
+
+# Helper function to make an MPO
+import ITensors: op
+op(::OpName"Zero", ::SiteType, s::Index) = ITensor(s', dag(s))
+
+function InfiniteMPOMatrix(model::Model, s::CelledVector; kwargs...)
+  N = length(s)
+  temp_H = InfiniteSum{MPO}(model, s; kwargs...)
+  range_H = nrange(temp_H)[1]
+  ls = CelledVector([Index(1, "Link,c=1,n=$n") for n in 1:N])
+  mpos = [Matrix{ITensor}(undef, 1, 1) for i in 1:N]
+  for j in 1:N
+    Hmat = fill(op("Zero", s[j]), range_H + 1, range_H + 1)
+    identity = op("Id", s[j])
+    Hmat[1, 1] = identity
+    Hmat[end, end] = identity
+    for n in 0:(range_H - 1)
+      idx = findfirst(x -> x == j, findsites(temp_H[j - n]; ncell=N))
+      if isnothing(idx)
+        Hmat[range_H + 1 - n, range_H - n] = identity
+      else
+        Hmat[range_H + 1 - n, range_H - n] = temp_H[j - n][idx]#replacetags(linkinds, temp_H[j - n][idx], "Link, l=$n", tags(ls[j-1]))
+      end
+    end
+    mpos[j] = Hmat
+    #mpos[j] += dense(Hmat) * setelt(ls[j-1] => total_dim) * setelt(ls[j] => total_dim)
+  end
+  #return mpos
+  return InfiniteMPOMatrix(mpos)
+end

src/orthogonalize.jl

@@ -145,3 +145,15 @@ function mixed_canonical(
 end
 
 ITensors.orthogonalize(ψ::InfiniteMPS, ::Colon; kwargs...) = mixed_canonical(ψ; kwargs...)
+
+#TODO put these functions somewhere else
+function ortho_overlap(AC, C)
+  AL, _ = polar(AC * dag(C), uniqueinds(AC, C))
+  return noprime(AL)
+end
+
+function ortho_polar(AC, C)
+  UAC, _ = polar(AC, uniqueinds(AC, C))
+  UC, _ = polar(C, commoninds(C, AC))
+  return noprime(UAC) * noprime(dag(UC))
+end