IDMRG updates (#294)

Improvements and changes to more easily extend IDMRG methods for `approximate`, `leading_boundary` and `find_groundstate`, with arbitrary number of physical legs.

Author: lkdvos

src/MPSKit.jl

@@ -173,6 +173,7 @@ include("algorithms/statmech/leading_boundary.jl")
 include("algorithms/statmech/vumps.jl")
 include("algorithms/statmech/vomps.jl")
 include("algorithms/statmech/gradient_grassmann.jl")
+include("algorithms/statmech/idmrg.jl")
 
 include("algorithms/fidelity_susceptibility.jl")
 

src/algorithms/approximate/idmrg.jl

@@ -21,7 +21,7 @@ function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
             end
 
             # right to left sweep
-            for col in size(ψ, 2):-1:1
+            for col in reverse(1:size(ψ, 2))
                 for row in 1:size(ψ, 1)
                     ψ.AC[row + 1, col] = AC_projection(CartesianIndex(row, col),
                                                        ψ, toapprox, envs)
@@ -71,38 +71,42 @@ function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
             C_current = ψ.C[:, 0]
 
             # sweep from left to right
-            for col in 1:size(ψ, 2)
+            for site in 1:size(ψ, 2)
                 for row in 1:size(ψ, 1)
-                    AC2′ = AC2_projection(CartesianIndex(row, col), ψ, toapprox, envs)
+                    AC2′ = AC2_projection(CartesianIndex(row, site), ψ, toapprox, envs;
+                                          kind=:ACAR)
                     al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme, alg=alg.alg_svd)
                     normalize!(c)
 
-                    ψ.AL[row + 1, col] = al
-                    ψ.C[row + 1, col] = complex(c)
-                    ψ.AR[row + 1, col + 1] = _transpose_front(ar)
-                    ψ.AC[row + 1, col + 1] = _transpose_front(c * ar)
+                    ψ.AL[row + 1, site] = al
+                    ψ.C[row + 1, site] = complex(c)
+                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+                    ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
                 end
-                transfer_leftenv!(envs, ψ, toapprox, col + 1)
-                transfer_rightenv!(envs, ψ, toapprox, col)
+
+                transfer_leftenv!(envs, ψ, toapprox, site + 1)
+                transfer_rightenv!(envs, ψ, toapprox, site)
             end
+
             normalize!(envs, ψ, toapprox)
 
             # sweep from right to left
-            for col in (size(ψ, 2) - 1):-1:0
+            for site in reverse(0:(size(ψ, 2) - 1))
                 for row in 1:size(ψ, 1)
-                    # TODO: also write this as AC2_projection?
-                    AC2 = ϕ.AL[row, col] * _transpose_tail(ϕ.AC[row, col + 1])
-                    AC2′ = AC2_hamiltonian(CartesianIndex(row, col), ψ, O, ϕ, envs) * AC2
+                    AC2′ = AC2_projection(CartesianIndex(row, site), ψ, toapprox, envs;
+                                          kind=:ALAC)
                     al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme, alg=alg.alg_svd)
                     normalize!(c)
 
-                    ψ.AL[row + 1, col] = al
-                    ψ.C[row + 1, col] = complex(c)
-                    ψ.AR[row + 1, col + 1] = _transpose_front(ar)
+                    ψ.AL[row + 1, site] = al
+                    ψ.C[row + 1, site] = complex(c)
+                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
                 end
-                transfer_leftenv!(envs, ψ, toapprox, col + 1)
-                transfer_rightenv!(envs, ψ, toapprox, col)
+
+                transfer_leftenv!(envs, ψ, toapprox, site + 1)
+                transfer_rightenv!(envs, ψ, toapprox, site)
             end
+
             normalize!(envs, ψ, toapprox)
 
             # update error
@@ -127,7 +131,7 @@ function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
 
     # TODO: immediately compute in-place
     alg_gauge = updatetol(alg.alg_gauge, iter, ϵ)
-    ψ′ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
+    ψ′ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
     copy!(ψ, ψ′) # ensure output destination is unchanged
 
     recalculate!(envs, ψ, toapprox)

src/algorithms/derivatives/derivatives.jl

@@ -164,20 +164,21 @@ for kind in (:C, :AC, :AC2)
     projection = Symbol(kind, "_projection")
     hamiltonian = Symbol(kind, "_hamiltonian")
 
-    @eval function $projection(site, below, above::Tuple, envs)
-        return $projection(site, below, above..., envs)
+    @eval function $projection(site, below, above::Tuple, envs; kwargs...)
+        return $projection(site, below, above..., envs; kwargs...)
     end
-    @eval function $projection(site, below, above::AbstractMPS, envs)
-        return $projection(site, below, nothing, above, envs)
+    @eval function $projection(site, below, above::AbstractMPS, envs; kwargs...)
+        return $projection(site, below, nothing, above, envs; kwargs...)
     end
     @eval function $projection(site::CartesianIndex{2}, below::MultilineMPS, operator,
-                               above::MultilineMPS, envs)
+                               above::MultilineMPS, envs; kwargs...)
         row, col = Tuple(site)
