diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index aa0075e72..76f9fc9fe 100644 --- a/docs/src/assets/mpskit.bib +++ b/docs/src/assets/mpskit.bib @@ -67,6 +67,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 05a4af08d..7472c9cfb 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -28,9 +28,19 @@ When releasing a new version, move the "Unreleased" changes to a new version sec truncation. The sweep direction is selected by the `left_to_right` keyword. Both `approximate((O, ϕ), alg)` and `approximate!(ψ, (O, ϕ), alg)` are supported, where the destination `ψ` is a write target rather than an initial guess and may alias `ϕ`; they return `(ψ, ϵ)`. +- `BUG` time-evolution algorithm: a 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 `trunc` 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 29a9bf6b5..02b078b61 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 ba4a3aeea..920a8fefa 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 NoiseSchedule, FunctionalSchedule, ExponentialDecay, Warmup, DMRG3S @@ -166,6 +166,7 @@ include("algorithms/post_expand/post_expand.jl") include("algorithms/post_expand/dmrg3s.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..428b6b50e --- /dev/null +++ b/src/algorithms/timestep/bug.jl @@ -0,0 +1,205 @@ +""" +$(TYPEDEF) + +Single-site time-evolution algorithm for **finite** MPS, based on the +Basis-Update & Galerkin (BUG) integrator, an unconventional robust integrator for dynamical low-rank approximation with an observed second-order convergence. + +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 more natural choice for imaginary-time (dissipative) evolution, where the backward core step of the conventional projector-splitting integrator [`TDVP`](@ref) can become unstable for large timesteps. + +## Fields + +$(TYPEDFIELDS) + +## Algorithm + +Each half-sweep visits every site in turn and, at each site, +(i) splits off the bond ahead of it (in the sweep direction) with `alg_gauge`, truncating it back to `trunc`, +(ii) evolves the connecting tensor over `dt/2`, and +(iii) augments the basis with the new directions discovered by the evolved tensor (old basis first, `[U₀ │ K₁]`, orthonormalized with `alg_orth`) + +Notably, this last step does not include any truncation, and is meant to truncate the *previous* half-sweep's augmentation. +As a result, a truncation scheme `truncrank(D)` will result in a final MPS of dimension `2D`. +To restore a maximal dimension of `D`, apply [`changebonds`](@ref) with an [`SvdCut`](@ref) algorithm. + +!!! note + By default the state is not renormalized, as the (loss of) norm might accumulate useful information, + such as 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. + +!!! tip + Pass a `TimerOutputs.TimerOutput` as `timeroutput` to `timestep`/`timestep!` to obtain a + breakdown of the time spent in the three steps above (`cut_bond`, `AC_integrate`, `augment`) + over all local updates. + +## 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 orthonormalize the augmented basis `[U₀ │ K₁]` after each local update" + alg_orth::O + + "factorization used to gauge and truncate the bond ahead of each local update" + 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(), + trunc = notrunc(), alg_svd = Defaults.alg_svd(), + finalize = Defaults._finalize + ) + alg_gauge = _build_inner_gauge(trunc, alg_svd, alg_orth) + return BUG(integrator, alg_orth, alg_gauge, finalize) +end + +# `ψ.AC[site]` first, the neighbour second: the lazy `CView` walk keys off what is already cached, and +# both reads must precede the installs, which move the gauge center across the bond. +function _cut_bond!(site::Int, ::Val{:right}, ψ, alg) + AC₀ = ψ.AC[site] + AR_next = ψ.AR[site + 1] + + # nothing changes without truncation, so `ψ` is left alone — which is also why the non-mutating + # `left_orth` is required here: `AC₀` is still the center that gets evolved next + if !_truncates(alg) + AL_old, C₀ = left_orth(AC₀; alg) + return AL_old, C₀, AR_next, zero(real(scalartype(ψ))) + end + + # `Vᴴ` goes into the neighbour rather than into `C₀`, so the local update happens in the truncated + # bond space; in-place, as the installs below invalidate the `AC` cache anyway + AL_old, C₀, Vᴴ, ϵ = svd_trunc!(AC₀, alg) + AR_next = _mul_front(Vᴴ, AR_next) # old right block, rotated into the cut bond + ψ.AC[site] = (AL_old, C₀) + ψ.AC[site + 1] = (C₀, AR_next) + return AL_old, C₀, AR_next, ϵ +end + +# mirror on the `_transpose_tail` form, which the augment step also works in +function _cut_bond!(site::Int, ::Val{:left}, ψ, alg) + AC₀ = ψ.AC[site] + AL_prev = ψ.AL[site - 1] + + if !_truncates(alg) + C₀, AR_old_tail = right_orth(_transpose_tail(AC₀); alg) + return AR_old_tail, C₀, AL_prev, zero(real(scalartype(ψ))) + end + + U, C₀, AR_old_tail, ϵ = svd_trunc!(_transpose_tail(AC₀), alg) + AL_prev = _mul_tail(AL_prev, U) # old left block, rotated into the cut bond + ψ.AC[site] = (C₀, _transpose_front(AR_old_tail)) + ψ.AC[site - 1] = (AL_prev, C₀) + return AR_old_tail, C₀, AL_prev, ϵ +end + +# the site at which a sweep runs out of chain: no bond ahead to cut, so the local update is a plain +# evolve-and-install +_sweep_end(ψ, ::Val{:right}) = length(ψ) +_sweep_end(ψ, ::Val{:left}) = 1 + +# augment the old canonical basis at `site` with the new directions discovered by the evolved center +# (old first, and without truncation), and install it — together with the old center `C₀` transported +# into that basis and absorbed into the neighbouring site, which becomes the new gauge center +function _augment_basis!(site::Int, dir::Val{:right}, ψ, AL_old, AC, C₀, AR_next, alg_orth) + AL_new, _ = left_orth!(catdomain(AL_old, AC); alg = alg_orth) + transport = AL_new' * AL_old # old bond, in the augmented basis + return set_canonical!(ψ, site, dir, AL_new, _mul_front(transport * C₀, AR_next)) +end +function _augment_basis!(site::Int, dir::Val{:left}, ψ, AR_old_tail, AC, C₀, AL_prev, alg_orth) + _, AR_new_tail = right_orth!( + catcodomain(AR_old_tail, _transpose_tail(AC)); alg = alg_orth + ) + transport = AR_old_tail * AR_new_tail' + return set_canonical!( + ψ, site, dir, _transpose_front(AR_new_tail), _mul_tail(AL_prev, C₀ * transport) + ) +end + +function _evolve_center(site, ψ, H, alg::BUG, envs, t, h; imaginary_evolution) + Heff = AC_hamiltonian(site, ψ, H, ψ, envs) + return integrate(Heff, ψ.AC[site], t, h, alg.integrator; imaginary_evolution) +end + +function local_update!( + site, direction::Val, ψ, H, alg::BUG, envs, t, h; + imaginary_evolution, normalize, timeroutput + ) + # at the far end of the sweep there is no bond ahead to cut: evolve and finalize + if site == _sweep_end(ψ, direction) + AC = @timeit timeroutput "AC_integrate" _evolve_center( + site, ψ, H, alg, envs, t, h; imaginary_evolution + ) + normalize && normalize!(AC) + ψ.AC[site] = AC + return ψ + end + + # 1. cut (and truncate) the bond ahead, before evolving + A_old, C₀, neighbour, _ = @timeit timeroutput "cut_bond" _cut_bond!( + site, direction, ψ, alg.alg_gauge + ) + + # 2. evolve the connecting tensor + AC = @timeit timeroutput "AC_integrate" _evolve_center( + site, ψ, H, alg, envs, t, h; imaginary_evolution + ) + + # 3. augment the basis (old first, no truncation here) and install it, together with the + # transported old center at the neighbouring site + return @timeit timeroutput "augment" _augment_basis!( + site, direction, ψ, A_old, AC, C₀, neighbour, alg.alg_orth + ) +end + +function timestep!( + ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::BUG, + envs::AbstractMPSEnvironments = environments(ψ, H, ψ); + imaginary_evolution::Bool = false, normalize::Bool = false, + timeroutput::TimerOutput = DISABLED_TIMER + ) + L = length(ψ) + h = dt / 2 + + # left→right half-sweep (root = last site): `t → t + dt / 2` + @timeit timeroutput "half-sweep" for site in 1:L + ψ = local_update!( + site, Val(:right), ψ, H, alg, envs, t, h; + imaginary_evolution, normalize, timeroutput + ) + end + + # right→left half-sweep (root = first site): `t + dt / 2 → t + dt` + @timeit timeroutput "half-sweep" for site in L:-1:1 + ψ = local_update!( + site, Val(:left), ψ, H, alg, envs, t + h, h; + imaginary_evolution, normalize, timeroutput + ) + 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 && scalartype(H) <: Real)) + ψ′ = 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 06f2ef603..9a4fbf533 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 not 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 @@ -44,7 +45,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(; @@ -65,11 +66,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) @@ -124,23 +128,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) @@ -158,7 +162,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) @@ -167,9 +171,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) @@ -222,14 +226,14 @@ Used as the `algorithm` argument of [`timestep`](@ref), [`timestep!`](@ref) and "algorithm used for truncation of the two-site update" trunc::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 @@ -239,6 +243,7 @@ function timestep!( ac2′ = integrate(Hac2, ac2, t, dt / 2, alg.integrator; imaginary_evolution) nal, nc, nar = svd_trunc!(ac2′; trunc = alg.trunc, alg = alg.alg_svd) + normalize && normalize!(nc) ψ.AC[i] = (nal, complex(nc)) ψ.AC[i + 1] = (complex(nc), _transpose_front(nar)) @@ -258,6 +263,7 @@ function timestep!( ac2′ = integrate(Hac2, ac2, t + dt / 2, dt / 2, alg.integrator; imaginary_evolution) nal, nc, nar = svd_trunc!(ac2′; trunc = alg.trunc, alg = alg.alg_svd) + normalize && normalize!(nc) ψ.AC[i - 1] = (nal, complex(nc)) ψ.AC[i] = (complex(nc), _transpose_front(nar)) @@ -277,7 +283,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(ψ) @@ -290,5 +296,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 5410011d3..85193791d 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -17,13 +17,14 @@ through each of the time points obtained by iterating t_span. - `verbosity::Int = 0`: verbosity level for logging - `imaginary_evolution::Bool = false`: if true, the time evolution is done with an imaginary time step - 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. + instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful 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, which can be useful + to retain numerical stability when the norm loss is not information that is needed. # Returns -- `ψ`: the time-evolved state +- `ψ`: the time-stepped state - `envs`: the updated environment manager """ function time_evolve end, function time_evolve! end @@ -32,7 +33,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 @@ -41,7 +42,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) @@ -71,9 +72,10 @@ solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``. # Keyword Arguments - `imaginary_evolution::Bool = false`: if true, the time evolution is done with an imaginary time step - 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. + instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful 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, which can be useful + to retain numerical stability when the norm loss is not information that is needed. # Returns @@ -111,8 +113,7 @@ Construct an `MPO` that approximates ``\\exp(-iHdt)``. # Keyword Arguments -- `imaginary_evolution::Bool = false`: if true, the time evolution operator is constructed - with an imaginary time step 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. +- `imaginary_evolution::Bool = false`: if true, the time evolution is done with an imaginary time step + instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful to compute the + ground state of a Hamiltonian, or to compute finite-temperature properties of a system. """ make_time_mpo diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index b0d3d467a..4bbddf2a1 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -182,6 +182,78 @@ function Base.setindex!( end end +@doc """ + set_AL_AC!(ψ, site, AL, AC) -> ψ + set_AC_AR!(ψ, site, AC, AR) -> ψ + +Install a canonical tensor at `site` together with the neighbouring center tensor, in a single +update: `set_AL_AC!` writes the left-isometric `AL` at `site` and the center tensor `AC` at +`site + 1`, `set_AC_AR!` writes the right-isometric `AR` at `site` and the center tensor `AC` at +`site - 1`. The gauge center therefore ends up at `site + 1` and `site - 1` respectively. + +These are the pendants of `ψ.AC[site] = (AL, C)` / `ψ.AC[site] = (C, AR)` for algorithms that +already know the *next* center tensor: they avoid materializing a bond tensor purely to keep `ψ` +well-defined, and — unlike installing the canonical tensor on its own — they never leave `ψ` +without a gauge center. + +The isometric nature of `AL`/`AR` is not verified. +""" +set_AL_AC! +@doc (@doc set_AL_AC!) set_AC_AR! + +function set_AL_AC!(ψ::FiniteMPS, site::Int, AL::GenericMPSTensor, AC::GenericMPSTensor) + site < length(ψ) || + throw(ArgumentError(lazy"cannot set the center beyond the end of the chain at site $site")) + # materialize the halves that survive this update + site > 1 && ψ.AL[site - 1] + site + 1 < length(ψ) && ψ.AR[site + 2] + + ψ.ACs .= missing + ψ.Cs .= missing + ψ.ALs[(site + 1):end] .= missing + ψ.ARs[1:(site + 1)] .= missing + + ψ.ALs[site] = AL + ψ.ACs[site + 1] = AC + return ψ +end +function set_AC_AR!(ψ::FiniteMPS, site::Int, AC::GenericMPSTensor, AR::GenericMPSTensor) + site > 1 || + throw(ArgumentError(lazy"cannot set the center before the start of the chain at site $site")) + site < length(ψ) && ψ.AR[site + 1] + site - 1 > 1 && ψ.AL[site - 2] + + ψ.ACs .= missing + ψ.Cs .= missing + ψ.ALs[(site - 1):end] .= missing + ψ.ARs[1:(site - 1)] .= missing + + ψ.ARs[site] = AR + ψ.ACs[site - 1] = AC + return ψ +end + +set_AL_AC!(ψ::WindowMPS, args...) = (set_AL_AC!(ψ.window, args...); ψ) +set_AC_AR!(ψ::WindowMPS, args...) = (set_AC_AR!(ψ.window, args...); ψ) + +""" + set_canonical!(ψ, site, direction, A, AC) -> ψ + +Direction-dispatching wrapper around [`set_AL_AC!`](@ref)/[`set_AC_AR!`](@ref): install the canonical +tensor `A` at `site` — left-isometric for `direction = Val(:right)`, right-isometric for `Val(:left)` +— together with the center tensor `AC` at the next site *in that direction*, `site + 1` or `site - 1`. + +This is the pendant of [`gauge!`](@ref) for a sweep that already knows its next center tensor, and +lets a sweep body be written once for both directions. +""" +function set_canonical!( + ψ::AbstractFiniteMPS, site::Int, ::Val{Dir}, A::GenericMPSTensor, AC::GenericMPSTensor + ) where {Dir} + Dir === :right && return set_AL_AC!(ψ, site, A, AC) + Dir === :left && return set_AC_AR!(ψ, site, AC, A) + return throw(ArgumentError(lazy"invalid direction `$Dir`")) +end + function Base.getindex(v::ACView{<:WindowMPS, E}, i::Int)::E where {E} (i >= 1 && i <= length(v.parent)) || throw(ArgumentError("out of bounds")) return ACView(v.parent.window)[i] diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index 81696fe8a..70a49ab2d 100644 --- a/test/algorithms/timestep.jl +++ b/test/algorithms/timestep.jl @@ -9,15 +9,22 @@ 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))) + +maxbond(ψ) = maximum(i -> dim(left_virtualspace(ψ, i)), 1:length(ψ)) + @testset "timestep" verbose = true begin dt = 0.1 - algs = [TDVP(), TDVP2(; trunc = truncrank(10))] + # every rank-adaptive algorithm gets the same cap: without one, BUG's augmentation doubles the + # bond dimension on every half-sweep with nothing ever cutting it back + algs = [TDVP(), TDVP2(; trunc = truncrank(10)), BUG(; trunc = truncrank(10))] L = 10 H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) @@ -26,7 +33,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 +42,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 +50,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 +61,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) @@ -93,6 +100,86 @@ verbosity_conv = 1 end end +# BUG-specific: unlike `TDVP` there is no backward-in-time substep, which makes it a natural +# imaginary-time integrator (cf. its docstring). `TDVP`'s imaginary-time behaviour is covered by the +# CBE block below. +@testset "Finite imaginary-time BUG" begin + L = 10 + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + + alg = BUG(; trunc = truncrank(4)) + + 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:8 + ψ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 + @test maxbond(ψi) ≤ 8 # the cap bounds the bond at `2D` throughout + 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 # `normalize = true` renormalizes each step +end + +# The truncating path, which is BUG's distinguishing feature. `BUG(; trunc = …)` cuts the bond +# *ahead of* every local update — the bond carrying the previous half-sweep's augmentation — and +# augments without truncating, so the bond dimension oscillates between the requested rank and (at +# most) twice it. +@testset "BUG with truncation" begin + L = 8 + H = force_planar(transverse_field_ising(ComplexF64, Trivial; L, g = 1.5)) + δt, nsteps, D = 0.05, 20, 4 + + Random.seed!(11) + ψ₀ = FiniteMPS(rand, ComplexF64, L, ℙ^2, ℙ^2) # low rank to start with, entanglement grows + normalize!(ψ₀) + + tovec(ψ) = (v = convert(TensorMap, ψ); v / norm(v)) + ref = exp(-im * convert(TensorMap, H) * (nsteps * δt)) * tovec(ψ₀) + + # the cap actually bites here: the exact state needs bond 16 + ψ = ψ₀ + for k in 0:(nsteps - 1) + ψ, = timestep(ψ, H, k * δt, δt, BUG(; trunc = truncrank(D))) + @test maxbond(ψ) ≤ 2D # `D` right after the cut, at most `2D` after the augment + end + @test 1 - abs(dot(tovec(ψ), ref)) < 1.0e-3 +end + +# Genuine symmetric-tensor coverage (no `force_planar`) for the single-site finite integrators. +# BUG's augment step stacks bases per sector, so it can add or drop sectors: both `TDVP` and `BUG` +# must conserve the energy, accrue only an eigenstate phase, and preserve the total boundary charge +# (the fixed `right` virtual space at site L). TDVP2 is excluded here: it requires a `trunc` and is +# the two-site variant; these are single-site conservation properties. +@testset "Finite symmetric-tensor time evolution" begin + dt = 0.1 + L = 6 + + 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)]) + ψ₀, = find_groundstate(ψ, H; verbosity = 0) + E₀ = expectation_value(ψ₀, H) + + # the exact state has Schmidt rank ≤ 8 here, so this cap bounds the bond without discarding + # anything — it only stops BUG's augmentation from doubling unchecked + algs = [TDVP(), BUG(; trunc = truncrank(8))] + + @testset "U(1) Heisenberg ($(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 + @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 +end + @testset "Finite CBE-TDVP" verbose = true begin L = 10 H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) @@ -122,8 +209,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 +220,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 +247,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 +256,14 @@ end @testset "time_evolve" verbose = true begin t_span = 0:0.1:0.1 - algs = [TDVP(), TDVP2(; trunc = truncrank(10))] + algs = [TDVP(), TDVP2(; trunc = truncrank(10)), BUG(; trunc = truncrank(10))] 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