Iterator-based sweep solvers

docs/src/changelog.md

@@ -24,6 +24,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - `ClusterExpansion` time-evolution algorithm for `make_time_mpo`, constructing evolution MPOs
   for finite and infinite nearest-neighbor Hamiltonians through a nonperturbative cluster expansion
   that matches exact exponentials on clusters of up to `N` sites. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252))
+- In-place `timestep!` for an `InfiniteMPS` with `TDVP`, which also enables `time_evolve!` on
+  infinite states. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
 - `changebonds(::InfiniteMPOHamiltonian, ::SvdCut)` compresses infinite Hamiltonians with the
   almost-Schmidt decomposition of [Parker et al.](https://arxiv.org/abs/1909.06341), and
   `left_canonicalize!`/`right_canonicalize!` bring an `InfiniteMPOHamiltonian` into canonical form. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539))
@@ -33,6 +35,17 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - `left_canonicalize!`/`right_canonicalize!(::FiniteMPOHamiltonian)` now sweep over the whole chain
   instead of a single site, and return the gauge transformations alongside the operator, matching
   the infinite versions. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539))
+- Internal: the finite `DMRG`/`DMRG2` ground-state search and approximation, `DynamicalDMRG`,
+  multiline `IDMRG`/`IDMRG2` approximation and `leading_boundary`, and `time_evolve` now run as
+  iterators whose iterations are a forward and a backward half-sweep (or a single time step),
+  sharing one structure across solvers. Finite `TDVP`, `TDVP2` and `BUG` steps share a single sweep
+  over per-site local updates, and infinite `TDVP` steps are split into local update, gauge and
+  environment steps on the same step state. Results and logging are unchanged. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+- `time_evolve` with finite `TDVP`, `TDVP2` or `BUG`, or infinite `TDVP`, copies (or promotes to
+  complex) the initial state once and then takes every step in place, reusing one step state and
+  its allocator across steps, rather than copying the state at every step. A `finalize` callback
+  that keeps the state of an earlier step therefore sees it updated by later steps, as with
+  `time_evolve!`. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
 
 ### Deprecated
 
@@ -44,6 +57,13 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - `show` methods for finite MPS/MPOs no longer check unit space condition through `oneunit`, which is not defined for multifusion spaces. ([#538](https://github.com/QuantumKitHub/MPSKit.jl/pull/538))
 - The canonical forms of `MPOHamiltonian`s, and therefore `SvdCut` compression, are now correct when
   interactions contain identity components or redundant channels. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539))
+- `InfiniteMPS(C, AR)` no longer throws a `MethodError`, which also fixes infinite `TDVP` with
+  `leftorthflag = false`. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+- Real-time infinite `TDVP` from a real `InfiniteMPS` with a real Hamiltonian no longer fails with
+  an `InexactError`: the environments are now rebuilt for the complex state, as for finite MPS.
+  ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+- `TDVP`, `TDVP2` and `BUG` promote a real state to complex for imaginary-time evolution with a
+  complex time step, rather than failing. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
 
 ### Performance
 

src/algorithms/approximate/fvomps.jl

@@ -1,74 +1,82 @@
-function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environments(ψ, _environment_args(Oϕ)...))
+# Internal state of the finite DMRG/DMRG2 approximation, where `ϵ` is the largest relative
+# local change of the last sweep
+struct ApproximateState{S, O, E, A}
+    mps::S
+    operator::O
+    envs::E
+    iter::Int
+    ϵ::Float64
+    allocator::A
+end
+
+function approximate!(
+        ψ::AbstractFiniteMPS, Oϕ, alg::Union{DMRG, DMRG2},
+        envs = environments(ψ, _environment_args(Oϕ)...)
+    )
     allocator = default_allocator(ψ, SerialScheduler())
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
-    log = IterLog("DMRG2")
+    log = IterLog(alg)
+    it = IterativeSolver(alg, ApproximateState(ψ, Oϕ, envs, 0, 2 * alg.tol, allocator))
 
     with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            ϵ = 0.0
-            for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1]
-                AC2′ = AC2_projection(pos, ψ, Oϕ, envs; alg.backend, allocator)
-                al, c, ar, = svd_trunc!(AC2′, inner_alg_gauge(alg))
-
-                AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
-                ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))
-
-                ψ.AC[pos] = (al, complex(c))
-                ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
-            end
-
-            # finalize
-            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ), typeof(envs)}
-
+        @log_initialization loginit!(log, it.ϵ)
+        for (_, _, ϵ) in Iterators.take(it, alg.maxiter)
             if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
+                @log_convergence logfinish!(log, it.iter, ϵ)
                 break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
+            elseif it.iter == alg.maxiter
+                @log_nonconvergence logcancel!(log, it.iter, ϵ)
             else
-                @log_iteration logiter!(log, iter, ϵ)
+                @log_iteration logiter!(log, it.iter, ϵ)
             end
         end
     end
 
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, localchange = ϵ, numiter = iter)
+    state = it.state
+    info = AlgorithmInfo(; converged = state.ϵ <= alg.tol, localchange = state.ϵ, numiter = state.iter)
+    return state.mps, state.envs, info
 end
 
-function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments(ψ, _environment_args(Oϕ)...))
-    allocator = default_allocator(ψ, SerialScheduler())
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
-    log = IterLog("DMRG")
-
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            ϵ = 0.0
-            for pos in [1:(length(ψ) - 1); length(ψ):-1:2]
-                AC′ = AC_projection(pos, ψ, Oϕ, envs; alg.backend, allocator)
-                AC = ψ.AC[pos]
-                ϵ = max(ϵ, norm(AC′ - AC) / norm(AC′))
+function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::ApproximateState)
+    iter = state.iter + 1
+    state = ApproximateState(
+        state.mps, state.operator, state.envs, state.iter, 0.0, state.allocator
+    )
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    ψ, envs = it.finalize(
+        iter, state.mps, state.operator, state.envs
+    )::Tuple{typeof(state.mps), typeof(state.envs)}
+    it.state = ApproximateState(ψ, state.operator, envs, iter, state.ϵ, state.allocator)
+    return (ψ, envs, state.ϵ), it.state
+end
 
-                ψ.AC[pos] = AC′
-            end
+function sweep!(it::IterativeSolver{<:DMRG2}, state::ApproximateState, direction, iter)
+    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ, Oϕ, envs, ϵ = state.mps, state.operator, state.envs, state.ϵ
+    for pos in sites
+        AC2′ = AC2_projection(pos, ψ, Oϕ, envs; it.alg.backend, state.allocator)
+        al, c, ar, = svd_trunc!(AC2′, inner_alg_gauge(it.alg))
 
-            # finalize
-            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ), typeof(envs)}
+        AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
+        ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))
 
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
-        end
+        ψ.AC[pos] = (al, complex(c))
+        ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
     end
+    return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator)
+end
 
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, localchange = ϵ, numiter = iter)
+function sweep!(it::IterativeSolver{<:DMRG}, state::ApproximateState, direction, iter)
+    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ, Oϕ, envs, ϵ = state.mps, state.operator, state.envs, state.ϵ
+    for pos in sites
+        AC′ = AC_projection(pos, ψ, Oϕ, envs; it.alg.backend, state.allocator)
+        AC = ψ.AC[pos]
+        ϵ = max(ϵ, norm(AC′ - AC) / norm(AC′))
+
+        ψ.AC[pos] = AC′
+    end
+    return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator)
 end