feat: BUG integrator

docs/src/assets/mpskit.bib

@@ -29,6 +29,21 @@ @article{capponi2025
   abstract = {We investigate the nature of the quantum phase transition in modulated Heisenberg spin chains. In the odd- case, the transition separates a trivial nondegenerate phase to a doubly degenerate gapped chiral symmetry-protected topological (SPT) phase which breaks spontaneously the inversion symmetry. The transition is not an Ising transition associated to the breaking of the inversion symmetry, but is governed by the delocalization of the edge states of the SPT phase. In this respect, a modulated Heisenberg spin chain provides a simple example in one dimension of a non-Landau phase transition which is described by the conformal field theory. We show that the chiral SPT phase exhibits fractionalized spinon excitations, which can be confined by slightly changing the model parameters.}
 }
 
+@article{ceruti2022,
+  title = {An Unconventional Robust Integrator for Dynamical Low-Rank Approximation},
+  author = {Ceruti, Gianluca and Lubich, Christian},
+  year = {2022},
+  month = mar,
+  journal = {BIT Numerical Mathematics},
+  volume = {62},
+  number = {1},
+  pages = {23--44},
+  publisher = {Springer},
+  doi = {10.1007/s10543-021-00873-0},
+  url = {https://doi.org/10.1007/s10543-021-00873-0},
+  abstract = {We propose and analyze a numerical integrator for computing the low-rank approximation to solutions of matrix differential equations. The proposed method is based on a variant of the projector-splitting integrator, but here the sub-steps are chosen such that the numerical integrator is robust to the presence of small singular values in the solution.}
+}
+
 @article{chepiga2017,
   title = {Excitation Spectrum and Density Matrix Renormalization Group Iterations},
   author = {Chepiga, Natalia and Mila, Fr{\'e}d{\'e}ric},

docs/src/changelog.md

@@ -21,8 +21,19 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 ### Added
 
+- `BUG` time-evolution algorithm: a symmetric second-order Basis-Update & Galerkin integrator for
+  finite MPS. Unlike `TDVP` it has no backward-in-time substep (stable for imaginary-time evolution),
+  and passing a truncating `trscheme` enables rank-adaptivity (the bond dimension grows and shrinks
+  automatically to track entanglement).
+
 ### Changed
 
+- Renormalization during time evolution is now controlled by an explicit `normalize` keyword on
+  `timestep`/`time_evolve` (default `false`), decoupled from `imaginary_evolution`. By default the
+  norm is preserved, so it retains useful information (the accumulated truncation error in real time,
+  or the decaying weight in imaginary time). Previously imaginary-time evolution always renormalized
+  every step; **to recover that behavior, pass `normalize = true`** (e.g. for ground-state or
+  thermal-state search via imaginary-time evolution).
 - `environments` now follows a single positional contract for every state and operator kind:
   `environments(below, operator, above, alg)`, where `alg` is the environment algorithm
   (slot 4). The operator form requires an explicit `above`. Auxiliary inputs are keyword-only:

docs/src/man/algorithms.md

@@ -104,6 +104,7 @@ This procedure is commonly referred to as the [`TDVP`](@ref) algorithm, which ag
 ```@docs; canonical=false
 TDVP
 TDVP2
+BUG
 ```
 
 ### Time evolution MPO

src/MPSKit.jl

@@ -35,7 +35,7 @@ export VUMPS, VOMPS, DMRG, DMRG2, IDMRG, IDMRG2, GradientGrassmann
 export excitations
 export FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2
 export time_evolve, timestep, timestep!, make_time_mpo
-export TDVP, TDVP2, WI, WII, TaylorCluster
+export TDVP, TDVP2, BUG, WI, WII, TaylorCluster
 export changebonds, changebonds!
 export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand
 export propagator
@@ -161,6 +161,7 @@ include("algorithms/changebonds/randexpand.jl")
 include("algorithms/changebonds/sketchedexpand.jl")
 
 include("algorithms/timestep/tdvp.jl")
+include("algorithms/timestep/bug.jl")
 include("algorithms/timestep/taylorcluster.jl")
 include("algorithms/timestep/wii.jl")
 include("algorithms/timestep/integrators.jl")

src/algorithms/timestep/bug.jl

@@ -0,0 +1,190 @@
+"""
+$(TYPEDEF)
+
+Single-site, symmetric second-order time-evolution algorithm for **finite** MPS, based on the
+Basis-Update & Galerkin (BUG) integrator, an unconventional robust integrator for dynamical low-rank
+approximation.
+
+Unlike [`TDVP`](@ref), BUG advances both the basis (K-step) and the core (Galerkin C-step) tensors
+*forward* in time, with no backward-in-time substep. This makes it a natural choice for
+imaginary-time / dissipative evolution, where the backward core step of TDVP can become unstable.
+
+Each half-sweep augments every bond with the new directions discovered by the evolved connecting
+tensor (old basis first, `[U₀ │ K₁]`) and truncates back down to the tolerance of `trscheme` in the
+same orthonormalization step. The truncation is folded into the augment orth (a truncated SVD of the
+stacked basis), so the old directions — whose orthonormal columns carry singular value `≥ 1` — are
+always kept and only newly-appended directions are cut. The bond therefore *grows* to track the
+entanglement under a `truncerror` tolerance; a hard rank cap (`truncrank`) can additionally shrink it.
+
+!!! warning
+    With `trscheme = notrunc()` the augmented basis is kept at full rank every half-sweep, so the
+    bond grows unboundedly (up to the local Hilbert-space dimension) — this is **not** a fixed-rank
+    integrator. To cap the bond dimension at `D`, pass `trscheme = truncrank(D)`: because the old
+    directions carry singular value `≥ 1`, a rank-`D` cut keeps them and drops the newly appended
+    directions, so a bond already at `D` stays there (a fixed-rank-`D` step). A bond currently below
+    `D` can still grow up to `D` as new directions are admitted.
+
+!!! note
+    By default the state is not renormalized, so the norm keeps useful information (the
+    accumulated truncation error in real time, or the decaying weight in imaginary time).
+    Pass `normalize = true` to `timestep`/`time_evolve` to renormalize after every half-sweep
+    instead. This is independent of `imaginary_evolution`.
+
+## Fields
+
+$(TYPEDFIELDS)
+
+## References
+
+* [Ceruti et al. BIT Numer. Math. 62 (2022)](@cite ceruti2022)
+"""
+struct BUG{A, O, G, F} <: Algorithm
+    "algorithm used in the exponential solvers"
+    integrator::A
+
+    "algorithm used to re-orthonormalize the basis after each local update"
+    alg_orth::O
+
+    "factorization used for the in-sweep augment gauge: a truncated SVD (`alg_svd` with `trscheme`)"
+    alg_gauge::G
+
+    "callback function applied after each iteration, of signature `finalize(t, ψ, H, envs) -> ψ, envs`"
+    finalize::F
+end
+function BUG(;
+        integrator = Defaults.alg_expsolve(), alg_orth = Defaults.alg_orth(),
+        trscheme = notrunc(), alg_svd = Defaults.alg_svd(),
+        finalize = Defaults._finalize
+    )
+    alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
+    return BUG(integrator, alg_orth, alg_gauge, finalize)
+end
+
+# left→right BUG update at `site`: evolve the connecting tensor, old-first augment + truncate the
+# basis, and return the `transport` and center `AC_old` advanced to `site + 1` (root = last site).
+function local_update!(
+        site, ::Val{:right}, ψ, H, alg::BUG, envs, t, h, transport, AC_old, ARs;
+        imaginary_evolution, normalize
+    )
+    # 1. Transport the old AC to the new frame and evolve
+    AC₀ = _mul_front(transport, AC_old)
+    AC = integrate(
+        AC_hamiltonian(site, ψ, H, ψ, envs), AC₀,
+        t, h, alg.integrator; imaginary_evolution
+    )
+
+    # If end of chain, finalize
+    if site == length(ψ)
+        normalize && normalize!(AC)
+        ψ.AC[site] = AC
+        return ψ, transport, AC_old
+    end
+
+    # 2. move gauge right
+    oldbasis, C₀ = left_gauge(AC₀, alg.alg_orth)                        # old AL - C in new frame
+    AC_next = _mul_front(C₀, ARs[site + 1])
+
+
+    # 3. reproject onto new basis
+    newbasis, _ = left_gauge(catdomain(oldbasis, AC), alg.alg_gauge)    # augmented basis
+    newbond = newbasis' * AC
+    ψ.AC[site] = (newbasis, newbond)
+    new_transport = newbasis' * oldbasis
+
+    return ψ, new_transport, AC_next
+end
+
+# mirror of the left→right update on the `_transpose_tail` form; `transport` and `AC_old` are returned
+# advanced to `site - 1` (root = first site).
+function local_update!(
+        site, ::Val{:left}, ψ, H, alg::BUG, envs, t, h, transport, AC_old, ALs;
+        imaginary_evolution, normalize
+    )
+    # 1. Transport the old AC to the new frame and evolve
+    AC₀ = AC_old * transport
+    AC = integrate(
+        AC_hamiltonian(site, ψ, H, ψ, envs), AC₀,
+        t, h, alg.integrator; imaginary_evolution
+    )
+
+    # If start of chain, finalize
+    if site == 1
+        normalize && normalize!(AC)
+        ψ.AC[site] = AC
+        return ψ, transport, AC_old
+    end
+
+    # 2. move gauge left
+    C₀, oldbasis = right_gauge(AC₀, alg.alg_orth)              # old C - AR in new frame
+    AC_next = _mul_tail(ALs[site - 1], C₀)
+
+    # 3. reproject onto new basis
+    _, newbasis_tail = right_orth(
+        catcodomain(_transpose_tail(oldbasis), _transpose_tail(AC));
+        alg = MatrixAlgebraKit.RightOrthViaSVD(alg.alg_gauge)   # augmented basis
+    )
+    newbasis = _transpose_front(newbasis_tail)
+    newbond = _transpose_tail(AC) * _transpose_tail(newbasis)'
+    ψ.AC[site] = (newbond, newbasis)
+    new_transport = _transpose_tail(oldbasis) * _transpose_tail(newbasis)'
+
+    return ψ, new_transport, AC_next
+end
+
+function timestep!(
+        ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::BUG,
+        envs::AbstractMPSEnvironments = environments(ψ, H, ψ);
+        imaginary_evolution::Bool = false, normalize::Bool = false
+    )
+    # symmetric 2nd-order: a dt/2 left→right half-sweep composed with its dt/2 mirror
+    L = length(ψ)
+    h = dt / 2
+
+    # left→right half-sweep (root = last site): freeze the bases as `[AC[1], AR[2], …, AR[L]]` and
+    # carry `AC_old`/`transport` forward
+    ARs = ψ.AR[2:end]
+    pushfirst!(ARs, ψ.AC[1])
+    transport = isomorphism(scalartype(ψ), left_virtualspace(ψ, 1) ← left_virtualspace(ψ, 1))
+    AC_old = ARs[1]
+    for site in 1:L
+        ψ, transport, AC_old = local_update!(
+            site, Val(:right), ψ, H, alg, envs, t, h, transport, AC_old, ARs;
+            imaginary_evolution, normalize
+        )
+    end
+
+    # right→left half-sweep (root = first site), the mirror: freeze `[AL[1], …, AL[L-1], AC[L]]` and
+    # carry `AC_old`/`transport` backward (starting from `t + h`)
+    ALs = ψ.AL[1:(L - 1)]
+    push!(ALs, ψ.AC[L])
+    transport = isomorphism(scalartype(ψ), right_virtualspace(ψ, L) ← right_virtualspace(ψ, L))
+    AC_old = ALs[L]
+    for site in L:-1:1
+        ψ, transport, AC_old = local_update!(
+            site, Val(:left), ψ, H, alg, envs, t + h, h, transport, AC_old, ALs;
+            imaginary_evolution, normalize
+        )
+    end
+
+    return ψ, envs
+end
+
+# copying version
+function timestep(
+        ψ::AbstractFiniteMPS, H, time::Number, timestep::Number,
+        alg::BUG, envs::AbstractMPSEnvironments...;
+        imaginary_evolution::Bool = false, normalize::Bool = false, kwargs...
+    )
+    isreal = (scalartype(ψ) <: Real && !imaginary_evolution)
+    ψ′ = isreal ? complex(ψ) : copy(ψ)
+    if length(envs) != 0 && isreal
+        @warn "Currently cannot reuse real environments for complex evolution"
+        envs′ = environments(ψ′, H, ψ′)
+    elseif length(envs) == 1
+        envs′ = only(envs)
+    else
+        @assert length(envs) == 0 "Invalid signature"
+        envs′ = environments(ψ′, H, ψ′)
+    end
+    return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...)
+end