-        return $projection(col, below[row + 1], operator[row], above[row], envs[row])
+        return $projection(col, below[row + 1], operator[row], above[row], envs[row];
+                           kwargs...)
     end
-    @eval function $projection(site, below, above::LazySum, envs)
+    @eval function $projection(site, below, above::LazySum, envs; kwargs...)
         return sum(zip(above.ops, envs.envs)) do x
-            return $projection(site, below, x...)
+            return $projection(site, below, x...; kwargs...)
         end
     end
 end
@@ -187,9 +188,8 @@ end
 function AC_projection(site, below, operator, above, envs)
     return AC_hamiltonian(site, below, operator, above, envs) * above.AC[site]
 end
-function AC2_projection(site::Int, below, operator, above, envs)
-    AC2 = above.AC[site] * _transpose_tail(above.AR[site + 1])
-    return AC2_hamiltonian(site, below, operator, above, envs) * AC2
+function AC2_projection(site::Int, below, operator, above, envs; kwargs...)
+    return AC2_hamiltonian(site, below, operator, above, envs) * AC2(above, site; kwargs...)
 end
 
 # Multiline

src/algorithms/groundstate/idmrg.jl

@@ -40,7 +40,12 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG, envs=environments(ost
             for pos in 1:length(ψ)
                 h = AC_hamiltonian(pos, ψ, H, ψ, envs)
                 _, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve)
-                ψ.AL[pos], ψ.C[pos] = leftorth!(ψ.AC[pos])
+                if pos == length(ψ)
+                    # AC needed in next sweep
+                    ψ.AL[pos], ψ.C[pos] = leftorth(ψ.AC[pos])
+                else
+                    ψ.AL[pos], ψ.C[pos] = leftorth!(ψ.AC[pos])
+                end
                 transfer_leftenv!(envs, ψ, H, ψ, pos + 1)
             end
 
@@ -124,7 +129,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os
 
             # sweep from left to right
             for pos in 1:(length(ψ) - 1)
-                ac2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
+                ac2 = AC2(ψ, pos; kind=:ACAR)
                 h_ac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs)
                 _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve)
 
@@ -141,20 +146,22 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os
             end
 
             # update the edge
-            @plansor ac2[-1 -2; -3 -4] := ψ.AC[end][-1 -2; 1] * inv(ψ.C[0])[1; 2] *
-                                          ψ.AL[1][2 -4; 3] * ψ.C[1][3; -3]
+            ψ.AL[end] = ψ.AC[end] / ψ.C[end]
+            ψ.AC[1] = _mul_tail(ψ.AL[1], ψ.C[1])
+            ac2 = AC2(ψ, 0; kind=:ALAC)
             h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs)
             _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve)
 
             al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=alg.alg_svd)
             normalize!(c)
 
-            ψ.AC[end] = al * c
             ψ.AL[end] = al
             ψ.C[end] = complex(c)
             ψ.AR[1] = _transpose_front(ar)
+
+            ψ.AC[end] = _mul_tail(al, c)
             ψ.AC[1] = _transpose_front(c * ar)
-            ψ.AL[1] = ψ.AC[1] * inv(ψ.C[1])
+            ψ.AL[1] = ψ.AC[1] / ψ.C[1]
 
             C_current = complex(c)
 
@@ -164,15 +171,15 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os
 
             # sweep from right to left
             for pos in (length(ψ) - 1):-1:1
-                ac2 = ψ.AL[pos] * _transpose_tail(ψ.AC[pos + 1])
+                ac2 = AC2(ψ, pos; kind=:ALAC)
                 h_ac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs)
                 _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve)
 
                 al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=alg.alg_svd)
                 normalize!(c)
 
                 ψ.AL[pos] = al
-                ψ.AC[pos] = al * c
+                ψ.AC[pos] = _mul_tail(al, c)
                 ψ.C[pos] = complex(c)
                 ψ.AR[pos + 1] = _transpose_front(ar)
                 ψ.AC[pos + 1] = _transpose_front(c * ar)
@@ -182,17 +189,19 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os
             end
 
             # update the edge
-            @plansor ac2[-1 -2; -3 -4] := ψ.C[end - 1][-1; 1] * ψ.AR[end][1 -2; 2] *
-                                          inv(ψ.C[end])[2; 3] * ψ.AC[1][3 -4; -3]
+            ψ.AC[end] = _mul_front(ψ.C[end - 1], ψ.AR[end])
+            ψ.AR[1] = _transpose_front(ψ.C[end] \ _transpose_tail(ψ.AC[1]))
+            ac2 = AC2(ψ, 0; kind=:ACAR)
             h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs)
             _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve)
             al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=alg.alg_svd)
             normalize!(c)
 
-            ψ.AR[end] = _transpose_front(inv(ψ.C[end - 1]) * _transpose_tail(al * c))
             ψ.AL[end] = al
             ψ.C[end] = complex(c)
             ψ.AR[1] = _transpose_front(ar)
+
+            ψ.AR[end] = _transpose_front(ψ.C[end - 1] \ _transpose_tail(al * c))
             ψ.AC[1] = _transpose_front(c * ar)
 
             transfer_leftenv!(envs, ψ, H, ψ, 1)