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..7cbcba8cc 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -21,8 +21,19 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Added +- `BUG` time-evolution algorithm: a symmetric second-order Basis-Update & Galerkin integrator for + finite MPS. Unlike `TDVP` it has no backward-in-time substep (stable for imaginary-time evolution), + and passing a truncating `trscheme` enables rank-adaptivity (the bond dimension grows and shrinks + automatically to track entanglement). + ### Changed +- Renormalization during time evolution is now controlled by an explicit `normalize` keyword on + `timestep`/`time_evolve` (default `false`), decoupled from `imaginary_evolution`. By default the + norm is preserved, so it retains useful information (the accumulated truncation error in real time, + or the decaying weight in imaginary time). Previously imaginary-time evolution always renormalized + every step; **to recover that behavior, pass `normalize = true`** (e.g. for ground-state or + thermal-state search via imaginary-time evolution). - `environments` now follows a single positional contract for every state and operator kind: `environments(below, operator, above, alg)`, where `alg` is the environment algorithm (slot 4). The operator form requires an explicit `above`. Auxiliary inputs are keyword-only: 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..35d42d637 --- /dev/null +++ b/src/algorithms/timestep/bug.jl @@ -0,0 +1,190 @@ +""" +$(TYPEDEF) + +Single-site, symmetric second-order time-evolution algorithm for **finite** MPS, based on the +Basis-Update & Galerkin (BUG) integrator, an unconventional robust integrator for dynamical low-rank +approximation. + +Unlike [`TDVP`](@ref), BUG advances both the basis (K-step) and the core (Galerkin C-step) tensors +*forward* in time, with no backward-in-time substep. This makes it a natural choice for +imaginary-time / dissipative evolution, where the backward core step of TDVP can become unstable. + +Each half-sweep augments every bond with the new directions discovered by the evolved connecting +tensor (old basis first, `[U₀ │ K₁]`) and truncates back down to the tolerance of `trscheme` in the +same orthonormalization step. The truncation is folded into the augment orth (a truncated SVD of the +stacked basis), so the old directions — whose orthonormal columns carry singular value `≥ 1` — are +always kept and only newly-appended directions are cut. The bond therefore *grows* to track the +entanglement under a `truncerror` tolerance; a hard rank cap (`truncrank`) can additionally shrink it. + +!!! warning + With `trscheme = notrunc()` the augmented basis is kept at full rank every half-sweep, so the + bond grows unboundedly (up to the local Hilbert-space dimension) — this is **not** a fixed-rank + integrator. To cap the bond dimension at `D`, pass `trscheme = truncrank(D)`: because the old + directions carry singular value `≥ 1`, a rank-`D` cut keeps them and drops the newly appended + directions, so a bond already at `D` stays there (a fixed-rank-`D` step). A bond currently below + `D` can still grow up to `D` as new directions are admitted. + +!!! note + By default the state is not renormalized, so the norm keeps useful information (the + accumulated truncation error in real time, or the decaying weight in imaginary time). + Pass `normalize = true` to `timestep`/`time_evolve` to renormalize after every half-sweep + instead. This is independent of `imaginary_evolution`. + +## Fields + +$(TYPEDFIELDS) + +## References + +* [Ceruti et al. BIT Numer. Math. 62 (2022)](@cite ceruti2022) +""" +struct BUG{A, O, G, F} <: Algorithm + "algorithm used in the exponential solvers" + integrator::A + + "algorithm used to re-orthonormalize the basis after each local update" + alg_orth::O + + "factorization used for the in-sweep augment gauge: a truncated SVD (`alg_svd` with `trscheme`)" + alg_gauge::G + + "callback function applied after each iteration, of signature `finalize(t, ψ, H, envs) -> ψ, envs`" + finalize::F +end +function BUG(; + integrator = Defaults.alg_expsolve(), alg_orth = Defaults.alg_orth(), + trscheme = notrunc(), alg_svd = Defaults.alg_svd(), + finalize = Defaults._finalize + ) + alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme) + return BUG(integrator, alg_orth, alg_gauge, finalize) +end + +# left→right BUG update at `site`: evolve the connecting tensor, old-first augment + truncate the +# basis, and return the `transport` and center `AC_old` advanced to `site + 1` (root = last site). +function local_update!( + site, ::Val{:right}, ψ, H, alg::BUG, envs, t, h, transport, AC_old, ARs; + imaginary_evolution, normalize + ) + # 1. Transport the old AC to the new frame and evolve + AC₀ = _mul_front(transport, AC_old) + AC = integrate( + AC_hamiltonian(site, ψ, H, ψ, envs), AC₀, + t, h, alg.integrator; imaginary_evolution + ) + + # If end of chain, finalize + if site == length(ψ) + normalize && normalize!(AC) + ψ.AC[site] = AC + return ψ, transport, AC_old + end + + # 2. move gauge right + oldbasis, C₀ = left_gauge(AC₀, alg.alg_orth) # old AL - C in new frame + AC_next = _mul_front(C₀, ARs[site + 1]) + + + # 3. reproject onto new basis + newbasis, _ = left_gauge(catdomain(oldbasis, AC), alg.alg_gauge) # augmented basis + newbond = newbasis' * AC + ψ.AC[site] = (newbasis, newbond) + new_transport = newbasis' * oldbasis + + return ψ, new_transport, AC_next +end + +# mirror of the left→right update on the `_transpose_tail` form; `transport` and `AC_old` are returned +# advanced to `site - 1` (root = first site). +function local_update!( + site, ::Val{:left}, ψ, H, alg::BUG, envs, t, h, transport, AC_old, ALs; + imaginary_evolution, normalize + ) + # 1. Transport the old AC to the new frame and evolve + AC₀ = AC_old * transport + AC = integrate( + AC_hamiltonian(site, ψ, H, ψ, envs), AC₀, + t, h, alg.integrator; imaginary_evolution + ) + + # If start of chain, finalize + if site == 1 + normalize && normalize!(AC) + ψ.AC[site] = AC + return ψ, transport, AC_old + end + + # 2. move gauge left + C₀, oldbasis = right_gauge(AC₀, alg.alg_orth) # old C - AR in new frame + AC_next = _mul_tail(ALs[site - 1], C₀) + + # 3. reproject onto new basis + _, newbasis_tail = right_orth( + catcodomain(_transpose_tail(oldbasis), _transpose_tail(AC)); + alg = MatrixAlgebraKit.RightOrthViaSVD(alg.alg_gauge) # augmented basis + ) + newbasis = _transpose_front(newbasis_tail) + newbond = _transpose_tail(AC) * _transpose_tail(newbasis)' + ψ.AC[site] = (newbond, newbasis) + new_transport = _transpose_tail(oldbasis) * _transpose_tail(newbasis)' + + return ψ, new_transport, AC_next +end + +function timestep!( + ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::BUG, + envs::AbstractMPSEnvironments = environments(ψ, H, ψ); + imaginary_evolution::Bool = false, normalize::Bool = false + ) + # symmetric 2nd-order: a dt/2 left→right half-sweep composed with its dt/2 mirror + L = length(ψ) + h = dt / 2 + + # left→right half-sweep (root = last site): freeze the bases as `[AC[1], AR[2], …, AR[L]]` and + # carry `AC_old`/`transport` forward + ARs = ψ.AR[2:end] + pushfirst!(ARs, ψ.AC[1]) + transport = isomorphism(scalartype(ψ), left_virtualspace(ψ, 1) ← left_virtualspace(ψ, 1)) + AC_old = ARs[1] + for site in 1:L + ψ, transport, AC_old = local_update!( + site, Val(:right), ψ, H, alg, envs, t, h, transport, AC_old, ARs; + imaginary_evolution, normalize + ) + end + + # right→left half-sweep (root = first site), the mirror: freeze `[AL[1], …, AL[L-1], AC[L]]` and + # carry `AC_old`/`transport` backward (starting from `t + h`) + ALs = ψ.AL[1:(L - 1)] + push!(ALs, ψ.AC[L]) + transport = isomorphism(scalartype(ψ), right_virtualspace(ψ, L) ← right_virtualspace(ψ, L)) + AC_old = ALs[L] + for site in L:-1:1 + ψ, transport, AC_old = local_update!( + site, Val(:left), ψ, H, alg, envs, t + h, h, transport, AC_old, ALs; + imaginary_evolution, normalize + ) + end + + return ψ, envs +end + +# copying version +function timestep( + ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, + alg::BUG, envs::AbstractMPSEnvironments...; + imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... + ) + isreal = (scalartype(ψ) <: Real && !imaginary_evolution) + ψ′ = isreal ? complex(ψ) : copy(ψ) + if length(envs) != 0 && isreal + @warn "Currently cannot reuse real environments for complex evolution" + envs′ = environments(ψ′, H, ψ′) + elseif length(envs) == 1 + envs′ = only(envs) + else + @assert length(envs) == 0 "Invalid signature" + envs′ = environments(ψ′, H, ψ′) + end + return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) +end diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 993782e54..44b39a54b 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 @@ -40,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(; @@ -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) @@ -214,14 +218,14 @@ $(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 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/bug_augment.jl b/test/algorithms/bug_augment.jl new file mode 100644 index 000000000..bc1996f8f --- /dev/null +++ b/test/algorithms/bug_augment.jl @@ -0,0 +1,144 @@ +println(" +------------------------------------ +| BUG basis-augmentation tests | +------------------------------------ +") + +using .TestSetup +using Test, TestExtras +using MPSKit +using MPSKit: _transpose_tail, _transpose_front, left_orth, right_orth, right_gauge, Defaults +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)): +# 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 diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index a846b2c5c..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)) @@ -122,8 +284,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 +295,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 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(ψit) ≈ 1 atol = 1.0e-6 # imaginary-time renormalizes each step + @test norm(ψn) ≈ 1 atol = 1.0e-6 end @testset "imaginary-time lowers energy" begin @@ -150,7 +322,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 @@ -159,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