From dc71343be9da016469535d07287cfb77abaf6204 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 14 Jul 2026 17:50:58 -0400 Subject: [PATCH 1/6] time_evolve: add `normalize` kwarg, decoupled from `imaginary_evolution` Renormalization after each step is now controlled by an explicit `normalize` kwarg on `timestep`/`time_evolve` (default `false`) instead of being tied to `imaginary_evolution`. By default the norm is preserved, so it retains useful information (accumulated truncation error in real time, decaying weight in imaginary time); pass `normalize = true` to renormalize, e.g. for ground-state or thermal-state search. Threaded through the driver and both TDVP/TDVP2 integrators; test norm/imaginary-time expectations updated accordingly. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/algorithms/timestep/tdvp.jl | 44 +++++++++++++++----------- src/algorithms/timestep/time_evolve.jl | 14 ++++++-- test/algorithms/timestep.jl | 22 +++++++++---- 3 files changed, 54 insertions(+), 26 deletions(-) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 993782e54..afaf3ba4e 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -11,10 +11,11 @@ the enlarged bond back down (selecting the truncated-SVD gauge). The expansion i state-preserving, as required for a consistent time evolution. !!! note - Real-time evolution preserves the norm: neither the bond expansion nor the truncation - renormalizes, so the state norm reflects the accumulated truncation error. Imaginary-time - evolution instead renormalizes at every step, like a ground-state search. CBE is only - available for finite MPS. + By default the norm is preserved: neither the bond expansion nor the truncation + renormalizes, so the state 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 at every step instead, like a ground-state + search. This is independent of `imaginary_evolution`. CBE is only available for finite MPS. ## Fields @@ -61,11 +62,14 @@ end function timestep( ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs::AbstractMPSEnvironments = environments(ψ, H, ψ); - leftorthflag = true, imaginary_evolution::Bool = false + leftorthflag = true, imaginary_evolution::Bool = false, normalize::Bool = false ) + # `normalize` is accepted for signature uniformity with the finite integrators, but an + # `InfiniteMPS` is always normalized to norm-1-per-site by the gauge/reconstruction below + # (a structural gauge requirement, not information erasure), so the flag has no effect here. # convert state to complex if necessary if scalartype(ψ) <: Real && (!imaginary_evolution || !isreal(dt)) - return timestep(complex(ψ), H, t, dt, alg, envs; leftorthflag, imaginary_evolution) + return timestep(complex(ψ), H, t, dt, alg, envs; leftorthflag, imaginary_evolution, normalize) end temp_ACs = similar(ψ.AC) @@ -120,23 +124,23 @@ end function timestep!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs::AbstractMPSEnvironments = environments(ψ, H, ψ); - imaginary_evolution::Bool = false + imaginary_evolution::Bool = false, normalize::Bool = false ) # sweep left to right for i in 1:(length(ψ) - 1) # 1. optionally expand the bond ahead of the local update (CBE) isnothing(alg.alg_expand) || - changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) + changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize) # 2. evolve the (possibly expanded) center tensor forward Hac = AC_hamiltonian(i, ψ, H, ψ, envs) AC = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution) # 3. gauge: split AC -> AL[i], C[i] (QR center-move, or truncated SVD cutting the - # enlarged bond back down) and move the center to i+1. Real-time evolution preserves - # the norm; imaginary-time evolution renormalizes. - left_gauge!(ψ, i, AC, alg.alg_gauge; normalize = imaginary_evolution) + # enlarged bond back down) and move the center to i+1. By default the norm is + # preserved; `normalize` renormalizes. + left_gauge!(ψ, i, AC, alg.alg_gauge; normalize) # 4. evolve the bond tensor backward Hc = C_hamiltonian(i, ψ, H, ψ, envs) @@ -154,7 +158,7 @@ function timestep!( for i in length(ψ):-1:2 # 1. optionally expand the bond ahead of the local update (CBE) isnothing(alg.alg_expand) || - changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize = imaginary_evolution) + changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize) # 2. evolve the (possibly expanded) center tensor forward Hac = AC_hamiltonian(i, ψ, H, ψ, envs) @@ -163,9 +167,9 @@ function timestep!( imaginary_evolution ) - # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 (real-time preserves the - # norm; imaginary-time renormalizes) - right_gauge!(ψ, i, AC, alg.alg_gauge; normalize = imaginary_evolution) + # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 (norm preserved by + # default; `normalize` renormalizes) + right_gauge!(ψ, i, AC, alg.alg_gauge; normalize) # 4. evolve the bond tensor backward Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs) @@ -221,7 +225,7 @@ end function timestep!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, envs::AbstractMPSEnvironments = environments(ψ, H, ψ); - imaginary_evolution::Bool = false + imaginary_evolution::Bool = false, normalize::Bool = false ) # sweep left to right @@ -231,6 +235,9 @@ function timestep!( ac2′ = integrate(Hac2, ac2, t, dt / 2, alg.integrator; imaginary_evolution) nal, nc, nar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + # `nc` is the norm-carrying bond tensor (`nal`/`nar` are isometries), so normalizing it + # normalizes the whole state, mirroring single-site TDVP's per-gauge renormalization + normalize && normalize!(nc) ψ.AC[i] = (nal, complex(nc)) ψ.AC[i + 1] = (complex(nc), _transpose_front(nar)) @@ -250,6 +257,7 @@ function timestep!( ac2′ = integrate(Hac2, ac2, t + dt / 2, dt / 2, alg.integrator; imaginary_evolution) nal, nc, nar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + normalize && normalize!(nc) ψ.AC[i - 1] = (nal, complex(nc)) ψ.AC[i] = (complex(nc), _transpose_front(nar)) @@ -269,7 +277,7 @@ end function timestep( ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, alg::Union{TDVP, TDVP2}, envs::AbstractMPSEnvironments...; - imaginary_evolution::Bool = false, kwargs... + imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... ) isreal = (scalartype(ψ) <: Real && !imaginary_evolution) ψ′ = isreal ? complex(ψ) : copy(ψ) @@ -282,5 +290,5 @@ function timestep( @assert length(envs) == 0 "Invalid signature" envs′ = environments(ψ′, H, ψ′) end - return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, kwargs...) + return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) end diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 51288abd7..a321520c0 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -20,6 +20,11 @@ through each of the time points obtained by iterating t_span. instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this function to compute the ground state of a Hamiltonian, or to compute finite-temperature properties of a system. +- `normalize::Bool=false`: if true, the state is renormalized after every step. This is + independent of `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). Pass `true` to renormalize each step, e.g. when computing a ground + state or thermal state via imaginary-time evolution. """ function time_evolve end, function time_evolve! end @@ -27,7 +32,7 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time @eval function $time_evolve( ψ, H, t_span::AbstractVector{<:Number}, alg, envs = environments(ψ, H, ψ); - verbosity::Int = 0, imaginary_evolution::Bool = false + verbosity::Int = 0, imaginary_evolution::Bool = false, normalize::Bool = false ) log = IterLog("TDVP") LoggingExtras.withlevel(; verbosity) do @@ -36,7 +41,7 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time t = t_span[iter] dt = t_span[iter + 1] - t - ψ, envs = $timestep(ψ, H, t, dt, alg, envs; imaginary_evolution) + ψ, envs = $timestep(ψ, H, t, dt, alg, envs; imaginary_evolution, normalize) ψ, envs = alg.finalize(t, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)} @infov 3 logiter!(log, iter, 0, t) @@ -69,6 +74,11 @@ solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``. instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this function to compute the ground state of a Hamiltonian, or to compute finite-temperature properties of a system. +- `normalize::Bool=false`: if true, the state is renormalized after the step. This is + independent of `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). Pass `true` to renormalize, e.g. when computing a ground state or + thermal state via imaginary-time evolution. """ function timestep end, function timestep! end diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index a846b2c5c..ede3f5b9c 100644 --- a/test/algorithms/timestep.jl +++ b/test/algorithms/timestep.jl @@ -122,8 +122,9 @@ end @test abs(dot(ref, cbe)) > abs(dot(ref, plain)) end - # the bond truncation must preserve the norm for real-time evolution (the norm reflects the - # discarded weight) and only renormalize for imaginary-time evolution + # by default (`normalize = false`) the bond truncation preserves the norm, so it reflects the + # discarded weight; `normalize = true` renormalizes each step. This is independent of + # `imaginary_evolution`. @testset "norm handling" begin Random.seed!(6) ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart)) @@ -132,15 +133,24 @@ end ψrt = ψ₀ for _ in 1:12 - ψrt, = timestep(ψrt, H, 0.0, 0.5, lossy) # real time + ψrt, = timestep(ψrt, H, 0.0, 0.5, lossy) # real time, norm preserved by default end @test norm(ψrt) < 1 - 1.0e-3 # truncation loss is not renormalized away + # imaginary-time, norm preserved by default: the weight is *not* pinned to unit norm + # (imaginary-time evolution rescales the state, so the norm drifts away from 1) ψit = ψ₀ for _ in 1:12 - ψit, = timestep(ψit, H, 0.0, 0.5, lossy; imaginary_evolution = true) # no external normalize! + ψit, = timestep(ψit, H, 0.0, 0.5, lossy; imaginary_evolution = true) end - @test norm(ψit) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + @test abs(norm(ψit) - 1) > 1.0e-3 + + # imaginary-time with `normalize = true`: renormalized to unit norm each step + ψn = ψ₀ + for _ in 1:12 + ψn, = timestep(ψn, H, 0.0, 0.5, lossy; imaginary_evolution = true, normalize = true) + end + @test norm(ψn) ≈ 1 atol = 1.0e-6 end @testset "imaginary-time lowers energy" begin @@ -150,7 +160,7 @@ end E₀ = real(expectation_value(ψ₀, H)) ψ = ψ₀ for _ in 1:8 - ψ, = timestep(ψ, H, 0.0, 0.1, alg; imaginary_evolution = true) # gauge renormalizes + ψ, = timestep(ψ, H, 0.0, 0.1, alg; imaginary_evolution = true, normalize = true) # gauge renormalizes end @test real(expectation_value(ψ, H)) < E₀ @test dim(left_virtualspace(ψ, L ÷ 2)) > Dstart From bbe9fd49feaf28dfa5a5792ff6248d90687f27e2 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 14 Jul 2026 17:51:08 -0400 Subject: [PATCH 2/6] BUG: sequential Basis-Update & Galerkin integrator for finite MPS Add the `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); passing a truncating `trscheme` enables rank-adaptivity, growing and shrinking the bond dimension to track the entanglement. Wires up the export, include, docs (`@docs` entry, changelog, Ceruti & Lubich 2022 reference). Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/assets/mpskit.bib | 15 ++++ docs/src/changelog.md | 5 ++ docs/src/man/algorithms.md | 1 + src/MPSKit.jl | 3 +- src/algorithms/timestep/bug.jl | 159 +++++++++++++++++++++++++++++++++ 5 files changed, 182 insertions(+), 1 deletion(-) create mode 100644 src/algorithms/timestep/bug.jl diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index 296f4194b..1c2e16e7a 100644 --- a/docs/src/assets/mpskit.bib +++ b/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}, diff --git a/docs/src/changelog.md b/docs/src/changelog.md index f7cd0741e..c1bc75339 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -21,6 +21,11 @@ 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 - `environments` now follows a single positional contract for every state and operator kind: diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 3cf71d52a..694e46fbb 100644 --- a/docs/src/man/algorithms.md +++ b/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 diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 2ef8cf900..d91e2ae55 100644 --- a/src/MPSKit.jl +++ b/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") diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl new file mode 100644 index 000000000..c5f5724c3 --- /dev/null +++ b/src/algorithms/timestep/bug.jl @@ -0,0 +1,159 @@ +""" +$(TYPEDEF) + +Single site MPS time-evolution algorithm 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. + +Passing a truncating `trscheme` (anything other than `notrunc()`) switches on **rank-adaptivity**: +each half-sweep augments every bond with the new directions discovered by the evolved connecting +tensor (old basis first, `[U₀ │ K₁]`) and then truncates back down to the tolerance of `trscheme`. +The bond dimension grows and shrinks automatically to track the entanglement; `notrunc()` recovers +the fixed-rank integrator. + +!!! 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, T, S, F} <: Algorithm + "algorithm used in the exponential solvers" + integrator::A + + "tolerance for gauging algorithm" + tolgauge::Float64 + + "maximal amount of iterations for gauging algorithm" + gaugemaxiter::Int + + "algorithm used to re-orthonormalize the basis after each local update" + alg_orth::O + + "truncation scheme used to cut the bond back down for rank-adaptive BUG" + trscheme::T + + "algorithm used for the singular value decomposition" + alg_svd::S + + "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" + finalize::F +end +function BUG(; + integrator = Defaults.alg_expsolve(), tolgauge = Defaults.tolgauge, + gaugemaxiter = Defaults.maxiter, alg_orth = Defaults.alg_orth(), + trscheme = notrunc(), alg_svd = Defaults.alg_svd(), + finalize = Defaults._finalize + ) + return BUG(integrator, tolgauge, gaugemaxiter, alg_orth, trscheme, alg_svd, finalize) +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 + truncates = !(alg.trscheme isa MatrixAlgebraKit.NoTruncation) + svdcut = SvdCut(; trscheme = alg.trscheme, alg_svd = alg.alg_svd) + + # left→right half-sweep (root = last site) + ψ.AC[1] # gauge center to site 1 + ψ_old = copy(ψ) # frozen bases / reprojection inputs + T = isomorphism(scalartype(ψ), left_virtualspace(ψ_old, 1) ← left_virtualspace(ψ_old, 1)) + for i in 1:(L - 1) + Ĉ = _mul_front(T, ψ_old.AC[i]) # reproject old connecting tensor + AC = integrate(AC_hamiltonian(i, ψ, H, ψ, envs), Ĉ, t, h, alg.integrator; imaginary_evolution) + U₀ = _mul_front(T, ψ_old.AL[i]) # old left isometry in the new frame + if truncates + U = _bug_augment_left(U₀, AC, alg.alg_orth) # old-first augment; cut deferred to changebonds! + C = U' * AC + else + U, C = left_gauge(AC, alg.alg_orth) + end + T = U' * U₀ # transport (new ← old) + ψ.AC[i] = (U, C) + end + AC = integrate( + AC_hamiltonian(L, ψ, H, ψ, envs), _mul_front(T, ψ_old.AC[L]), + t, h, alg.integrator; imaginary_evolution + ) + normalize && normalize!(AC) + ψ.AC[L] = AC + truncates && changebonds!(ψ, svdcut; normalize) + + # right→left half-sweep (root = first site), the mirror + ψ.AC[L] # gauge center to site L + ψ_old = copy(ψ) + T = isomorphism(scalartype(ψ), right_virtualspace(ψ_old, L) ← right_virtualspace(ψ_old, L)) + for i in L:-1:2 + Ĉ = ψ_old.AC[i] * T # reproject old connecting tensor + AC = integrate(AC_hamiltonian(i, ψ, H, ψ, envs), Ĉ, t + h, h, alg.integrator; imaginary_evolution) + U₀ = ψ_old.AR[i] * T # old right isometry in the new frame + if truncates + U = _bug_augment_right(U₀, AC, alg.alg_orth) + C = _transpose_tail(AC) * _transpose_tail(U)' + else + C, U = right_gauge(AC, alg.alg_orth) + end + T = _transpose_tail(U₀) * _transpose_tail(U)' # transport (old ← new) + ψ.AC[i] = (C, U) + end + AC = integrate( + AC_hamiltonian(1, ψ, H, ψ, envs), ψ_old.AC[1] * T, + t + h, h, alg.integrator; imaginary_evolution + ) + normalize && normalize!(AC) + ψ.AC[1] = AC + truncates && changebonds!(ψ, svdcut; normalize) + + return ψ, envs +end + +# augment the RIGHT bond (left→right sweep): orthonormalize the stacked `[U₀ │ K₁]` (old isometry +# first) so `U₀` stays the leading per-sector block and only the new directions of `K₁` are appended. +function _bug_augment_left(U₀, K₁, alg_orth = Defaults.alg_orth()) + Û, _ = left_orth(catdomain(U₀, K₁); alg = alg_orth) + return Û +end + +# mirror of `_bug_augment_left` for the right→left sweep, on the `_transpose_tail` form (right-isometry +# with orthonormal rows): an LQ orthonormalizes the stacked rows `[U₀; K₁]`, keeping `U₀` leading. +function _bug_augment_right(U₀, K₁, alg_orth = Defaults.alg_orth()) + stacked = catcodomain(_transpose_tail(U₀), _transpose_tail(K₁)) + _, Û = right_orth(stacked; alg = alg_orth) + return _transpose_front(Û) +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 From 8cd8d5cfcc27dba03936db152c3e5717d000d0c8 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 14 Jul 2026 17:51:09 -0400 Subject: [PATCH 3/6] BUG: tests at TDVP-test depth Add tests for the sequential `BUG` integrator: energy conservation, eigenstate phase, agreement with TDVP, second-order convergence vs dense propagation, imaginary-time cooling, LazySum/time-dependent operators, rank-adaptivity, and symmetric tensors (U(1)/Z2/SU(2), incl. rank-adaptive-symmetric). Plus unit tests for the `_bug_augment_*` basis-augmentation helpers. Co-Authored-By: Claude Opus 4.8 (1M context) --- test/algorithms/bug.jl | 527 +++++++++++++++++++++++++++++++++ test/algorithms/bug_augment.jl | 136 +++++++++ 2 files changed, 663 insertions(+) create mode 100644 test/algorithms/bug.jl create mode 100644 test/algorithms/bug_augment.jl diff --git a/test/algorithms/bug.jl b/test/algorithms/bug.jl new file mode 100644 index 000000000..d635ebd66 --- /dev/null +++ b/test/algorithms/bug.jl @@ -0,0 +1,527 @@ +println(" +----------------------------- +| BUG time-stepping tests | +----------------------------- +") + +using .TestSetup +using Test, TestExtras +using MPSKit +using TensorKit +using TensorKit: ℙ +using LinearAlgebra: dot, norm +using Random + +@testset "BUG time evolution" verbose = true begin + dt = 0.1 + L = 10 + + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + ψ = FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + # 1. energy conservation + eigenstate phase + @testset "energy conservation" begin + ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + end + + # 2. agreement with TDVP over a few real-time steps of a random MPS + @testset "agreement with TDVP" begin + Random.seed!(1234) + ψr = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + δt = 0.01 + ψ_bug, ψ_tdvp = ψr, ψr + for k in 0:4 + ψ_bug, = timestep(ψ_bug, H, k * δt, δt, BUG()) + ψ_tdvp, = timestep(ψ_tdvp, H, k * δt, δt, TDVP()) + end + @test expectation_value(ψ_bug, H) ≈ expectation_value(ψ_tdvp, H) atol = 1.0e-3 + @test abs(dot(ψ_bug, ψ_tdvp)) ≈ 1 atol = 1.0e-3 + end + + # 3. second-order convergence on a small full-rank system (isolates the temporal order) + @testset "second-order convergence" begin + Random.seed!(2) + Lc = 4 + Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lc)) + ψ_full = FiniteMPS(rand, ComplexF64, Lc, ℙ^2, ℙ^4) # full-rank: 1,2,4,2,1 + + Hmat = convert(TensorMap, Hc) + ψvec = convert(TensorMap, ψ_full) + ψvec /= norm(ψvec) + + T = 0.5 + dts = [0.1, 0.05, 0.025] + errs = map(dts) do δt + n = round(Int, T / δt) + ref = exp(-im * Hmat * (n * δt)) * ψvec + ψ = copy(ψ_full) + envs = environments(ψ, Hc, ψ) + for k in 0:(n - 1) + timestep!(ψ, Hc, k * δt, δt, BUG(), envs) + end + ψout = convert(TensorMap, ψ) + ψout /= norm(ψout) + return 1 - abs(dot(ψout, ref)) + end + + slopes = [ + (log(errs[i + 1]) - log(errs[i])) / (log(dts[i + 1]) - log(dts[i])) + for i in 1:(length(dts) - 1) + ] + @info "BUG convergence" errs slopes + for s in slopes + @test s ≈ 2 atol = 0.3 + end + end + + # 4. imaginary-time evolution lowers the energy toward the ground state (and, having no + # backward substep, stays norm-preserving/stable) + @testset "imaginary-time lowers energy" begin + Random.seed!(5) + ψi = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + E_start = real(expectation_value(ψi, H)) + E_prev = E_start + for _ in 1:20 + ψi, = timestep(ψi, H, 0.0, 0.1, BUG(); imaginary_evolution = true, normalize = true) + E_now = real(expectation_value(ψi, H)) + @test E_now ≤ E_prev + 1.0e-6 # monotone (non-increasing) energy + E_prev = E_now + end + @test E_prev < E_start - 1.0 # substantial lowering toward the ground state + @test norm(ψi) ≈ 1 atol = 1.0e-6 # imaginary-time BUG renormalizes each step + end + + # 5. LazySum / MultipliedOperator smoke tests + @testset "LazySum" begin + Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H]) + ψl, envs = timestep(ψ₀, Hlazy, 0.0, dt, BUG()) + E = expectation_value(ψl, Hlazy, envs) + @test (3 + 1.55 - 0.1) * E₀ ≈ E atol = 1.0e-2 + end + + @testset "TimeDependent LazySum" begin + Ht = MultipliedOperator(H, t -> 4) + MultipliedOperator(H, 1.45) + ψa, envsa = timestep(ψ₀, Ht(1.0), 0.0, dt, BUG()) + Ea = expectation_value(ψa, Ht(1.0), envsa) + + ψt, envst = timestep(ψ₀, Ht, 1.0, dt, BUG()) + Et = expectation_value(ψt, Ht(1.0), envst) + @test Ea ≈ Et atol = 1.0e-8 + end +end + +# Rank-adaptive BUG (Stage 2): a truncating `trscheme` enables basis augmentation + truncation, so +# the bond dimension grows and shrinks to track the entanglement of the evolving state. Trivial +# tensors only (symmetry stress is Chunk 2.3). The default `notrunc()` regression is covered by the +# fixed-rank testsets above (they must all still pass unchanged). +@testset "BUG rank-adaptive" verbose = true begin + # 1. bond growth: a tight tolerance grows a low-bond-dim state, a looser tolerance keeps it + # smaller. + @testset "bond growth" begin + Random.seed!(101) + L = 6 + H = force_planar(transverse_field_ising(ComplexF64, Trivial; L)) + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^2)) # low bond dim (2) + normalize!(ψ₀) + Dstart = maximum(dim(left_virtualspace(ψ₀, k)) for k in 1:L) + + ψtight = ψ₀ + for k in 0:2 + ψtight, = timestep(ψtight, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) + end + Dtight = maximum(dim(left_virtualspace(ψtight, k)) for k in 1:L) + + ψloose = ψ₀ + for k in 0:2 + ψloose, = timestep(ψloose, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-2))) + end + Dloose = maximum(dim(left_virtualspace(ψloose, k)) for k in 1:L) + + @info "BUG rank-adaptive bond growth" Dstart Dloose Dtight + @test Dtight > Dstart # rank-adaptivity grows the bond + @test Dloose < Dtight # a looser tolerance keeps a smaller bond + end + + # 2. THE HARD GATE: overlap-error vs the dense exp(-iHT) reference decreases (monotonically, up + # to the plateau at the fixed-dt floor) as the truncation tolerance ϑ shrinks. This proves the + # augmentation actually captures the true dynamics. + @testset "accuracy improves as ϑ decreases" begin + Random.seed!(202) + L = 6 + H = force_planar(transverse_field_ising(ComplexF64, Trivial; L)) + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^2)) # low-rank start the dynamics grows + normalize!(ψ₀) + + Hmat = convert(TensorMap, H) + ψvec = convert(TensorMap, ψ₀) + ψvec /= norm(ψvec) + + T = 0.2 + dt = 0.05 + n = round(Int, T / dt) + ref = exp(-im * Hmat * (n * dt)) * ψvec + + ϑs = [1.0e-2, 1.0e-4, 1.0e-6] + errs = map(ϑs) do ϑ + alg = BUG(; trscheme = truncerror(; atol = ϑ)) + ψ = copy(ψ₀) + envs = environments(ψ, H, ψ) + for k in 0:(n - 1) + timestep!(ψ, H, k * dt, dt, alg, envs) + end + ψout = convert(TensorMap, ψ) + ψout /= norm(ψout) + return 1 - abs(dot(ψout, ref)) + end + + @info "BUG rank-adaptive accuracy vs ϑ" ϑs errs + for i in 1:(length(ϑs) - 1) + @test errs[i + 1] ≤ 1.5 * errs[i] # monotone within plateau noise near the dt-floor + end + @test errs[end] < errs[1] / 10 # clear net improvement toward the dt-floor + end + + # 3. CBE-style comparison: from a low-rank state, rank-adaptive BUG tracks a bond-adaptive TDVP2 + # reference better than fixed-rank `BUG()` does (mirrors the CBE-TDVP test). + @testset "tracks TDVP2 better than fixed-rank BUG" begin + Random.seed!(303) + L = 8 + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + Dstart, Dcap, dt = 2, 16, 0.05 + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart)) + + ref, adaptive, fixed = ψ₀, ψ₀, ψ₀ + for _ in 1:6 + ref, = timestep(ref, H, 0.0, dt, TDVP2(; trscheme = truncrank(Dcap))) + adaptive, = timestep(adaptive, H, 0.0, dt, BUG(; trscheme = truncerror(; atol = 1.0e-8))) + fixed, = timestep(fixed, H, 0.0, dt, BUG()) + end + + @test dim(left_virtualspace(adaptive, L ÷ 2)) > Dstart # adaptive grew the bond + @test dim(left_virtualspace(fixed, L ÷ 2)) == Dstart # fixed-rank stuck at Dstart + @test abs(dot(ref, adaptive)) > abs(dot(ref, fixed)) # and tracks the reference better + end + + # 4. imaginary-time ground-state search: from a low bond dim, rank-adaptive imaginary-time BUG + # grows the bond and lowers the energy toward the true ground state. + @testset "imaginary-time grows bond and lowers energy" begin + Random.seed!(404) + L = 8 + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + ψgs, = find_groundstate(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^16), H; verbosity = 0) + Egs = real(expectation_value(ψgs, H)) + + ψ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^2)) # low bond dim + Dstart = maximum(dim(left_virtualspace(ψ, k)) for k in 1:L) + E_start = real(expectation_value(ψ, H)) + for _ in 1:30 + ψ, = timestep(ψ, H, 0.0, 0.1, BUG(; trscheme = truncerror(; atol = 1.0e-8)); imaginary_evolution = true, normalize = true) + end + Dend = maximum(dim(left_virtualspace(ψ, k)) for k in 1:L) + E_end = real(expectation_value(ψ, H)) + + @info "BUG rank-adaptive imaginary-time" Dstart Dend E_start E_end Egs + @test Dend > Dstart # the bond grew as entanglement built up + @test E_end < E_start - 1.0 # substantial lowering + @test E_end ≈ Egs atol = 0.6 # toward the true ground state (loose) + @test norm(ψ) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + end + + # 5. real-time energy conservation with truncation (loose atol) and norm behaviour: for a small + # tolerance the truncation is negligible, so the norm stays ≈ 1. + @testset "real-time energy conservation and norm" begin + Random.seed!(505) + L = 6 + H = force_planar(transverse_field_ising(ComplexF64, Trivial; L)) + ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + normalize!(ψ₀) + E₀ = real(expectation_value(ψ₀, H)) + + ψ = ψ₀ + for k in 0:4 + ψ, = timestep(ψ, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-8))) + end + @test real(expectation_value(ψ, H)) ≈ E₀ atol = 1.0e-2 # energy conserved (loose) + @test norm(ψ) ≈ 1 atol = 1.0e-6 # tiny ϑ ⇒ norm preserved + end +end + +# Charge-sector (symmetric-tensor) coverage for the fixed-rank BUG. These use *genuine* +# symmetric tensors (no `force_planar`), exercising the graded-bond paths flagged in the design +# doc's hsector risk register (H1/H6/H7): the transport-tensor seed `isomorphism(V ← V)`, the +# (co)domain/dual conventions of the inlined transport contractions, and the adjoints carrying +# sector duals. A fixed-rank step must preserve the total charge and the graded structure of every bond. +@testset "BUG symmetric tensors" verbose = true begin + dt = 0.1 + L = 6 + + # 1. U(1)-symmetric Heisenberg, both in the natural total-Sz = 0 sector and in a fixed nonzero + # total-charge (Sz = 1) sector: energy conservation + eigenstate phase + sector preservation. + @testset "U(1) Heisenberg (total Sz = $label)" for (label, right) in + (("0", U1Space(0 => 1)), ("1", U1Space(1 => 1))) + Random.seed!(2718) + H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψ = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]; right) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + Vl₀ = left_virtualspace.(Ref(ψ₀), 1:L) + Vr₀ = right_virtualspace.(Ref(ψ₀), 1:L) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) + E1 = expectation_value(ψ1, H, envs) + + @test E₀ ≈ E1 atol = 1.0e-2 + @test imag(E1) ≈ 0 atol = 1.0e-8 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + # the fixed-rank step preserves the graded structure (sector content) of every bond + @test left_virtualspace.(Ref(ψ1), 1:L) == Vl₀ + @test right_virtualspace.(Ref(ψ1), 1:L) == Vr₀ + end + + # 2. A second symmetry group. Z2 (transverse-field Ising) and SU2 (Heisenberg) both stress the + # graded transport tensor; same assertions (energy conservation + eigenstate phase). + @testset "Z2 transverse-field Ising" begin + Random.seed!(161803) + H = transverse_field_ising(ComplexF64, Z2Irrep; g = 1.0, L) + ψ = FiniteMPS(physicalspace(H), Z2Space(0 => 4, 1 => 4)) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + Vl₀ = left_virtualspace.(Ref(ψ₀), 1:L) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + @test left_virtualspace.(Ref(ψ1), 1:L) == Vl₀ + end + + @testset "SU(2) Heisenberg" begin + Random.seed!(577215) + H = heisenberg_XXX(ComplexF64, SU2Irrep; spin = 1 // 2, L) + # SU(2) spin-1/2 bonds alternate between integer / half-integer spins, so use the + # model's own full-rank virtual spaces rather than a hand-picked (integer-only) space. + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψ = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + Vl₀ = left_virtualspace.(Ref(ψ₀), 1:L) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + @test left_virtualspace.(Ref(ψ1), 1:L) == Vl₀ + end + + # 3. Imaginary-time symmetric evolution lowers the energy while preserving the sector + norm. + @testset "imaginary-time (U(1))" begin + Random.seed!(141421) + H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψi = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]) + Vl₀ = left_virtualspace.(Ref(ψi), 1:L) + Vr₀ = right_virtualspace.(Ref(ψi), 1:L) + E_start = real(expectation_value(ψi, H)) + E_prev = E_start + for _ in 1:15 + ψi, = timestep(ψi, H, 0.0, 0.1, BUG(); imaginary_evolution = true, normalize = true) + E_now = real(expectation_value(ψi, H)) + @test E_now ≤ E_prev + 1.0e-6 # monotone (non-increasing) energy + E_prev = E_now + end + @test E_prev < E_start - 0.5 # substantial lowering toward the ground state + @test norm(ψi) ≈ 1 atol = 1.0e-6 # imaginary-time BUG renormalizes each step + # the sector content of every bond is preserved throughout the imaginary-time sweep + @test left_virtualspace.(Ref(ψi), 1:L) == Vl₀ + @test right_virtualspace.(Ref(ψi), 1:L) == Vr₀ + end +end + +# Charge-sector RANK-ADAPTIVE coverage (Chunk 2.3): the "sector-adaptivity" path — genuine +# symmetric tensors under a *truncating* `trscheme`, so bonds grow and shrink per sector. This is +# the hardest part of the design doc; it stresses the H2–H5/H10 pitfalls of the risk register: +# augmentation as a per-sector direct sum (H2/H10), the "old basis first" invariant per sector +# (H3), a global-ϑ truncation dropping a sector to dimension 0 ⇒ dynamic bond grading (H4), and +# total-boundary-charge conservation through the zero-block embeddings (H5). All states start +# genuinely low-rank so the dynamics *must* grow the bonds; no `force_planar`. +@testset "BUG rank-adaptive symmetric" verbose = true begin + # A genuinely low-rank U(1) start. For spin-1/2 the virtual bonds alternate integer / + # half-integer parity, so a single-sector cap collapses the state (no fusion channels); use a + # mixed-parity cap with multiplicity 1 per sector. The full-rank Sz=0 profile is [1,2,4,8,4,2]; + # this cap gives [1,2,3,2,3,2], leaving room (esp. on the middle bond) for per-sector growth. + u1cap = U1Space(-1 // 2 => 1, 1 // 2 => 1, 0 => 1, 1 => 1, -1 => 1) + function low_rank_u1(L; seed = 2718, right = U1Space(0 => 1)) + Random.seed!(seed) + H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) + ψ = FiniteMPS(physicalspace(H), u1cap; right) + normalize!(ψ) + return H, ψ + end + maxbond(ψ) = maximum(dim(left_virtualspace(ψ, k)) for k in 1:length(ψ)) + secmults(V) = [c => dim(V, c) for c in sectors(V)] + + # 1. THE CORE GATE: rank-adaptivity under symmetry grows the bond per sector for a tight ϑ and + # keeps it small for a loose ϑ, without any `SpaceMismatch`, while the total boundary charge + # (the fixed `right` virtual space) is preserved. Checked in the natural total-Sz=0 sector and + # in a fixed nonzero-charge (Sz=1) sector. + @testset "rank-adaptivity + total-charge preservation (Sz = $label)" for (label, right) in + (("0", U1Space(0 => 1)), ("1", U1Space(1 => 1))) + L = 6 + H, ψ₀ = low_rank_u1(L; right) + Rtot = right_virtualspace(ψ₀, L) + Dstart = maxbond(ψ₀) + + ψtight = ψ₀ + for k in 0:2 + ψtight, = timestep(ψtight, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) + end + Dtight = maxbond(ψtight) + + ψloose = ψ₀ + for k in 0:2 + ψloose, = timestep(ψloose, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-2))) + end + Dloose = maxbond(ψloose) + + @info "BUG rank-adaptive symmetric bond growth (Sz=$label)" Dstart Dloose Dtight + @test Dtight > Dstart # per-sector augmentation grows the bond + @test Dloose < Dtight # a looser tolerance keeps a smaller bond + @test right_virtualspace(ψtight, L) == Rtot # H5: total boundary charge preserved + @test right_virtualspace(ψloose, L) == Rtot + @test norm(ψtight) ≈ 1 atol = 1.0e-6 # tiny ϑ ⇒ negligible truncation + end + + # 2. THE HARD GATE: overlap-error vs a *dense* `exp(-iH·T)` reference decreases as ϑ shrinks + # (monotonically up to the fixed-dt plateau). This proves per-sector augmentation+truncation + # actually captures the true dynamics under symmetry, not just that some bond grows. + @testset "accuracy improves as ϑ decreases (U(1))" begin + L = 6 + H, ψ₀ = low_rank_u1(L; seed = 202) + + Hmat = convert(TensorMap, H) + ψvec = convert(TensorMap, ψ₀) + ψvec /= norm(ψvec) + + T = 0.2 + dt = 0.05 + n = round(Int, T / dt) + ref = exp(-im * Hmat * (n * dt)) * ψvec + + ϑs = [1.0e-2, 1.0e-4, 1.0e-6] + errs = map(ϑs) do ϑ + ψ = copy(ψ₀) + for k in 0:(n - 1) + ψ, = timestep(ψ, H, k * dt, dt, BUG(; trscheme = truncerror(; atol = ϑ))) + end + ψout = convert(TensorMap, ψ) + ψout /= norm(ψout) + return 1 - abs(dot(ψout, ref)) + end + + @info "BUG rank-adaptive symmetric accuracy vs ϑ" ϑs errs + for i in 1:(length(ϑs) - 1) + @test errs[i + 1] ≤ 1.5 * errs[i] # monotone within plateau noise near the dt-floor + end + # Clear net improvement toward the dt-floor. This charge-capped low-rank system saturates + # its (small) per-sector bonds quickly: at ϑ ≲ 1e-4 no directions are truncated, so the two + # tightest tolerances land on the *identical* pure-dt-discretization floor (the ϑ=1e-4 and + # 1e-6 errors coincide to the digit). The net truncation-limited improvement is therefore + # modest (≈2.5×) rather than the orders of magnitude a full-rank trivial system shows. + @test errs[end] < errs[1] / 2 + end + + # 3. DYNAMIC BOND GRADING: an interior bond's per-sector multiplicities + # `[dim(V, c) for c in sectors(V)]` are time-dependent — they change across a single + # rank-adaptive step (the graded structure is not fixed), while the total charge is. + @testset "dynamic bond grading (per-sector multiplicities are time-dependent)" begin + L = 6 + H, ψ₀ = low_rank_u1(L) + mid = 4 + before = secmults(left_virtualspace(ψ₀, mid)) + ψ, = timestep(ψ₀, H, 0.0, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) + after = secmults(left_virtualspace(ψ, mid)) + @info "BUG rank-adaptive symmetric dynamic grading" before after + @test before != after # graded structure evolved + @test right_virtualspace(ψ, L) == right_virtualspace(ψ₀, L) # ... at fixed total charge + end + + # 4. H4 — a global-ϑ cut can truncate a whole sector to dimension 0 (removing it from the bond), + # and downstream rank-adaptive steps must tolerate the changed / asymmetric grading. + # + # NOTE on the sharpest H4 test (deterministic *drop-and-re-add* in one run): this is NOT + # achievable with single-site BUG, and asserting it would be wrong. The augmentation candidate + # `ACᵢ` in each half-sweep always carries the (already-truncated) bond as its *domain*, so + # `_bug_augment_left`/`_bug_augment_right` can only append directions whose charge sectors are + # already present on that bond — they grow multiplicity within existing sectors but can never + # re-introduce a sector that was truncated away (the graded analog of single-site TDVP's + # inability to change bond quantum numbers; re-adding a sector needs a two-site / CBE-style + # candidate, cf. `OptimalExpand`). We verified empirically that a sector dropped from an + # interior bond does not reappear over many subsequent tight steps. We therefore assert the + # deterministic *drop* and the H4 tolerance requirement (subsequent steps run without + # `SpaceMismatch`, yield a valid normalizable state, and preserve the total charge). + @testset "sector drop-to-zero + dynamic-grading tolerance (H4)" begin + L = 6 + H, ψ₀ = low_rank_u1(L) + # grow a rich interior bond that carries the subdominant ±1 sectors + ψrich = ψ₀ + for k in 0:5 + ψrich, = timestep(ψrich, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) + end + mid = 3 + Rtot = right_virtualspace(ψrich, L) + @test U1Irrep(-1) in sectors(left_virtualspace(ψrich, mid)) # subdominant sector present + + # a global rank cut pools singular values across all sectors under one threshold, so the + # subdominant sector is truncated to dimension 0 and removed from the bond (H4). + ψdrop, = timestep(ψrich, H, 0.0, 0.05, BUG(; trscheme = truncrank(2))) + @test !(U1Irrep(-1) in sectors(left_virtualspace(ψdrop, mid))) # dropped to dim 0 + @test right_virtualspace(ψdrop, L) == Rtot # ... charge still preserved + + # subsequent rank-adaptive steps must tolerate the reduced / asymmetric grading: no + # SpaceMismatch, a valid normalizable state, and a conserved total charge. + ψcont = ψdrop + for k in 0:3 + ψcont, = timestep(ψcont, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) + end + @info "BUG rank-adaptive symmetric H4" dropped = sectors(left_virtualspace(ψdrop, mid)) continued = sectors(left_virtualspace(ψcont, mid)) + @test isfinite(real(expectation_value(ψcont, H))) + @test norm(ψcont) > 0 + @test right_virtualspace(ψcont, L) == Rtot + end + + # 5. IMAGINARY-TIME symmetric rank-adaptive ground-state search: from a low bond dim it grows the + # per-sector bonds, lowers the energy toward `find_groundstate`, and preserves both the total + # charge and the (renormalized) norm. + @testset "imaginary-time symmetric grows bond and lowers energy" begin + L = 6 + H, ψ₀ = low_rank_u1(L) + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψgs, = find_groundstate(FiniteMPS(physicalspace(H), maxV[2:(end - 1)]), H; verbosity = 0) + Egs = real(expectation_value(ψgs, H)) + + Rtot = right_virtualspace(ψ₀, L) + Dstart = maxbond(ψ₀) + E_start = real(expectation_value(ψ₀, H)) + + ψ = ψ₀ + for _ in 1:30 + ψ, = timestep(ψ, H, 0.0, 0.1, BUG(; trscheme = truncerror(; atol = 1.0e-8)); imaginary_evolution = true, normalize = true) + end + Dend = maxbond(ψ) + E_end = real(expectation_value(ψ, H)) + + @info "BUG rank-adaptive symmetric imaginary-time" Dstart Dend E_start E_end Egs + @test Dend > Dstart # per-sector bonds grew as entanglement built up + @test E_end < E_start - 1.0 # substantial lowering + @test E_end ≈ Egs atol = 0.1 # toward the true ground state + @test norm(ψ) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + @test right_virtualspace(ψ, L) == Rtot # total charge preserved throughout + end +end diff --git a/test/algorithms/bug_augment.jl b/test/algorithms/bug_augment.jl new file mode 100644 index 000000000..947d2c378 --- /dev/null +++ b/test/algorithms/bug_augment.jl @@ -0,0 +1,136 @@ +println(" +------------------------------------ +| BUG basis-augmentation tests | +------------------------------------ +") + +using .TestSetup +using Test, TestExtras +using MPSKit +using MPSKit: _bug_augment_left, _bug_augment_right, _transpose_tail, _transpose_front, + left_orth, right_gauge +using TensorKit +using TensorKit: ℙ +using LinearAlgebra: I, norm +using Random + +# The augment helpers keep the OLD isometry `U₀` as the leading per-sector block and append the +# component of the evolved candidate `K₁` that is orthogonal to it (no truncation). The four core +# properties (checked below for both sweep directions, trivial + U(1)): +# 1. isometry `Û' Û ≈ 𝟙` +# 2. old-first, per sector `M = Û' U₀ ≈ [𝟙; 0]` block-by-block +# 3. range ⊇ old, candidate `Û (Û' U₀) ≈ U₀` and `Û (Û' K₁) ≈ K₁` +# 4. rank growth `dim(Vr₀) < dim(V̂) ≤ 2·dim(Vr₀)` + +# Assert the per-sector "old-first" invariant on the overlap `M` (codomain = augmented bond, +# domain = old bond `V_old`): every sector block must be `[𝟙; 0]`, i.e. the leading +# `dim(V_old, c)` rows are the identity and the rest vanish. Iterating `blocks(M)` visits exactly +# the sectors common to `M`'s (co)domain, so purely-new sectors (absent from `V_old`) are correctly +# skipped here. +function _check_old_first(M, V_old; tol = 1.0e-10) + for (c, b) in blocks(M) + r0 = dim(V_old, c) + @test b[1:r0, :] ≈ I + if r0 < size(b, 1) + @test norm(b[(r0 + 1):end, :]) < tol + end + end + return nothing +end + +function check_augment_left(U₀, K₁; tol = 1.0e-10) + Û = _bug_augment_left(U₀, K₁) + M = Û' * U₀ # old bond's coordinates in the augmented basis (V̂ ← Vr₀) + # 1. isometry + @test Û' * Û ≈ one(Û' * Û) + # 2. old-first, per sector (M : V̂ ← Vr₀) + _check_old_first(M, domain(U₀)[1]; tol) + # 3. range captures both the old basis and the evolved candidate + @test Û * (Û' * U₀) ≈ U₀ + @test Û * (Û' * K₁) ≈ K₁ + @test Û * M ≈ U₀ # consistency of the returned overlap + # 4. rank growth: strictly bigger than the old bond, at most doubled + r = dim(domain(U₀)) + @test r < dim(domain(Û)) ≤ 2r + return Û, M +end + +function check_augment_right(U₀, K₁; tol = 1.0e-10) + Û = _bug_augment_right(U₀, K₁) + ût = _transpose_tail(Û) # V̂ ← P ⊗ Vr, right-isometric (row space) + u0t = _transpose_tail(U₀) + k1t = _transpose_tail(K₁) + M = ût * u0t' # old bond's coordinates in the augmented basis (V̂ ← Vl₀) + # 1. isometry (right-canonical ⇒ tail has orthonormal rows) + @test ût * ût' ≈ one(ût * ût') + # 2. old-first, per sector (M : V̂ ← Vl₀) + _check_old_first(M, codomain(u0t)[1]; tol) + # 3. row space captures both the old basis and the evolved candidate + @test (u0t * ût') * ût ≈ u0t + @test (k1t * ût') * ût ≈ k1t + @test M' * ût ≈ u0t # consistency of the returned overlap + # 4. rank growth on the left bond + r = dim(codomain(u0t)) + @test r < dim(codomain(ût)) ≤ 2r + return Û, M +end + +@testset "BUG basis augmentation" verbose = true begin + # ------------------------------------------------------------------------------------------- + # Trivial (dense) tensors, ComplexF64 to exercise the adjoints. + # ------------------------------------------------------------------------------------------- + @testset "trivial tensors" begin + Random.seed!(20260707) + Vl = ℂ^2 + P = ℂ^3 + Vr = ℂ^2 + + # left→right: augment the RIGHT bond (domain of a left-isometry) + U₀_L, _ = left_orth(randn(ComplexF64, Vl ⊗ P ← Vr)) # Vl⊗P ← Vr, left-isometric + K₁_L = randn(ComplexF64, Vl ⊗ P ← Vr) # generic evolved candidate + Û_L, _ = check_augment_left(U₀_L, K₁_L) + + # right→left: augment the LEFT bond (codomain of a right-isometry) + _, U₀_R = right_gauge(randn(ComplexF64, Vl ⊗ P ← Vr)) # V ⊗ P ← Vr, right-isometric + K₁_R = randn(ComplexF64, Vl ⊗ P ← Vr) + Û_R, _ = check_augment_right(U₀_R, K₁_R) + end + + # ------------------------------------------------------------------------------------------- + # U(1)-symmetric tensors: per-sector direct sums, old-first sector-by-sector, and a genuinely + # new sector introduced by the candidate on the left augment. + # ------------------------------------------------------------------------------------------- + @testset "U(1) symmetric tensors" begin + Random.seed!(31415926) + + # left→right, WITH a new sector: `Vr₀` deliberately omits sector 0, which `Vl ⊗ P` + # contains and the candidate `K₁` populates ⇒ augmentation must add it to `V̂`. + Vl = U1Space(0 => 1, 1 => 1) + P = U1Space(0 => 1, 1 => 1) # Vl⊗P : 0=>1, 1=>2, 2=>1 + Vr0 = U1Space(1 => 1, 2 => 1) # omits sector 0 + Vr0K = U1Space(0 => 1, 1 => 1, 2 => 1) # candidate populates sector 0 + + U₀_L, _ = left_orth(randn(ComplexF64, Vl ⊗ P ← Vr0)) + @test domain(U₀_L)[1] == Vr0 # left_orth kept the full old bond + K₁_L = randn(ComplexF64, Vl ⊗ P ← Vr0K) + Û_L, _ = check_augment_left(U₀_L, K₁_L) + + # the augmented bond is a per-sector direct sum that gained sector 0 + V̂_L = domain(Û_L)[1] + @test !(U1Irrep(0) in sectors(Vr0)) + @test U1Irrep(0) in sectors(V̂_L) + # every old sector survives with at least its old multiplicity (old-first ⇒ nothing dropped) + for c in sectors(Vr0) + @test dim(V̂_L, c) ≥ dim(Vr0, c) + end + + # right→left augment on a compatible symmetric configuration + Plr = U1Space(0 => 1, 1 => 1) + Vrr = U1Space(0 => 1, 1 => 1) + Vl0 = U1Space(0 => 1, 1 => 1) + _, U₀_R = right_gauge(randn(ComplexF64, Vl0 ⊗ Plr ← Vrr)) + Vl_K = U1Space(0 => 1, 1 => 1) + K₁_R = randn(ComplexF64, Vl_K ⊗ Plr ← Vrr) + Û_R, _ = check_augment_right(U₀_R, K₁_R) + end +end From bdf7cada23e20ddacc12e46f8d57584b42a7af25 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 14 Jul 2026 18:56:08 -0400 Subject: [PATCH 4/6] =?UTF-8?q?BUG:=20review=20cleanup=20=E2=80=94=20drop?= =?UTF-8?q?=20dead=20fields,=20doc=20fixes,=20changelog?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Remove unused `tolgauge`/`gaugemaxiter` fields from the `BUG` struct and constructor: they are read only for `InfiniteMPS` gauge reconstruction in TDVP, and BUG is finite-only, so they were silently ignored. - Clarify the `BUG` docstring (single-site, symmetric second-order, finite MPS) and correct the `finalize` signature docs (`finalize(t, …)`, not `iter`) in BUG, TDVP and TDVP2. - Document the `normalize` behavior change under changelog "Changed": it is now an explicit kwarg (default false) decoupled from `imaginary_evolution`; imaginary-time evolution no longer renormalizes unless `normalize = true`. Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/changelog.md | 6 ++++++ src/algorithms/timestep/bug.jl | 18 ++++++------------ src/algorithms/timestep/tdvp.jl | 4 ++-- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index c1bc75339..7cbcba8cc 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -28,6 +28,12 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### 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: diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl index c5f5724c3..4a83a62b1 100644 --- a/src/algorithms/timestep/bug.jl +++ b/src/algorithms/timestep/bug.jl @@ -1,8 +1,9 @@ """ $(TYPEDEF) -Single site MPS time-evolution algorithm based on the Basis-Update & Galerkin (BUG) integrator, -an unconventional robust integrator for dynamical low-rank approximation. +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 @@ -32,12 +33,6 @@ struct BUG{A, O, T, S, F} <: Algorithm "algorithm used in the exponential solvers" integrator::A - "tolerance for gauging algorithm" - tolgauge::Float64 - - "maximal amount of iterations for gauging algorithm" - gaugemaxiter::Int - "algorithm used to re-orthonormalize the basis after each local update" alg_orth::O @@ -47,16 +42,15 @@ struct BUG{A, O, T, S, F} <: Algorithm "algorithm used for the singular value decomposition" alg_svd::S - "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" + "callback function applied after each iteration, of signature `finalize(t, ψ, H, envs) -> ψ, envs`" finalize::F end function BUG(; - integrator = Defaults.alg_expsolve(), tolgauge = Defaults.tolgauge, - gaugemaxiter = Defaults.maxiter, alg_orth = Defaults.alg_orth(), + integrator = Defaults.alg_expsolve(), alg_orth = Defaults.alg_orth(), trscheme = notrunc(), alg_svd = Defaults.alg_svd(), finalize = Defaults._finalize ) - return BUG(integrator, tolgauge, gaugemaxiter, alg_orth, trscheme, alg_svd, finalize) + return BUG(integrator, alg_orth, trscheme, alg_svd, finalize) end function timestep!( diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index afaf3ba4e..44b39a54b 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -41,7 +41,7 @@ struct TDVP{A, E, G, F} <: Algorithm "factorization used for the post-update gauge: a QR algorithm (no truncation) or a truncated SVD" alg_gauge::G - "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" + "callback function applied after each iteration, of signature `finalize(t, ψ, H, envs) -> ψ, envs`" finalize::F end function TDVP(; @@ -218,7 +218,7 @@ $(TYPEDFIELDS) "algorithm used for truncation of the two-site update" trscheme::TruncationStrategy - "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" + "callback function applied after each iteration, of signature `finalize(t, ψ, H, envs) -> ψ, envs`" finalize::F = Defaults._finalize end From e312c277f7187b55fab6a1dce5c306b9d0de137e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 20 Jul 2026 19:44:17 -0400 Subject: [PATCH 5/6] rework tests --- test/algorithms/bug.jl | 527 --------------------------------- test/algorithms/bug_augment.jl | 14 +- test/algorithms/timestep.jl | 178 ++++++++++- 3 files changed, 181 insertions(+), 538 deletions(-) delete mode 100644 test/algorithms/bug.jl diff --git a/test/algorithms/bug.jl b/test/algorithms/bug.jl deleted file mode 100644 index d635ebd66..000000000 --- a/test/algorithms/bug.jl +++ /dev/null @@ -1,527 +0,0 @@ -println(" ------------------------------ -| BUG time-stepping tests | ------------------------------ -") - -using .TestSetup -using Test, TestExtras -using MPSKit -using TensorKit -using TensorKit: ℙ -using LinearAlgebra: dot, norm -using Random - -@testset "BUG time evolution" verbose = true begin - dt = 0.1 - L = 10 - - H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) - ψ = FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4) - ψ₀, = find_groundstate(ψ, H; verbosity = 0) - E₀ = expectation_value(ψ₀, H) - - # 1. energy conservation + eigenstate phase - @testset "energy conservation" begin - ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) - E1 = expectation_value(ψ1, H, envs) - @test E₀ ≈ E1 atol = 1.0e-2 - @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 - end - - # 2. agreement with TDVP over a few real-time steps of a random MPS - @testset "agreement with TDVP" begin - Random.seed!(1234) - ψr = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) - δt = 0.01 - ψ_bug, ψ_tdvp = ψr, ψr - for k in 0:4 - ψ_bug, = timestep(ψ_bug, H, k * δt, δt, BUG()) - ψ_tdvp, = timestep(ψ_tdvp, H, k * δt, δt, TDVP()) - end - @test expectation_value(ψ_bug, H) ≈ expectation_value(ψ_tdvp, H) atol = 1.0e-3 - @test abs(dot(ψ_bug, ψ_tdvp)) ≈ 1 atol = 1.0e-3 - end - - # 3. second-order convergence on a small full-rank system (isolates the temporal order) - @testset "second-order convergence" begin - Random.seed!(2) - Lc = 4 - Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lc)) - ψ_full = FiniteMPS(rand, ComplexF64, Lc, ℙ^2, ℙ^4) # full-rank: 1,2,4,2,1 - - Hmat = convert(TensorMap, Hc) - ψvec = convert(TensorMap, ψ_full) - ψvec /= norm(ψvec) - - T = 0.5 - dts = [0.1, 0.05, 0.025] - errs = map(dts) do δt - n = round(Int, T / δt) - ref = exp(-im * Hmat * (n * δt)) * ψvec - ψ = copy(ψ_full) - envs = environments(ψ, Hc, ψ) - for k in 0:(n - 1) - timestep!(ψ, Hc, k * δt, δt, BUG(), envs) - end - ψout = convert(TensorMap, ψ) - ψout /= norm(ψout) - return 1 - abs(dot(ψout, ref)) - end - - slopes = [ - (log(errs[i + 1]) - log(errs[i])) / (log(dts[i + 1]) - log(dts[i])) - for i in 1:(length(dts) - 1) - ] - @info "BUG convergence" errs slopes - for s in slopes - @test s ≈ 2 atol = 0.3 - end - end - - # 4. imaginary-time evolution lowers the energy toward the ground state (and, having no - # backward substep, stays norm-preserving/stable) - @testset "imaginary-time lowers energy" begin - Random.seed!(5) - ψi = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) - E_start = real(expectation_value(ψi, H)) - E_prev = E_start - for _ in 1:20 - ψi, = timestep(ψi, H, 0.0, 0.1, BUG(); imaginary_evolution = true, normalize = true) - E_now = real(expectation_value(ψi, H)) - @test E_now ≤ E_prev + 1.0e-6 # monotone (non-increasing) energy - E_prev = E_now - end - @test E_prev < E_start - 1.0 # substantial lowering toward the ground state - @test norm(ψi) ≈ 1 atol = 1.0e-6 # imaginary-time BUG renormalizes each step - end - - # 5. LazySum / MultipliedOperator smoke tests - @testset "LazySum" begin - Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H]) - ψl, envs = timestep(ψ₀, Hlazy, 0.0, dt, BUG()) - E = expectation_value(ψl, Hlazy, envs) - @test (3 + 1.55 - 0.1) * E₀ ≈ E atol = 1.0e-2 - end - - @testset "TimeDependent LazySum" begin - Ht = MultipliedOperator(H, t -> 4) + MultipliedOperator(H, 1.45) - ψa, envsa = timestep(ψ₀, Ht(1.0), 0.0, dt, BUG()) - Ea = expectation_value(ψa, Ht(1.0), envsa) - - ψt, envst = timestep(ψ₀, Ht, 1.0, dt, BUG()) - Et = expectation_value(ψt, Ht(1.0), envst) - @test Ea ≈ Et atol = 1.0e-8 - end -end - -# Rank-adaptive BUG (Stage 2): a truncating `trscheme` enables basis augmentation + truncation, so -# the bond dimension grows and shrinks to track the entanglement of the evolving state. Trivial -# tensors only (symmetry stress is Chunk 2.3). The default `notrunc()` regression is covered by the -# fixed-rank testsets above (they must all still pass unchanged). -@testset "BUG rank-adaptive" verbose = true begin - # 1. bond growth: a tight tolerance grows a low-bond-dim state, a looser tolerance keeps it - # smaller. - @testset "bond growth" begin - Random.seed!(101) - L = 6 - H = force_planar(transverse_field_ising(ComplexF64, Trivial; L)) - ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^2)) # low bond dim (2) - normalize!(ψ₀) - Dstart = maximum(dim(left_virtualspace(ψ₀, k)) for k in 1:L) - - ψtight = ψ₀ - for k in 0:2 - ψtight, = timestep(ψtight, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) - end - Dtight = maximum(dim(left_virtualspace(ψtight, k)) for k in 1:L) - - ψloose = ψ₀ - for k in 0:2 - ψloose, = timestep(ψloose, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-2))) - end - Dloose = maximum(dim(left_virtualspace(ψloose, k)) for k in 1:L) - - @info "BUG rank-adaptive bond growth" Dstart Dloose Dtight - @test Dtight > Dstart # rank-adaptivity grows the bond - @test Dloose < Dtight # a looser tolerance keeps a smaller bond - end - - # 2. THE HARD GATE: overlap-error vs the dense exp(-iHT) reference decreases (monotonically, up - # to the plateau at the fixed-dt floor) as the truncation tolerance ϑ shrinks. This proves the - # augmentation actually captures the true dynamics. - @testset "accuracy improves as ϑ decreases" begin - Random.seed!(202) - L = 6 - H = force_planar(transverse_field_ising(ComplexF64, Trivial; L)) - ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^2)) # low-rank start the dynamics grows - normalize!(ψ₀) - - Hmat = convert(TensorMap, H) - ψvec = convert(TensorMap, ψ₀) - ψvec /= norm(ψvec) - - T = 0.2 - dt = 0.05 - n = round(Int, T / dt) - ref = exp(-im * Hmat * (n * dt)) * ψvec - - ϑs = [1.0e-2, 1.0e-4, 1.0e-6] - errs = map(ϑs) do ϑ - alg = BUG(; trscheme = truncerror(; atol = ϑ)) - ψ = copy(ψ₀) - envs = environments(ψ, H, ψ) - for k in 0:(n - 1) - timestep!(ψ, H, k * dt, dt, alg, envs) - end - ψout = convert(TensorMap, ψ) - ψout /= norm(ψout) - return 1 - abs(dot(ψout, ref)) - end - - @info "BUG rank-adaptive accuracy vs ϑ" ϑs errs - for i in 1:(length(ϑs) - 1) - @test errs[i + 1] ≤ 1.5 * errs[i] # monotone within plateau noise near the dt-floor - end - @test errs[end] < errs[1] / 10 # clear net improvement toward the dt-floor - end - - # 3. CBE-style comparison: from a low-rank state, rank-adaptive BUG tracks a bond-adaptive TDVP2 - # reference better than fixed-rank `BUG()` does (mirrors the CBE-TDVP test). - @testset "tracks TDVP2 better than fixed-rank BUG" begin - Random.seed!(303) - L = 8 - H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) - Dstart, Dcap, dt = 2, 16, 0.05 - ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^Dstart)) - - ref, adaptive, fixed = ψ₀, ψ₀, ψ₀ - for _ in 1:6 - ref, = timestep(ref, H, 0.0, dt, TDVP2(; trscheme = truncrank(Dcap))) - adaptive, = timestep(adaptive, H, 0.0, dt, BUG(; trscheme = truncerror(; atol = 1.0e-8))) - fixed, = timestep(fixed, H, 0.0, dt, BUG()) - end - - @test dim(left_virtualspace(adaptive, L ÷ 2)) > Dstart # adaptive grew the bond - @test dim(left_virtualspace(fixed, L ÷ 2)) == Dstart # fixed-rank stuck at Dstart - @test abs(dot(ref, adaptive)) > abs(dot(ref, fixed)) # and tracks the reference better - end - - # 4. imaginary-time ground-state search: from a low bond dim, rank-adaptive imaginary-time BUG - # grows the bond and lowers the energy toward the true ground state. - @testset "imaginary-time grows bond and lowers energy" begin - Random.seed!(404) - L = 8 - H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) - ψgs, = find_groundstate(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^16), H; verbosity = 0) - Egs = real(expectation_value(ψgs, H)) - - ψ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^2)) # low bond dim - Dstart = maximum(dim(left_virtualspace(ψ, k)) for k in 1:L) - E_start = real(expectation_value(ψ, H)) - for _ in 1:30 - ψ, = timestep(ψ, H, 0.0, 0.1, BUG(; trscheme = truncerror(; atol = 1.0e-8)); imaginary_evolution = true, normalize = true) - end - Dend = maximum(dim(left_virtualspace(ψ, k)) for k in 1:L) - E_end = real(expectation_value(ψ, H)) - - @info "BUG rank-adaptive imaginary-time" Dstart Dend E_start E_end Egs - @test Dend > Dstart # the bond grew as entanglement built up - @test E_end < E_start - 1.0 # substantial lowering - @test E_end ≈ Egs atol = 0.6 # toward the true ground state (loose) - @test norm(ψ) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step - end - - # 5. real-time energy conservation with truncation (loose atol) and norm behaviour: for a small - # tolerance the truncation is negligible, so the norm stays ≈ 1. - @testset "real-time energy conservation and norm" begin - Random.seed!(505) - L = 6 - H = force_planar(transverse_field_ising(ComplexF64, Trivial; L)) - ψ₀ = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) - normalize!(ψ₀) - E₀ = real(expectation_value(ψ₀, H)) - - ψ = ψ₀ - for k in 0:4 - ψ, = timestep(ψ, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-8))) - end - @test real(expectation_value(ψ, H)) ≈ E₀ atol = 1.0e-2 # energy conserved (loose) - @test norm(ψ) ≈ 1 atol = 1.0e-6 # tiny ϑ ⇒ norm preserved - end -end - -# Charge-sector (symmetric-tensor) coverage for the fixed-rank BUG. These use *genuine* -# symmetric tensors (no `force_planar`), exercising the graded-bond paths flagged in the design -# doc's hsector risk register (H1/H6/H7): the transport-tensor seed `isomorphism(V ← V)`, the -# (co)domain/dual conventions of the inlined transport contractions, and the adjoints carrying -# sector duals. A fixed-rank step must preserve the total charge and the graded structure of every bond. -@testset "BUG symmetric tensors" verbose = true begin - dt = 0.1 - L = 6 - - # 1. U(1)-symmetric Heisenberg, both in the natural total-Sz = 0 sector and in a fixed nonzero - # total-charge (Sz = 1) sector: energy conservation + eigenstate phase + sector preservation. - @testset "U(1) Heisenberg (total Sz = $label)" for (label, right) in - (("0", U1Space(0 => 1)), ("1", U1Space(1 => 1))) - Random.seed!(2718) - H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) - maxV = MPSKit.max_virtualspaces(physicalspace(H)) - ψ = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]; right) - ψ₀, = find_groundstate(ψ, H; verbosity = 0) - E₀ = expectation_value(ψ₀, H) - - Vl₀ = left_virtualspace.(Ref(ψ₀), 1:L) - Vr₀ = right_virtualspace.(Ref(ψ₀), 1:L) - - ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) - E1 = expectation_value(ψ1, H, envs) - - @test E₀ ≈ E1 atol = 1.0e-2 - @test imag(E1) ≈ 0 atol = 1.0e-8 - @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 - # the fixed-rank step preserves the graded structure (sector content) of every bond - @test left_virtualspace.(Ref(ψ1), 1:L) == Vl₀ - @test right_virtualspace.(Ref(ψ1), 1:L) == Vr₀ - end - - # 2. A second symmetry group. Z2 (transverse-field Ising) and SU2 (Heisenberg) both stress the - # graded transport tensor; same assertions (energy conservation + eigenstate phase). - @testset "Z2 transverse-field Ising" begin - Random.seed!(161803) - H = transverse_field_ising(ComplexF64, Z2Irrep; g = 1.0, L) - ψ = FiniteMPS(physicalspace(H), Z2Space(0 => 4, 1 => 4)) - ψ₀, = find_groundstate(ψ, H; verbosity = 0) - E₀ = expectation_value(ψ₀, H) - Vl₀ = left_virtualspace.(Ref(ψ₀), 1:L) - - ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) - E1 = expectation_value(ψ1, H, envs) - @test E₀ ≈ E1 atol = 1.0e-2 - @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 - @test left_virtualspace.(Ref(ψ1), 1:L) == Vl₀ - end - - @testset "SU(2) Heisenberg" begin - Random.seed!(577215) - H = heisenberg_XXX(ComplexF64, SU2Irrep; spin = 1 // 2, L) - # SU(2) spin-1/2 bonds alternate between integer / half-integer spins, so use the - # model's own full-rank virtual spaces rather than a hand-picked (integer-only) space. - maxV = MPSKit.max_virtualspaces(physicalspace(H)) - ψ = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]) - ψ₀, = find_groundstate(ψ, H; verbosity = 0) - E₀ = expectation_value(ψ₀, H) - Vl₀ = left_virtualspace.(Ref(ψ₀), 1:L) - - ψ1, envs = timestep(ψ₀, H, 0.0, dt, BUG()) - E1 = expectation_value(ψ1, H, envs) - @test E₀ ≈ E1 atol = 1.0e-2 - @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 - @test left_virtualspace.(Ref(ψ1), 1:L) == Vl₀ - end - - # 3. Imaginary-time symmetric evolution lowers the energy while preserving the sector + norm. - @testset "imaginary-time (U(1))" begin - Random.seed!(141421) - H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) - maxV = MPSKit.max_virtualspaces(physicalspace(H)) - ψi = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]) - Vl₀ = left_virtualspace.(Ref(ψi), 1:L) - Vr₀ = right_virtualspace.(Ref(ψi), 1:L) - E_start = real(expectation_value(ψi, H)) - E_prev = E_start - for _ in 1:15 - ψi, = timestep(ψi, H, 0.0, 0.1, BUG(); imaginary_evolution = true, normalize = true) - E_now = real(expectation_value(ψi, H)) - @test E_now ≤ E_prev + 1.0e-6 # monotone (non-increasing) energy - E_prev = E_now - end - @test E_prev < E_start - 0.5 # substantial lowering toward the ground state - @test norm(ψi) ≈ 1 atol = 1.0e-6 # imaginary-time BUG renormalizes each step - # the sector content of every bond is preserved throughout the imaginary-time sweep - @test left_virtualspace.(Ref(ψi), 1:L) == Vl₀ - @test right_virtualspace.(Ref(ψi), 1:L) == Vr₀ - end -end - -# Charge-sector RANK-ADAPTIVE coverage (Chunk 2.3): the "sector-adaptivity" path — genuine -# symmetric tensors under a *truncating* `trscheme`, so bonds grow and shrink per sector. This is -# the hardest part of the design doc; it stresses the H2–H5/H10 pitfalls of the risk register: -# augmentation as a per-sector direct sum (H2/H10), the "old basis first" invariant per sector -# (H3), a global-ϑ truncation dropping a sector to dimension 0 ⇒ dynamic bond grading (H4), and -# total-boundary-charge conservation through the zero-block embeddings (H5). All states start -# genuinely low-rank so the dynamics *must* grow the bonds; no `force_planar`. -@testset "BUG rank-adaptive symmetric" verbose = true begin - # A genuinely low-rank U(1) start. For spin-1/2 the virtual bonds alternate integer / - # half-integer parity, so a single-sector cap collapses the state (no fusion channels); use a - # mixed-parity cap with multiplicity 1 per sector. The full-rank Sz=0 profile is [1,2,4,8,4,2]; - # this cap gives [1,2,3,2,3,2], leaving room (esp. on the middle bond) for per-sector growth. - u1cap = U1Space(-1 // 2 => 1, 1 // 2 => 1, 0 => 1, 1 => 1, -1 => 1) - function low_rank_u1(L; seed = 2718, right = U1Space(0 => 1)) - Random.seed!(seed) - H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) - ψ = FiniteMPS(physicalspace(H), u1cap; right) - normalize!(ψ) - return H, ψ - end - maxbond(ψ) = maximum(dim(left_virtualspace(ψ, k)) for k in 1:length(ψ)) - secmults(V) = [c => dim(V, c) for c in sectors(V)] - - # 1. THE CORE GATE: rank-adaptivity under symmetry grows the bond per sector for a tight ϑ and - # keeps it small for a loose ϑ, without any `SpaceMismatch`, while the total boundary charge - # (the fixed `right` virtual space) is preserved. Checked in the natural total-Sz=0 sector and - # in a fixed nonzero-charge (Sz=1) sector. - @testset "rank-adaptivity + total-charge preservation (Sz = $label)" for (label, right) in - (("0", U1Space(0 => 1)), ("1", U1Space(1 => 1))) - L = 6 - H, ψ₀ = low_rank_u1(L; right) - Rtot = right_virtualspace(ψ₀, L) - Dstart = maxbond(ψ₀) - - ψtight = ψ₀ - for k in 0:2 - ψtight, = timestep(ψtight, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) - end - Dtight = maxbond(ψtight) - - ψloose = ψ₀ - for k in 0:2 - ψloose, = timestep(ψloose, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-2))) - end - Dloose = maxbond(ψloose) - - @info "BUG rank-adaptive symmetric bond growth (Sz=$label)" Dstart Dloose Dtight - @test Dtight > Dstart # per-sector augmentation grows the bond - @test Dloose < Dtight # a looser tolerance keeps a smaller bond - @test right_virtualspace(ψtight, L) == Rtot # H5: total boundary charge preserved - @test right_virtualspace(ψloose, L) == Rtot - @test norm(ψtight) ≈ 1 atol = 1.0e-6 # tiny ϑ ⇒ negligible truncation - end - - # 2. THE HARD GATE: overlap-error vs a *dense* `exp(-iH·T)` reference decreases as ϑ shrinks - # (monotonically up to the fixed-dt plateau). This proves per-sector augmentation+truncation - # actually captures the true dynamics under symmetry, not just that some bond grows. - @testset "accuracy improves as ϑ decreases (U(1))" begin - L = 6 - H, ψ₀ = low_rank_u1(L; seed = 202) - - Hmat = convert(TensorMap, H) - ψvec = convert(TensorMap, ψ₀) - ψvec /= norm(ψvec) - - T = 0.2 - dt = 0.05 - n = round(Int, T / dt) - ref = exp(-im * Hmat * (n * dt)) * ψvec - - ϑs = [1.0e-2, 1.0e-4, 1.0e-6] - errs = map(ϑs) do ϑ - ψ = copy(ψ₀) - for k in 0:(n - 1) - ψ, = timestep(ψ, H, k * dt, dt, BUG(; trscheme = truncerror(; atol = ϑ))) - end - ψout = convert(TensorMap, ψ) - ψout /= norm(ψout) - return 1 - abs(dot(ψout, ref)) - end - - @info "BUG rank-adaptive symmetric accuracy vs ϑ" ϑs errs - for i in 1:(length(ϑs) - 1) - @test errs[i + 1] ≤ 1.5 * errs[i] # monotone within plateau noise near the dt-floor - end - # Clear net improvement toward the dt-floor. This charge-capped low-rank system saturates - # its (small) per-sector bonds quickly: at ϑ ≲ 1e-4 no directions are truncated, so the two - # tightest tolerances land on the *identical* pure-dt-discretization floor (the ϑ=1e-4 and - # 1e-6 errors coincide to the digit). The net truncation-limited improvement is therefore - # modest (≈2.5×) rather than the orders of magnitude a full-rank trivial system shows. - @test errs[end] < errs[1] / 2 - end - - # 3. DYNAMIC BOND GRADING: an interior bond's per-sector multiplicities - # `[dim(V, c) for c in sectors(V)]` are time-dependent — they change across a single - # rank-adaptive step (the graded structure is not fixed), while the total charge is. - @testset "dynamic bond grading (per-sector multiplicities are time-dependent)" begin - L = 6 - H, ψ₀ = low_rank_u1(L) - mid = 4 - before = secmults(left_virtualspace(ψ₀, mid)) - ψ, = timestep(ψ₀, H, 0.0, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) - after = secmults(left_virtualspace(ψ, mid)) - @info "BUG rank-adaptive symmetric dynamic grading" before after - @test before != after # graded structure evolved - @test right_virtualspace(ψ, L) == right_virtualspace(ψ₀, L) # ... at fixed total charge - end - - # 4. H4 — a global-ϑ cut can truncate a whole sector to dimension 0 (removing it from the bond), - # and downstream rank-adaptive steps must tolerate the changed / asymmetric grading. - # - # NOTE on the sharpest H4 test (deterministic *drop-and-re-add* in one run): this is NOT - # achievable with single-site BUG, and asserting it would be wrong. The augmentation candidate - # `ACᵢ` in each half-sweep always carries the (already-truncated) bond as its *domain*, so - # `_bug_augment_left`/`_bug_augment_right` can only append directions whose charge sectors are - # already present on that bond — they grow multiplicity within existing sectors but can never - # re-introduce a sector that was truncated away (the graded analog of single-site TDVP's - # inability to change bond quantum numbers; re-adding a sector needs a two-site / CBE-style - # candidate, cf. `OptimalExpand`). We verified empirically that a sector dropped from an - # interior bond does not reappear over many subsequent tight steps. We therefore assert the - # deterministic *drop* and the H4 tolerance requirement (subsequent steps run without - # `SpaceMismatch`, yield a valid normalizable state, and preserve the total charge). - @testset "sector drop-to-zero + dynamic-grading tolerance (H4)" begin - L = 6 - H, ψ₀ = low_rank_u1(L) - # grow a rich interior bond that carries the subdominant ±1 sectors - ψrich = ψ₀ - for k in 0:5 - ψrich, = timestep(ψrich, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) - end - mid = 3 - Rtot = right_virtualspace(ψrich, L) - @test U1Irrep(-1) in sectors(left_virtualspace(ψrich, mid)) # subdominant sector present - - # a global rank cut pools singular values across all sectors under one threshold, so the - # subdominant sector is truncated to dimension 0 and removed from the bond (H4). - ψdrop, = timestep(ψrich, H, 0.0, 0.05, BUG(; trscheme = truncrank(2))) - @test !(U1Irrep(-1) in sectors(left_virtualspace(ψdrop, mid))) # dropped to dim 0 - @test right_virtualspace(ψdrop, L) == Rtot # ... charge still preserved - - # subsequent rank-adaptive steps must tolerate the reduced / asymmetric grading: no - # SpaceMismatch, a valid normalizable state, and a conserved total charge. - ψcont = ψdrop - for k in 0:3 - ψcont, = timestep(ψcont, H, k * 0.05, 0.05, BUG(; trscheme = truncerror(; atol = 1.0e-10))) - end - @info "BUG rank-adaptive symmetric H4" dropped = sectors(left_virtualspace(ψdrop, mid)) continued = sectors(left_virtualspace(ψcont, mid)) - @test isfinite(real(expectation_value(ψcont, H))) - @test norm(ψcont) > 0 - @test right_virtualspace(ψcont, L) == Rtot - end - - # 5. IMAGINARY-TIME symmetric rank-adaptive ground-state search: from a low bond dim it grows the - # per-sector bonds, lowers the energy toward `find_groundstate`, and preserves both the total - # charge and the (renormalized) norm. - @testset "imaginary-time symmetric grows bond and lowers energy" begin - L = 6 - H, ψ₀ = low_rank_u1(L) - maxV = MPSKit.max_virtualspaces(physicalspace(H)) - ψgs, = find_groundstate(FiniteMPS(physicalspace(H), maxV[2:(end - 1)]), H; verbosity = 0) - Egs = real(expectation_value(ψgs, H)) - - Rtot = right_virtualspace(ψ₀, L) - Dstart = maxbond(ψ₀) - E_start = real(expectation_value(ψ₀, H)) - - ψ = ψ₀ - for _ in 1:30 - ψ, = timestep(ψ, H, 0.0, 0.1, BUG(; trscheme = truncerror(; atol = 1.0e-8)); imaginary_evolution = true, normalize = true) - end - Dend = maxbond(ψ) - E_end = real(expectation_value(ψ, H)) - - @info "BUG rank-adaptive symmetric imaginary-time" Dstart Dend E_start E_end Egs - @test Dend > Dstart # per-sector bonds grew as entanglement built up - @test E_end < E_start - 1.0 # substantial lowering - @test E_end ≈ Egs atol = 0.1 # toward the true ground state - @test norm(ψ) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step - @test right_virtualspace(ψ, L) == Rtot # total charge preserved throughout - end -end diff --git a/test/algorithms/bug_augment.jl b/test/algorithms/bug_augment.jl index 947d2c378..bc1996f8f 100644 --- a/test/algorithms/bug_augment.jl +++ b/test/algorithms/bug_augment.jl @@ -7,13 +7,21 @@ println(" using .TestSetup using Test, TestExtras using MPSKit -using MPSKit: _bug_augment_left, _bug_augment_right, _transpose_tail, _transpose_front, - left_orth, right_gauge +using MPSKit: _transpose_tail, _transpose_front, left_orth, right_orth, right_gauge, Defaults using TensorKit -using TensorKit: ℙ +using TensorKit: ℙ, catdomain, catcodomain using LinearAlgebra: I, norm using Random +# The basis-augmentation step is inlined into `timestep!`; these local one-liners mirror it so the +# invariants can be exercised directly. Left: orthonormalize the stacked `[U₀ │ K₁]` (old first); +# right: the `_transpose_tail` mirror (LQ on the stacked rows `[U₀; K₁]`). +_bug_augment_left(U₀, K₁) = first(left_orth(catdomain(U₀, K₁); alg = Defaults.alg_orth())) +function _bug_augment_right(U₀, K₁) + _, Û = right_orth(catcodomain(_transpose_tail(U₀), _transpose_tail(K₁)); alg = Defaults.alg_orth()) + return _transpose_front(Û) +end + # The augment helpers keep the OLD isometry `U₀` as the leading per-sector block and append the # component of the evolved candidate `K₁` that is orthogonal to it (no truncation). The four core # properties (checked below for both sweep directions, trivial + U(1)): diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index ede3f5b9c..d72063e4b 100644 --- a/test/algorithms/timestep.jl +++ b/test/algorithms/timestep.jl @@ -9,15 +9,18 @@ using Test, TestExtras using MPSKit using TensorKit using TensorKit: ℙ -using LinearAlgebra: norm +using LinearAlgebra: dot, norm using Random verbosity_full = 5 verbosity_conv = 1 +# name a time-evolution algorithm for @testset labels ("TDVP", "TDVP2", "BUG") +algname(alg) = string(nameof(typeof(alg))) + @testset "timestep" verbose = true begin dt = 0.1 - algs = [TDVP(), TDVP2(; trscheme = truncrank(10))] + algs = [TDVP(), TDVP2(; trscheme = truncrank(10)), BUG()] L = 10 H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) @@ -26,7 +29,7 @@ verbosity_conv = 1 ψ₀, = find_groundstate(ψ, H) E₀ = expectation_value(ψ₀, H) - @testset "Finite $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + @testset "Finite $(algname(alg))" for alg in algs ψ1, envs = timestep(ψ₀, H, 0.0, dt, alg) E1 = expectation_value(ψ1, H, envs) @test E₀ ≈ E1 atol = 1.0e-2 @@ -35,7 +38,7 @@ verbosity_conv = 1 Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H]) - @testset "Finite LazySum $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + @testset "Finite LazySum $(algname(alg))" for alg in algs ψ, envs = timestep(ψ₀, Hlazy, 0.0, dt, alg) E = expectation_value(ψ, Hlazy, envs) @test (3 + 1.55 - 0.1) * E₀ ≈ E atol = 1.0e-2 @@ -43,7 +46,7 @@ verbosity_conv = 1 Ht = MultipliedOperator(H, t -> 4) + MultipliedOperator(H, 1.45) - @testset "Finite TimeDependent LazySum $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + @testset "Finite TimeDependent LazySum $(algname(alg))" for alg in algs ψ, envs = timestep(ψ₀, Ht(1.0), 0.0, dt, alg) E = expectation_value(ψ, Ht(1.0), envs) @@ -54,7 +57,7 @@ verbosity_conv = 1 Ht2 = MultipliedOperator(H, t -> t < 0 ? error("t < 0!") : 4) + MultipliedOperator(H, 1.45) - @testset "Finite TimeDependent LazySum (fix negative t issue) $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + @testset "Finite TimeDependent LazySum (fix negative t issue) $(algname(alg))" for alg in algs ψ, envs = timestep(ψ₀, Ht2, 0.0, dt, alg) E = expectation_value(ψ, Ht2(0.0), envs) @@ -63,6 +66,65 @@ verbosity_conv = 1 @test E ≈ Et atol = 1.0e-8 end + # imaginary-time evolution lowers the energy toward the ground state. This is generic to + # symmetric integrators; BUG in particular has no backward substep (cf. its docstring), so it is + # a natural imaginary-time integrator. + @testset "Finite imaginary-time lowers energy $(algname(alg))" for alg in [TDVP(), BUG()] + Random.seed!(5) + ψi = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + E_start = real(expectation_value(ψi, H)) + E_prev = E_start + for _ in 1:20 + ψi, = timestep(ψi, H, 0.0, 0.1, alg; imaginary_evolution = true, normalize = true) + E_now = real(expectation_value(ψi, H)) + @test E_now ≤ E_prev + 1.0e-6 # monotone (non-increasing) energy + E_prev = E_now + end + @test E_prev < E_start - 1.0 # substantial lowering toward the ground state + @test norm(ψi) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + end + + # second-order (temporal) convergence on a small full-rank system, which isolates the temporal + # order. This is a BUG-only test: at full rank single-site TDVP / two-site TDVP2 integrate the + # *exact* dynamics (the tangent-space projector is the identity, cf. Lubich/Haegeman), so their + # overlap error sits at the numerical floor and a convergence slope is undefined. BUG's + # augment-and-Galerkin step has a genuine temporal error, so it shows a clean ≥ 2nd-order slope + # (often ≈ 4 on full-rank systems); we assert the order floor. + @testset "second-order convergence (BUG)" begin + Random.seed!(2) + Lc = 4 + Hc = force_planar(transverse_field_ising(ComplexF64, Trivial; L = Lc)) + ψ_full = FiniteMPS(rand, ComplexF64, Lc, ℙ^2, ℙ^4) # full-rank: 1,2,4,2,1 + + Hmat = convert(TensorMap, Hc) + ψvec = convert(TensorMap, ψ_full) + ψvec /= norm(ψvec) + + T = 0.5 + dts = [0.1, 0.05, 0.025] + errs = map(dts) do δt + n = round(Int, T / δt) + ref = exp(-im * Hmat * (n * δt)) * ψvec + ψ = copy(ψ_full) + envs = environments(ψ, Hc, ψ) + for k in 0:(n - 1) + timestep!(ψ, Hc, k * δt, δt, BUG(), envs) + end + ψout = convert(TensorMap, ψ) + ψout /= norm(ψout) + return 1 - abs(dot(ψout, ref)) + end + + slopes = [ + (log(errs[i + 1]) - log(errs[i])) / (log(dts[i + 1]) - log(dts[i])) + for i in 1:(length(dts) - 1) + ] + @info "BUG convergence" errs slopes + for s in slopes + @test s ≥ 1.7 + end + end + H = repeat(force_planar(heisenberg_XXX(; spin = 1)), 2) ψ₀ = InfiniteMPS([ℙ^3, ℙ^3], [ℙ^50, ℙ^50]) E₀ = expectation_value(ψ₀, H) @@ -93,6 +155,106 @@ verbosity_conv = 1 end end +# BUG is a drop-in for TDVP on the finite interface: over a few real-time steps of a generic +# (random) state the two integrators must agree, up to their shared temporal-discretization error. +# This is a cross-check between integrators, so it has no meaning as a TDVP-vs-TDVP test. +@testset "BUG agrees with TDVP" begin + L = 10 + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + Random.seed!(1234) + ψr = complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)) + δt = 0.01 + ψ_bug, ψ_tdvp = ψr, ψr + for k in 0:4 + ψ_bug, = timestep(ψ_bug, H, k * δt, δt, BUG()) + ψ_tdvp, = timestep(ψ_tdvp, H, k * δt, δt, TDVP()) + end + @test expectation_value(ψ_bug, H) ≈ expectation_value(ψ_tdvp, H) atol = 1.0e-3 + @test abs(dot(ψ_bug, ψ_tdvp)) ≈ 1 atol = 1.0e-3 +end + +# Genuine symmetric-tensor coverage (no `force_planar`) for the single-site finite integrators. +# Both TDVP and BUG must conserve the energy and accrue only an eigenstate phase, and must preserve +# the total boundary charge (the fixed `right` virtual space at site L). TDVP2 is excluded here: it +# requires a `trscheme` and is the two-site variant; these are single-site conservation properties. +@testset "Finite symmetric-tensor time evolution" verbose = true begin + dt = 0.1 + L = 6 + algs = [TDVP(), BUG()] + + # U(1)-symmetric Heisenberg, both in the natural total-Sz = 0 sector and in a fixed nonzero + # total-charge (Sz = 1) sector. + @testset "U(1) Heisenberg (total Sz = $label, $(algname(alg)))" for (label, right) in + (("0", U1Space(0 => 1)), ("1", U1Space(1 => 1))), alg in algs + Random.seed!(2718) + H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψ = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]; right) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, alg) + E1 = expectation_value(ψ1, H, envs) + + @test E₀ ≈ E1 atol = 1.0e-2 + @test imag(E1) ≈ 0 atol = 1.0e-8 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + @test right_virtualspace(ψ1, L) == right_virtualspace(ψ₀, L) + end + + # A second symmetry group. Z2 (transverse-field Ising) and SU2 (Heisenberg); same assertions. + @testset "Z2 transverse-field Ising ($(algname(alg)))" for alg in algs + Random.seed!(161803) + H = transverse_field_ising(ComplexF64, Z2Irrep; g = 1.0, L) + ψ = FiniteMPS(physicalspace(H), Z2Space(0 => 4, 1 => 4)) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, alg) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + @test right_virtualspace(ψ1, L) == right_virtualspace(ψ₀, L) + end + + @testset "SU(2) Heisenberg ($(algname(alg)))" for alg in algs + Random.seed!(577215) + H = heisenberg_XXX(ComplexF64, SU2Irrep; spin = 1 // 2, L) + # SU(2) spin-1/2 bonds alternate between integer / half-integer spins, so use the + # model's own full-rank virtual spaces rather than a hand-picked (integer-only) space. + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψ = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + ψ1, envs = timestep(ψ₀, H, 0.0, dt, alg) + E1 = expectation_value(ψ1, H, envs) + @test E₀ ≈ E1 atol = 1.0e-2 + @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1.0e-4 + @test right_virtualspace(ψ1, L) == right_virtualspace(ψ₀, L) + end + + # Imaginary-time symmetric evolution lowers the energy while preserving the sector + norm. + @testset "imaginary-time (U(1), $(algname(alg)))" for alg in algs + Random.seed!(141421) + H = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L) + maxV = MPSKit.max_virtualspaces(physicalspace(H)) + ψi = FiniteMPS(physicalspace(H), maxV[2:(end - 1)]) + Rtot = right_virtualspace(ψi, L) # total boundary charge, fixed throughout + E_start = real(expectation_value(ψi, H)) + E_prev = E_start + for _ in 1:15 + ψi, = timestep(ψi, H, 0.0, 0.1, alg; imaginary_evolution = true, normalize = true) + E_now = real(expectation_value(ψi, H)) + @test E_now ≤ E_prev + 1.0e-6 # monotone (non-increasing) energy + E_prev = E_now + end + @test E_prev < E_start - 0.5 # substantial lowering toward the ground state + @test norm(ψi) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + @test right_virtualspace(ψi, L) == Rtot # total boundary charge conserved throughout + end +end + @testset "Finite CBE-TDVP" verbose = true begin L = 10 H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) @@ -169,14 +331,14 @@ end @testset "time_evolve" verbose = true begin t_span = 0:0.1:0.1 - algs = [TDVP(), TDVP2(; trscheme = truncrank(10))] + algs = [TDVP(), TDVP2(; trscheme = truncrank(10)), BUG()] L = 10 H = force_planar(heisenberg_XXX(; spin = 1 // 2, L)) ψ₀ = FiniteMPS(L, ℙ^2, ℙ^1) E₀ = expectation_value(ψ₀, H) - @testset "Finite $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + @testset "Finite $(algname(alg))" for alg in algs ψ, envs = time_evolve(ψ₀, H, t_span, alg) E = expectation_value(ψ, H, envs) @test E₀ ≈ E atol = 1.0e-2 From f8a54a22271ea097d5e9d154f0e0d461fa93e49e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 20 Jul 2026 19:44:21 -0400 Subject: [PATCH 6/6] update implementation --- src/algorithms/timestep/bug.jl | 191 ++++++++++++++++++++------------- 1 file changed, 114 insertions(+), 77 deletions(-) diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl index 4a83a62b1..35d42d637 100644 --- a/src/algorithms/timestep/bug.jl +++ b/src/algorithms/timestep/bug.jl @@ -9,11 +9,20 @@ Unlike [`TDVP`](@ref), BUG advances both the basis (K-step) and the core (Galerk *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. -Passing a truncating `trscheme` (anything other than `notrunc()`) switches on **rank-adaptivity**: -each half-sweep augments every bond with the new directions discovered by the evolved connecting -tensor (old basis first, `[U₀ │ K₁]`) and then truncates back down to the tolerance of `trscheme`. -The bond dimension grows and shrinks automatically to track the entanglement; `notrunc()` recovers -the fixed-rank integrator. +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 @@ -29,18 +38,15 @@ $(TYPEDFIELDS) * [Ceruti et al. BIT Numer. Math. 62 (2022)](@cite ceruti2022) """ -struct BUG{A, O, T, S, F} <: Algorithm +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 - "truncation scheme used to cut the bond back down for rank-adaptive BUG" - trscheme::T - - "algorithm used for the singular value decomposition" - alg_svd::S + "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 @@ -50,7 +56,79 @@ function BUG(; trscheme = notrunc(), alg_svd = Defaults.alg_svd(), finalize = Defaults._finalize ) - return BUG(integrator, alg_orth, trscheme, alg_svd, 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!( @@ -61,77 +139,36 @@ function timestep!( # symmetric 2nd-order: a dt/2 left→right half-sweep composed with its dt/2 mirror L = length(ψ) h = dt / 2 - truncates = !(alg.trscheme isa MatrixAlgebraKit.NoTruncation) - svdcut = SvdCut(; trscheme = alg.trscheme, alg_svd = alg.alg_svd) - - # left→right half-sweep (root = last site) - ψ.AC[1] # gauge center to site 1 - ψ_old = copy(ψ) # frozen bases / reprojection inputs - T = isomorphism(scalartype(ψ), left_virtualspace(ψ_old, 1) ← left_virtualspace(ψ_old, 1)) - for i in 1:(L - 1) - Ĉ = _mul_front(T, ψ_old.AC[i]) # reproject old connecting tensor - AC = integrate(AC_hamiltonian(i, ψ, H, ψ, envs), Ĉ, t, h, alg.integrator; imaginary_evolution) - U₀ = _mul_front(T, ψ_old.AL[i]) # old left isometry in the new frame - if truncates - U = _bug_augment_left(U₀, AC, alg.alg_orth) # old-first augment; cut deferred to changebonds! - C = U' * AC - else - U, C = left_gauge(AC, alg.alg_orth) - end - T = U' * U₀ # transport (new ← old) - ψ.AC[i] = (U, C) + + # 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 - AC = integrate( - AC_hamiltonian(L, ψ, H, ψ, envs), _mul_front(T, ψ_old.AC[L]), - t, h, alg.integrator; imaginary_evolution - ) - normalize && normalize!(AC) - ψ.AC[L] = AC - truncates && changebonds!(ψ, svdcut; normalize) - - # right→left half-sweep (root = first site), the mirror - ψ.AC[L] # gauge center to site L - ψ_old = copy(ψ) - T = isomorphism(scalartype(ψ), right_virtualspace(ψ_old, L) ← right_virtualspace(ψ_old, L)) - for i in L:-1:2 - Ĉ = ψ_old.AC[i] * T # reproject old connecting tensor - AC = integrate(AC_hamiltonian(i, ψ, H, ψ, envs), Ĉ, t + h, h, alg.integrator; imaginary_evolution) - U₀ = ψ_old.AR[i] * T # old right isometry in the new frame - if truncates - U = _bug_augment_right(U₀, AC, alg.alg_orth) - C = _transpose_tail(AC) * _transpose_tail(U)' - else - C, U = right_gauge(AC, alg.alg_orth) - end - T = _transpose_tail(U₀) * _transpose_tail(U)' # transport (old ← new) - ψ.AC[i] = (C, U) + + # 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 - AC = integrate( - AC_hamiltonian(1, ψ, H, ψ, envs), ψ_old.AC[1] * T, - t + h, h, alg.integrator; imaginary_evolution - ) - normalize && normalize!(AC) - ψ.AC[1] = AC - truncates && changebonds!(ψ, svdcut; normalize) return ψ, envs end -# augment the RIGHT bond (left→right sweep): orthonormalize the stacked `[U₀ │ K₁]` (old isometry -# first) so `U₀` stays the leading per-sector block and only the new directions of `K₁` are appended. -function _bug_augment_left(U₀, K₁, alg_orth = Defaults.alg_orth()) - Û, _ = left_orth(catdomain(U₀, K₁); alg = alg_orth) - return Û -end - -# mirror of `_bug_augment_left` for the right→left sweep, on the `_transpose_tail` form (right-isometry -# with orthonormal rows): an LQ orthonormalizes the stacked rows `[U₀; K₁]`, keeping `U₀` leading. -function _bug_augment_right(U₀, K₁, alg_orth = Defaults.alg_orth()) - stacked = catcodomain(_transpose_tail(U₀), _transpose_tail(K₁)) - _, Û = right_orth(stacked; alg = alg_orth) - return _transpose_front(Û) -end - # copying version function timestep( ψ::AbstractFiniteMPS, H, time::Number, timestep::Number,