From cff38f46c0719e8b8b9210c2c30c86285ca5b25f Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Thu, 13 Aug 2026 11:56:29 +0200 Subject: [PATCH 01/29] track truncation error in one-site tdvp --- src/algorithms/timestep/tdvp.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index a3da18c52..352211c89 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -129,7 +129,9 @@ function _timestep_infinite( end recalculate!(envs, ψ′, H) - return ψ′, envs + # infinite one-site TDVP has a fixed bond dimension and never truncates so nothing is discarded + # the gauge-fixing residual is controlled by `tolgauge`, not reported here + return ψ′, envs, zero(real(scalartype(ψ′))) end function timestep!( @@ -148,6 +150,10 @@ function _timestep_finite!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs, allocator; imaginary_evolution::Bool, normalize::Bool ) + # discarded weight of the step, accumulated in squares over the local gauges + # stays exactly zero for QR gauging + ϵ² = zero(real(scalartype(ψ))) + # sweep left to right for i in 1:(length(ψ) - 1) # 1. optionally expand the bond ahead of the local update (CBE) @@ -161,7 +167,8 @@ function _timestep_finite!( # 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. By default the norm is # preserved; `normalize` renormalizes. - left_gauge!(ψ, i, AC, alg.alg_gauge; normalize) + _, ϵ = left_gauge!(ψ, i, AC, alg.alg_gauge; normalize) + ϵ² += ϵ^2 # 4. evolve the bond tensor backward Hc = C_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator) @@ -190,7 +197,8 @@ function _timestep_finite!( # 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) + _, ϵ = right_gauge!(ψ, i, AC, alg.alg_gauge; normalize) + ϵ² += ϵ^2 # 4. evolve the bond tensor backward Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator) @@ -207,7 +215,7 @@ function _timestep_finite!( imaginary_evolution ) - return ψ, envs + return ψ, envs, sqrt(ϵ²) end """ From 7f2c4bf0cfa0167d1c5b79fa1d1c8a5f18b24e02 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Thu, 13 Aug 2026 12:37:09 +0200 Subject: [PATCH 02/29] two-site tdvp error via gauge2! --- src/algorithms/timestep/tdvp.jl | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 352211c89..17ea533ed 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -274,16 +274,21 @@ function _timestep2_finite!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, envs, allocator; imaginary_evolution::Bool, normalize::Bool ) + # the two-site center always has to be split back up, so the gauge is always a truncated SVD + alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trunc) + + # discarded weight of the step, accumulated in squares over the local gauges + ϵ² = zero(real(scalartype(ψ))) + # sweep left to right for i in 1:(length(ψ) - 1) ac2 = _transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1]) Hac2 = AC2_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator) 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)) + # the discarded weight is the truncation error + _, ϵ = gauge2!(ψ, i, Val(:right), ac2′, alg_gauge; normalize) + ϵ² += ϵ^2 if i != (length(ψ) - 1) Hac = AC_hamiltonian(i + 1, ψ, H, ψ, envs; alg.backend, allocator) @@ -300,10 +305,8 @@ function _timestep2_finite!( Hac2 = AC2_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator) 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)) + _, ϵ = gauge2!(ψ, i - 1, Val(:left), ac2′, alg_gauge; normalize) + ϵ² += ϵ^2 if i != 2 Hac = AC_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator) @@ -314,7 +317,7 @@ function _timestep2_finite!( end end - return ψ, envs + return ψ, envs, sqrt(ϵ²) end # copying version From 37713fee9ce4cc1207a35d0c43fc394dbef253ec Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Thu, 13 Aug 2026 12:55:04 +0200 Subject: [PATCH 03/29] return errors in timestep and time_evolve --- src/algorithms/timestep/time_evolve.jl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index b4a1f8043..3f421c367 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -35,21 +35,25 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time envs = environments(ψ, H, ψ); verbosity::Int = 0, imaginary_evolution::Bool = false, normalize::Bool = false ) - log = IterLog("TDVP") + log = IterLog(string(nameof(typeof(alg)))) + # discarded weight over the whole evolution, accumulated in squares over the steps + ϵ² = zero(real(scalartype(ψ))) LoggingExtras.withlevel(; verbosity) do - @infov 2 loginit!(log, 0, first(t_span)) + @infov 2 loginit!(log, 0.0, first(t_span)) for iter in 1:(length(t_span) - 1) t = t_span[iter] dt = t_span[iter + 1] - t - ψ, envs = $timestep(ψ, H, t, dt, alg, envs; imaginary_evolution, normalize) + ψ, envs, ϵ = $timestep(ψ, H, t, dt, alg, envs; imaginary_evolution, normalize) ψ, envs = alg.finalize(t, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)} + ϵ² += ϵ^2 - @infov 3 logiter!(log, iter, 0, t) + # per-step truncation error; the running total is what is reported at the end + @infov 3 logiter!(log, iter, convert(Float64, ϵ), t) end - @infov 2 logfinish!(log, length(t_span), 0, t_span[end]) + @infov 2 logfinish!(log, length(t_span), convert(Float64, sqrt(ϵ²)), t_span[end]) end - return ψ, envs + return ψ, envs, sqrt(ϵ²) end end From 3b4fed6c981d811acff5e227057a049e3530e63d Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Thu, 13 Aug 2026 16:51:07 +0200 Subject: [PATCH 04/29] update docstrings + be more precise about what the returned error means in time evolution --- src/algorithms/timestep/tdvp.jl | 18 +++++++++----- src/algorithms/timestep/time_evolve.jl | 33 ++++++++++++++++++++++---- 2 files changed, 41 insertions(+), 10 deletions(-) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 17ea533ed..532528517 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -11,11 +11,16 @@ the enlarged bond back down (selecting the truncated-SVD gauge). The expansion i state-preserving, as required for a consistent time evolution. !!! note - 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. + By default the norm is not preserved: neither the bond expansion nor the truncation renormalizes, + so the state norm keeps useful information. In real time this is exact, namely the squared norm + drops by precisely the truncated ("discarded") weight, + ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon^2``, + with `ϵ` the truncation error returned by [`timestep`](@ref). In imaginary time the norm also carries + the physical decay of the weight, so it no longer isolates the truncation. Without `trunc` nothing is + discarded at all and the norm is conserved exactly in real 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 @@ -222,6 +227,7 @@ end $(TYPEDEF) Two-site MPS time-evolution algorithm based on the Time-Dependent Variational Principle. +See [`TDVP`](@ref) for more information. # Fields @@ -286,7 +292,7 @@ function _timestep2_finite!( Hac2 = AC2_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator) ac2′ = integrate(Hac2, ac2, t, dt / 2, alg.integrator; imaginary_evolution) - # the discarded weight is the truncation error + # the norm of the discarded singular values is the truncation error _, ϵ = gauge2!(ψ, i, Val(:right), ac2′, alg_gauge; normalize) ϵ² += ϵ^2 diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 3f421c367..e520d7b9d 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -1,6 +1,6 @@ """ - time_evolve(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ, envs) - time_evolve!(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ₀, envs) + time_evolve(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ, envs, ϵ) + time_evolve!(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ₀, envs, ϵ) Time-evolve the initial state `ψ₀` with Hamiltonian `H` over a given time span by stepping through each of the time points obtained by iterating t_span. @@ -26,6 +26,13 @@ through each of the time points obtained by iterating t_span. - `ψ`: the time-stepped state - `envs`: the updated environment manager +- `ϵ`: the truncation error accumulated over the whole evolution, i.e. the per-step errors of + [`timestep`](@ref), ``\\epsilon = \\sqrt{\\sum_{\\text{steps}} \\epsilon_i^2}``. + In real time with `normalize = false` this is exactly the norm lost to truncation over the whole + evolution, ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon^2``. This is zero for + algorithms that never truncate. See [`timestep`](@ref) for what it does and does not measure. + +The per-step error is logged at `verbosity ≥ 3` and the accumulated total at `verbosity ≥ 2`. """ function time_evolve end, function time_evolve! end @@ -58,8 +65,8 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time end """ - timestep(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ, envs) - timestep!(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ₀, envs) + timestep(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ, envs, ϵ) + timestep!(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ₀, envs, ϵ) Time-step the state `ψ₀` with Hamiltonian `H` over a given time step `dt` at time `t`, solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``. @@ -85,6 +92,24 @@ solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``. - `ψ`: the time-stepped state - `envs`: the updated environment manager +- `ϵ`: the truncation error of the step (see below) + +# Truncation error + +`ϵ` is the truncation error of the step, i.e. the 2-norm of the sum of squared singular values discarded +by the local gauge factorizations. In other words, `ϵ²` is the truncated ("discarded") weight. + +It is nonzero only for algorithms that truncate such as [`TDVP2`](@ref), [`BUG`](@ref) with a `trunc`, and +[`TDVP`](@ref) with a bond expansion, while it is exactly `0` for one-site [`TDVP`](@ref), +which runs at fixed bond dimension. A zero `ϵ` does not mean the step was exact, but that this +particular error channel is absent. In particular, the projection and time-discretisation errors are never +included in `ϵ`. + +In real time with `normalize = false`, `ϵ` is exactly the norm lost to truncation, +``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon^2``. + +See [Errors and accuracy](@ref) in the manual for what the other error sources are, why the +per-bond errors combine in a squared manner, and the precise statement and caveats of the norm identity. # Examples From 154be132a9f4b00901816a84b22d411be6d0f12a Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Thu, 13 Aug 2026 17:02:07 +0200 Subject: [PATCH 05/29] do the same for BUG --- src/algorithms/timestep/bug.jl | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl index 1ea9edd4d..4f53916c4 100644 --- a/src/algorithms/timestep/bug.jl +++ b/src/algorithms/timestep/bug.jl @@ -23,10 +23,18 @@ As a result, a truncation scheme `truncrank(D)` will result in a final MPS of di 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. + By default the state is not renormalized, as the (loss of) norm accumulates useful information. + In real time the squared norm lost is exactly the weight discarded by the bond cuts, + ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon^2``, with `ϵ` the truncation error + returned by [`timestep`](@ref). In imaginary time the norm also carries the physical decay of + the weight and no longer isolates the truncation. `ϵ` is exactly zero when not truncating. Pass `normalize = true` to `timestep`/`time_evolve` to renormalize after every half-sweep instead. +!!! warning + `ϵ` counts only the cuts in step (i). Step (iii) augments the basis without truncating, so within + a sweep the bond dimension temporarily exceeds `trunc` and the norm identity above applies to the + state returned at the end of the step. + !!! 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`) @@ -139,11 +147,11 @@ function local_update!( ) normalize && normalize!(AC) ψ.AC[site] = AC - return ψ + return ψ, zero(real(scalartype(ψ))) end # 1. cut (and truncate) the bond ahead, before evolving - A_old, C₀, neighbour, _ = @timeit timeroutput "cut_bond" _cut_bond!( + A_old, C₀, neighbour, ϵ = @timeit timeroutput "cut_bond" _cut_bond!( site, direction, ψ, alg.alg_gauge ) @@ -154,9 +162,10 @@ function local_update!( # 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!( + ψ = @timeit timeroutput "augment" _augment_basis!( site, direction, ψ, A_old, AC, C₀, neighbour, alg.alg_orth ) + return ψ, ϵ end function timestep!( @@ -171,23 +180,28 @@ function timestep!( # the sweep is serial, so a single allocator serves all local updates allocator = default_allocator(ψ, SerialScheduler()) + # discarded weight of the step, accumulated in squares over the bond cuts + ϵ² = zero(real(scalartype(ψ))) + # left→right half-sweep (root = last site): `t → t + dt / 2` @timeit timeroutput "half-sweep" for site in 1:L - ψ = local_update!( + ψ, ϵ = local_update!( site, Val(:right), ψ, H, alg, envs, t, h, allocator; imaginary_evolution, normalize, timeroutput ) + ϵ² += ϵ^2 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!( + ψ, ϵ = local_update!( site, Val(:left), ψ, H, alg, envs, t + h, h, allocator; imaginary_evolution, normalize, timeroutput ) + ϵ² += ϵ^2 end - return ψ, envs + return ψ, envs, sqrt(ϵ²) end # copying version From 08e1f16b8072eeaa4ccf253c86f6ea299643e1ef Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 14 Aug 2026 11:44:56 +0200 Subject: [PATCH 06/29] be elaborate on what the errors mean per algorithm --- src/algorithms/approximate/approximate.jl | 17 ++++++++++++++++- src/algorithms/approximate/zipup.jl | 2 +- src/algorithms/groundstate/dmrg.jl | 11 +++++++---- src/algorithms/groundstate/find_groundstate.jl | 8 +++++++- .../groundstate/gradient_grassmann.jl | 4 +++- src/algorithms/groundstate/idmrg.jl | 6 ++++-- src/algorithms/groundstate/vumps.jl | 3 ++- src/algorithms/propagator/corvector.jl | 4 +++- src/algorithms/statmech/leading_boundary.jl | 6 ++++-- src/algorithms/statmech/vomps.jl | 3 ++- src/algorithms/timestep/wii.jl | 4 ++-- src/states/ortho.jl | 4 ++-- src/states/orthoview.jl | 4 +++- 13 files changed, 56 insertions(+), 20 deletions(-) diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index 7cfb5dc12..f48942ecd 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -29,12 +29,27 @@ algorithm for you based on the type of `ψ₀` (`DMRG`/`DMRG2` for a finite MPS, `IDMRG2` for an infinite MPS) and only accepts the `(O, ψ)` tuple form of `toapprox`. Once you pass an explicit `algorithm`, keywords are no longer accepted here — configure the algorithm struct itself instead (e.g. `DMRG(; tol, maxiter, verbosity)`). -- `tol::Float64`: tolerance for convergence criterium +- `tol::Float64`: convergence tolerance, compared against the Galerkin error (see Returns below) - `maxiter::Int`: maximum amount of iterations - `verbosity::Int`: display progress information - `trunc`: if supplied, a truncated two-site sweep (`DMRG2`/`IDMRG2`) is prepended to refine the bond dimension before the single-site algorithm polishes the result. +# Returns + +- `ψ`: the approximated state +- `environments`: environments corresponding to the result (not returned by `Zipup`, which uses none) +- `ϵ::Float64`: an error measure whose meaning depends on the algorithm: + - for the iterative algorithms (`DMRG`, `DMRG2`, `IDMRG`, `IDMRG2`, `VOMPS`) it is the final + convergence error, i.e. the Galerkin error compared against `tol`. This is the same quantity + [`find_groundstate`](@ref) returns, measuring distance from the variational fixed point. + - for [`Zipup`](@ref) it is instead a truncation error, i.e. the largest 2-norm of the discarded + singular values over all bonds and sweeps. `Zipup` is a single non-iterative sweep, so there is + no convergence measure to report and no `tol` to compare against. + + The two are not comparable, and a small `ϵ` means different things in each case. See the manual on + the `ϵ` convention under [The error convention](@ref). + # Algorithms Each algorithm below only supports a subset of the general interface. Check this table before diff --git a/src/algorithms/approximate/zipup.jl b/src/algorithms/approximate/zipup.jl index 06b6ec401..c178c9db5 100644 --- a/src/algorithms/approximate/zipup.jl +++ b/src/algorithms/approximate/zipup.jl @@ -97,7 +97,7 @@ opposite direction imposes a final truncation with `alg_zipdown` in a locally ga the gauge center of `ψ` at the far end. The destination may alias `ϕ`. Also returns the truncation error `ϵ`, the largest 2-norm of the discarded singular values over all -bonds and both sweeps. +bonds and both sweeps, i.e. the worst single bond. """ zip_left_right! @doc (@doc zip_left_right!) zip_right_left! diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 9c996dbcb..940e8e503 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -60,7 +60,8 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`approximate`](@ref). """ struct DMRG{A, F, E, G, B} <: Algorithm - "tolerance for convergence criterium" + "convergence tolerance, compared against the Galerkin error (the tangent-space gradient + norm)" tol::Float64 "maximal amount of iterations" @@ -171,7 +172,8 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`approximate`](@ref). """ struct DMRG2{A, G, F, B} <: Algorithm - "tolerance for convergence criterium" + "convergence tolerance, compared against the Galerkin error (the tangent-space gradient + norm)" tol::Float64 "maximal amount of iterations" @@ -232,7 +234,7 @@ function local_update!( alg_gauge = _update_alg_gauge(alg.alg_gauge, iter, ϵ_global) # 2. gauge: truncated SVD split back into single-site tensors and install; - # the discarded weight is the truncation error + # the norm of the discarded singular values is the truncation error ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(ψ, pos, direction, O, envs, newA2center, alg_gauge; normalize = true) # 3. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1) @@ -273,7 +275,8 @@ Currently supported for the finite-system algorithms [`DMRG`](@ref) and [`DMRG2` - `ψ::AbstractFiniteMPS`: converged ground state - `environments`: environments corresponding to the converged state -- `ϵ::Float64`: final convergence error upon terminating the algorithm +- `ϵ::Float64`: final convergence error upon terminating the algorithm, + i.e. the Galerkin error (see [`find_groundstate`](@ref)) """ function find_groundstate!( ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs = environments(ψ, H, ψ) diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index 29d37271f..e4bd2af0c 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/src/algorithms/groundstate/find_groundstate.jl @@ -38,7 +38,13 @@ low-bond-dimension initial guess such as a product state. - `ψ::AbstractMPS`: converged ground state - `environments`: environments corresponding to the converged state -- `ϵ::Float64`: final convergence error upon terminating the algorithm +- `ϵ::Float64`: final convergence error upon terminating the algorithm, i.e. the quantity compared + against the algorithm's `tol`. It measures distance from a variational fixed point and is not + a truncation error. See the manual under [The error convention](@ref)) for more information. + Which measure it is depends on the algorithm: the sweeping algorithms + ([`DMRG`](@ref), [`DMRG2`](@ref), [`VUMPS`](@ref), [`IDMRG`](@ref), [`IDMRG2`](@ref)) + report the Galerkin error, whereas [`GradientGrassmann`](@ref) reports + the norm of the Riemannian gradient from its optimizer. # Examples diff --git a/src/algorithms/groundstate/gradient_grassmann.jl b/src/algorithms/groundstate/gradient_grassmann.jl index 0159a0ef7..ca323a5c7 100644 --- a/src/algorithms/groundstate/gradient_grassmann.jl +++ b/src/algorithms/groundstate/gradient_grassmann.jl @@ -14,7 +14,9 @@ with a preconditioner to induce the metric from the Hilbert space inner product. - `method = ConjugateGradient`: instance of optimization algorithm, or type of optimization algorithm to construct - `finalize!`: finalizer algorithm -- `tol = Defaults.tol`: tolerance for convergence criterium +- `tol = Defaults.tol`: convergence tolerance, compared against the norm of the Riemannian + (Grassmann) gradient reported by the optimizer. This is also the `ϵ` returned by + [`find_groundstate`](@ref) for this algorithm. - `maxiter = Defaults.maxiter`: maximum amount of iterations - `verbosity = Defaults.verbosity - 1`: level of information display - `hasconverged = OptimKit.DefaultHasConverged(tol)`: convergence criterium diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index e57d11de2..b8732417d 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -12,7 +12,8 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_boundary`](@ref), and [`approximate`](@ref). """ @kwdef struct IDMRG{A, B} <: Algorithm - "tolerance for convergence criterium" + "convergence tolerance, compared against the Galerkin error (the tangent-space gradient + norm)" tol::Float64 = Defaults.tol "maximal amount of iterations" @@ -45,7 +46,8 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_boundary`](@ref), and [`approximate`](@ref). """ @kwdef struct IDMRG2{A, S, B} <: Algorithm - "tolerance for convergence criterium" + "convergence tolerance, compared against the Galerkin error (the tangent-space gradient + norm)" tol::Float64 = Defaults.tol "maximal amount of iterations" diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index 716d2e21d..a7c6d5419 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -17,7 +17,8 @@ Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`leading_bou * [Vanderstraeten et al. SciPost Phys. Lect. Notes 7 (2019)](@cite vanderstraeten2019) """ @kwdef struct VUMPS{F, B} <: Algorithm - "tolerance for convergence criterium" + "convergence tolerance, compared against the Galerkin error (the tangent-space gradient + norm)" tol::Float64 = Defaults.tol "maximal amount of iterations" diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index 6b4f3066b..2f2fb3118 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -28,7 +28,9 @@ Used as the `algorithm` argument of [`propagator`](@ref). flavour::F = NaiveInvert() "algorithm used for the linear solvers" solver::S = Defaults.linearsolver - "tolerance for convergence criterium" + "convergence tolerance, compared against the largest change in a center tensor over a sweep, + `maxᵢ ‖ACᵢ′ - ACᵢ‖`. This represents a measure of how much the sweep still moves the state, + not a residual of the linear system (that is controlled by `solver`)" tol::Float64 = Defaults.tol * 10 "maximal amount of iterations" maxiter::Int = Defaults.maxiter diff --git a/src/algorithms/statmech/leading_boundary.jl b/src/algorithms/statmech/leading_boundary.jl index 218a3f8f7..de1694cdf 100644 --- a/src/algorithms/statmech/leading_boundary.jl +++ b/src/algorithms/statmech/leading_boundary.jl @@ -14,7 +14,7 @@ optimization algorithm will be attempted based on the supplied keywords. # Keyword Arguments -- `tol::Float64`: tolerance for convergence criterium +- `tol::Float64`: convergence tolerance, compared against the Galerkin error (see Returns below) - `maxiter::Int`: maximum amount of iterations - `verbosity::Int`: display progress information @@ -22,7 +22,9 @@ optimization algorithm will be attempted based on the supplied keywords. - `ψ::AbstractMPS`: converged leading boundary MPS - `environments`: environments corresponding to the converged boundary -- `ϵ::Float64`: final convergence error upon terminating the algorithm +- `ϵ::Float64`: final convergence error upon terminating the algorithm, i.e. the Galerkin error. + It is not a truncation error; see [`find_groundstate`](@ref) and the manual on the + `ϵ` convention under [The error convention](@ref), and [Ground-state accuracy](@ref). """ leading_boundary # TODO: alg selector diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl index 9427f5488..dfbd81d1b 100644 --- a/src/algorithms/statmech/vomps.jl +++ b/src/algorithms/statmech/vomps.jl @@ -18,7 +18,8 @@ Used as the `algorithm` argument of [`leading_boundary`](@ref) and [`approximate * [Vanhecke et al. SciPost Phys. Core 4 (2021)](@cite vanhecke2021) """ @kwdef struct VOMPS{F, B} <: Algorithm - "tolerance for convergence criterium" + "convergence tolerance, compared against the Galerkin error (the tangent-space gradient + norm)" tol::Float64 = Defaults.tol "maximal amount of iterations" diff --git a/src/algorithms/timestep/wii.jl b/src/algorithms/timestep/wii.jl index e41405fd2..101da56d3 100644 --- a/src/algorithms/timestep/wii.jl +++ b/src/algorithms/timestep/wii.jl @@ -17,9 +17,9 @@ Used as the `algorithm` argument of [`make_time_mpo`](@ref). * [Paeckel et al. Ann. of Phys. 411 (2019)](@cite paeckel2019) """ @kwdef struct WII <: Algorithm - "tolerance for convergence criterium" + "tolerance of the Arnoldi exponentiation used to build each local block" tol::Float64 = Defaults.tol - "maximal number of iterations" + "maximal number of iterations of that exponentiation" maxiter::Int = Defaults.maxiter end diff --git a/src/states/ortho.jl b/src/states/ortho.jl index 9dfc803e4..ed03e41bc 100644 --- a/src/states/ortho.jl +++ b/src/states/ortho.jl @@ -17,7 +17,7 @@ $(TYPEDFIELDS) Used as the `alg` argument of [`gaugefix!`](@ref). """ @kwdef struct LeftCanonical <: Algorithm - "tolerance for convergence criterium" + "convergence tolerance, compared against the residual of the gauge fixed-point iteration" tol::Float64 = Defaults.tolgauge "maximal amount of iterations" maxiter::Int = Defaults.maxiter @@ -46,7 +46,7 @@ $(TYPEDFIELDS) Used as the `alg` argument of [`gaugefix!`](@ref). """ @kwdef struct RightCanonical <: Algorithm - "tolerance for convergence criterium" + "convergence tolerance, compared against the residual of the gauge fixed-point iteration" tol::Float64 = Defaults.tolgauge "maximal amount of iterations" maxiter::Int = Defaults.maxiter diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index 4bbddf2a1..0e990d889 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -313,7 +313,9 @@ standard MPS-tensor form. Passing a [`TruncatedAlgorithm`](@extref MatrixAlgebraKit.TruncatedAlgorithm) instead performs a truncated SVD may shrink the bond. Also returns the truncation error `ϵ`: the 2-norm of the discarded singular values from the -truncated SVD, or `0` for a norm-preserving QR/LQ gauge. +truncated SVD, or `0` for a norm-preserving QR/LQ gauge. The package-wide convention has `ϵ` +represent an amplitude, so that `ϵ²` is the truncated ("discarded") weight and the squared norm of the +factorized tensor drops by exactly `ϵ²`. """ left_gauge @doc (@doc left_gauge) right_gauge From bc0c7e239b397aa50b467df17279116ba11c2d0c Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 14 Aug 2026 11:49:03 +0200 Subject: [PATCH 07/29] documentation on errors and accuracy for various algorithms --- docs/src/assets/mpskit.bib | 55 +++++++++ docs/src/man/algorithms.md | 153 ++++++++++++++++++++++++- src/algorithms/timestep/time_evolve.jl | 2 +- 3 files changed, 206 insertions(+), 4 deletions(-) diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index 44e762933..28d9eb56e 100644 --- a/docs/src/assets/mpskit.bib +++ b/docs/src/assets/mpskit.bib @@ -999,3 +999,58 @@ @article{hubig2015 doi = {10.1103/PhysRevB.91.155115}, url = {https://link.aps.org/doi/10.1103/PhysRevB.91.155115} } + +@article{hubig2018, + title = {Error estimates for extrapolations with matrix-product states}, + author = {Hubig, C. and Haegeman, J. and Schollw\"ock, U.}, + journal = {Phys. Rev. B}, + volume = {97}, + issue = {4}, + pages = {045125}, + year = {2018}, + month = {Jan}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevB.97.045125}, + url = {https://link.aps.org/doi/10.1103/PhysRevB.97.045125} +} + +@article{lubich2015, + title = {Time Integration of Tensor Trains}, + author = {Lubich, Christian and Oseledets, Ivan V. and Vandereycken, Bart}, + journal = {SIAM Journal on Numerical Analysis}, + volume = {53}, + number = {2}, + pages = {917--941}, + year = {2015}, + publisher = {Society for Industrial and Applied Mathematics}, + doi = {10.1137/140976546}, + url = {https://doi.org/10.1137/140976546} +} + +@article{li2024, + title = {Time-{{Dependent Variational Principle}} with {{Controlled Bond Expansion}} for {{Matrix Product States}}}, + author = {Li, Jheng-Wei and Gleis, Andreas and {von Delft}, Jan}, + year = {2024}, + month = jul, + journal = {Physical Review Letters}, + volume = {133}, + number = {2}, + pages = {026401}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevLett.133.026401}, + url = {https://link.aps.org/doi/10.1103/PhysRevLett.133.026401} +} + +@article{schollwoeck2011, + title = {The density-matrix renormalization group in the age of matrix product states}, + author = {Schollw{\"o}ck, Ulrich}, + year = {2011}, + month = jan, + journal = {Annals of Physics}, + volume = {326}, + number = {1}, + pages = {96--192}, + issn = {0003-4916}, + doi = {10.1016/j.aop.2010.09.012}, + url = {https://www.sciencedirect.com/science/article/pii/S0003491610001752} +} diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 02b078b61..60a76ae2e 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -16,6 +16,8 @@ In MPSKit.jl, this can be achieved through `find_groundstate`: find_groundstate ``` +The returned error measures convergence to a variational fixed point, which is not the same as accuracy; see [Ground-state accuracy](@ref). + There are a variety of algorithms that have been developed over the years, and many of them have been implemented in MPSKit. Keep in mind that some of them are exclusive to finite or infinite systems, while others may work for both. Many of these algorithms have different advantages and disadvantages, and figuring out the optimal algorithm is not always straightforward, since this may strongly depend on the model. @@ -49,7 +51,7 @@ DMRG2 For infinite systems, a similar approach can be used by dynamically adding new sites to the middle of the system and optimizing over them. This gradually increases the system size until the boundary effects are no longer felt. However, because of this approach, for critical systems this algorithm can be quite slow to converge, since the number of steps needs to be larger than the correlation length of the system. -Again, both a single-site and a two-site version are implemented, to have the option to dynamically increase the bonddimension at a higher cost. +Again, both a single-site and a two-site version are implemented, to have the option to dynamically increase the bond dimension at a higher cost. ```@docs; canonical=false IDMRG @@ -101,6 +103,11 @@ The first is focused around approximately solving the equation for a small times This can be achieved by projecting the equation onto the tangent space of the MPS, and then solving the results. This procedure is commonly referred to as the [`TDVP`](@ref) algorithm, which again has a two-site variant to allow for dynamically altering the bond dimension. +There are three ways to let the bond dimension follow the entanglement rather than fixing it up front: +- [`TDVP2`](@ref) evolves two sites at a time and splits the result back apart with a truncated SVD. +- [`TDVP`](@ref) with an `alg_expand` keeps the cheaper single-site update and instead expands the bond with directions orthogonal to the current state before each local update, recovering controlled bond expansion (CBE). +- [`BUG`](@ref) is a different integrator altogether: it advances basis and core tensors forward in time with no backward substep, which makes it better behaved for imaginary-time evolution, and it is rank-adaptive when given a `trunc`. + ```@docs; canonical=false TDVP TDVP2 @@ -119,10 +126,15 @@ WII TaylorCluster ``` +Time evolution has three distinct error sources, only one of which is reported back to the user. +See [Time evolution accuracy](@ref). + ## Excitations It might also be desirable to obtain information beyond the lowest energy state of a given system, and study the dispersion relation. While it is typically not feasible to resolve states in the middle of the energy spectrum, there are several ways to target a few of the lowest-lying energy states. +None of these report an error. +For what limits their accuracy, see [Excitation accuracy](@ref). ```@docs; canonical=false excitations @@ -263,6 +275,138 @@ Es, ϕs = excitations(H, ChepigaAnsatz2(), ψ, envs; num=1) isapprox(Es[1] - E₀, 2(g - 1); rtol=1e-2) # infinite analytical result ``` +## Errors and accuracy + +Most algorithms in MPSKit report an error alongside their result, and the manual pages above refer to it as `ϵ` throughout. +That single name covers genuinely different quantities, and the differences matter. +This section clarifies the differences, and explains what they do and, just as important, what they don't measure. + +### The error convention + +Where `ϵ` is a *truncation* error, it is the 2-norm of the discarded singular values of a single factorisation, so that `ϵ²` is the discarded weight and the squared norm of the factorised tensor drops by exactly `ϵ²`. +The name "discarded weight" refers to it truly representing a probability: the ``\sigma_\alpha^2`` are the eigenvalues of the reduced density matrix across the cut, i.e. the statistical weights of the Schmidt states, summing to 1 for a normalised state. +So `ϵ²` is the probability weight thrown away and `ϵ` is the corresponding amplitude, which is why squares appear wherever these errors are combined. +Note that `ϵ` is absolute rather than relative: it is not divided by the norm of the state, which under a non-renormalising algorithm drifts away from 1 precisely as truncation accumulates. + +What differs between algorithms is how the per-bond values are *aggregated*. +For this reason, values between algorithms are not directly comparable. + +!!! warning + Not every `ϵ` is a truncation error. + The `ϵ` returned by [`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms is a *convergence* measure, with no truncation interpretation at all. + The two are unrelated quantities that happen to share a name; see the two sections below. + +### Ground-state accuracy + +[`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms return the quantity their `tol` is compared against. +For the sweeping algorithms ([`DMRG`](@ref), [`DMRG2`](@ref), [`VUMPS`](@ref), [`IDMRG`](@ref), [`IDMRG2`](@ref)) this is the Galerkin error: the norm of the local gradient projected orthogonally to the current state. +It vanishes exactly at a variational fixed point. +[`GradientGrassmann`](@ref) instead reports the norm of the Riemannian gradient from its optimizer. +Both vanish at a fixed point and both are gradient norms. +However, since they are taken in different metrics, their magnitudes are not directly comparable. +In particular, a `tol` tuned for one is not a `tol` tuned for the other. + +In other words, convergence is only defined relative to the manifold you are optimising over. +A single-site algorithm at a fixed bond dimension can drive its `ϵ` to machine precision and still be far from the true ground state, because the error that remains is the bond dimension itself, which no amount of further sweeping can address. +A small `ϵ` certifies a fixed point, not an accurate state. +Growing the bond dimension is the job of the two-site algorithms ([`DMRG2`](@ref), [`IDMRG2`](@ref)) or of a bond expansion ([`DMRG`](@ref) with an `alg_expand`, or an expanding `alg_gauge` such as [`DMRG3S`](@ref)); see also [`changebonds`](@ref). + +Once an algorithm does truncate, the two error notions interact. +The Galerkin error cannot fall below the level set by the weight being discarded each sweep, so a truncating scheme converges once `ϵ` reaches the truncation error rather than the (unreachable) bare `tol`. +[`DMRG`](@ref)/[`DMRG2`](@ref) account for this: their stopping test is `ϵ ≤ max(tol, maximum(ϵ_trunc))`, which reduces to the plain `ϵ ≤ tol` when nothing is truncated. + +Neither measure is an error bar on an observable. +For that, the standard route is the energy variance ``\langle H^2 \rangle - \langle H \rangle^2`` and extrapolation of observables towards zero variance ([Hubig et al.](@cite hubig2018)). + +### Time evolution accuracy + +Unlike a ground-state search, a time evolution has no convergence criterion to run to. +There is no fixed point, and the error is made at every step. +Three sources behave differently and only one of them is reported. + +* **Truncation error.** + Whenever a bond is cut back down, the discarded singular values are lost from the state. + [`timestep`](@ref) and [`time_evolve`](@ref) return this as their third value `ϵ`, the norm of the discarded component, so that `ϵ²` is the discarded weight. + It is the error you control through the algorithm's `trunc`, and the only one that is free to compute, since the truncating SVD produces it anyway. + It is non-zero for [`TDVP2`](@ref), for [`BUG`](@ref) with a `trunc`, and for [`TDVP`](@ref) with a bond expansion. + Plain single-site [`TDVP`](@ref) runs at fixed bond dimension and returns exactly `0`. + +* **Projection error.** + Single-site [`TDVP`](@ref) confines the evolution to the tangent space of a fixed-bond-dimension manifold, ``\lVert (1 - P_{T_\psi}) H \psi \rVert``. + The component of the exact evolution pointing off that manifold is simply dropped, and this happens even with no truncation and exact local solves. + It is not reported as measuring it costs an extra effective-Hamiltonian application per site. + This is what a bond expansion (CBE) exists to reduce ([Li et al.](@cite li2024)). + +* **Time-discretization error.** + The projector splitting is globally ``O(dt^2)`` for the symmetric back-and-forth sweep ([Lubich et al.](@cite lubich2015), [Paeckel et al.](@cite paeckel2019)), so it is controlled by `dt` alone. + This can only be estimated by comparing one step of `dt` against two of `dt / 2`. + +A trustworthy run needs all three under control, not just a small `ϵ`. +In practice: pick `dt` from a convergence check, pick `trunc` from the accumulated `ϵ`, and use a bond-adaptive scheme ([`TDVP2`](@ref), [`BUG`](@ref), or [`TDVP`](@ref) with `alg_expand`) whenever entanglement grows during the evolution, since a fixed bond dimension silently converts entanglement growth into projection error. + +#### Summing local errors in squares and relation to the norm + +``\epsilon^2`` is the quantity that adds exactly. +Each local truncation is an orthogonal projection, so it removes exactly ``\epsilon_k^2`` from the squared norm, and the substeps between truncations preserve the norm. +Summing the squares therefore tracks a conserved "cost". +The alternative is to sum ``\epsilon_k``: ``\lVert \psi_{\text{untruncated}} - \psi \rVert \le \sum_k \epsilon_k``, a distance to the untruncated solution, which is a different and always larger quantity. + +By default none of the time evolution algorithms renormalize (`normalize = false`), which is deliberate. +In real time the local exponentials are unitary, so truncation is the only thing that changes the norm and it becomes a running record of what truncation has cost, + +```math +\lVert \psi \rVert^2 = \lVert \psi_0 \rVert^2 - \epsilon^2 . +``` + +The reported `ϵ` is the norm deficit, and this composes across steps. +This follows from the following two facts put together, one per half of a local update. +1) An SVD truncation is an orthogonal projection onto the kept Schmidt vectors and is 2-norm optimal at that rank ([Schollwöck](@cite schollwoeck2011)), so the kept and discarded parts are orthogonal. +By Pythagoras the squared norm drops by exactly the discarded weight, the usual way of quantifying truncation during a time evolution ([Paeckel et al.](@cite paeckel2019)). +2) The local exponentials of the projector-splitting sweep are unitary, so TDVP conserves the norm and the energy exactly when the local equations are solved exactly ([Paeckel et al.](@cite paeckel2019)), contributing nothing to the norm change. + +These two hold for the local updates of a step, so composing over all steps gives the identity. + +!!! note + "Exactly" in the second fact is up to the tolerance of the local exponentials, and is thus in practice only approximate due to integrator tolerance. + +It is also specific to real time with `normalize = false`: + +* **Imaginary time** evolves with the non-unitary ``\exp(-H dt)``, which rescales the state on its own. + The norm then moves for two independent reasons, namely the physical decay of the weight and the truncation loss. One cannot separate them from each other. `ϵ` still counts only the truncation. +* **`normalize = true`** renormalizes at every local update, destroying the identity by construction. + This is usually what you want for imaginary-time evolution used as a ground-state or thermal-state search. + `ϵ` is still reported and is unaffected. +* **No truncation at all** (plain single-site [`TDVP`](@ref), or [`BUG`](@ref) with a QR gauge) gives ``\epsilon = 0``, and in real time the norm is then conserved exactly. +* An **`InfiniteMPS`** is regauged to norm 1 per site structurally, so its norm carries no such information and `normalize` has no effect. + +### Excitation accuracy + +[`excitations`](@ref) returns only `(energies, states)`: there is no error term, and none of the sources below is reported back to you. +They are worth knowing about, because the dominant one is usually not the one the algorithm is working on. + +* **Inherited ground-state error.** + Every method here builds on a ground state you supply and treats it as exact. + Its error propagates straight into the gap, and since a gap is a difference of two large energies, it is typically the limiting factor. + A well-converged ground state (in the sense of `ϵ` and bond dimension) is necessary for a meaningful excitation calculation. + +* **Ansatz limitation.** + [`QuasiparticleAnsatz`](@ref) varies over the single-quasiparticle tangent space on top of a fixed ground state. + It is variational within that space and well suited to isolated quasiparticle branches, but multi-particle continua are not representable in it, so results there are not to be trusted. + For infinite systems the momentum superposition itself is exact, so momentum is a good quantum number and no error enters through it. + +* **Eigensolver convergence.** + The local eigenvalue problem is solved with KrylovKit, and a run that fails to converge `num` states emits a warning on the residual when the verbosity is set high enough. + This residual is neither returned nor thrown, so it is worth not running with warnings suppressed. + +* **Penalty-based orthogonality** ([`FiniteExcited`](@ref)). + Higher states are found by minimising ``H + \lambda \sum_i |\psi_i\rangle\langle\psi_i|`` against the previously converged states, with ``\lambda`` the `weight` field. + A finite `weight` enforces orthogonality only approximately, so a residual overlap with a lower state biases the reported energy downwards. + Since the reported value is the expectation value of the bare `H`, this bias is invisible in the output. + Raising `weight` suppresses it at the cost of stretching the spectrum and slowing down the eigensolver's per-gap eigensolves. + +* **Truncation** ([`ChepigaAnsatz2`](@ref)). + The two-site excited state is split back to single-site tensors with a truncated SVD governed by `trunc`, and the resulting discarded weight is not reported. + ## `changebonds` Many of the previously mentioned algorithms do not possess a way to dynamically change to @@ -275,15 +419,18 @@ state. changebonds ``` +All of these are controlled by a `trunc`, and the weight they discard is measured the same way as explained in the `ϵ` convention under [The error convention](@ref). +`changebonds` does not report it, since every algorithm has its own interpretation of the discarded singular values. + There are several different algorithms implemented, each having their own advantages and disadvantages: -* [`SvdCut`](@ref): The simplest method for changing the bonddimension is found by simply +* [`SvdCut`](@ref): The simplest method for changing the bond dimension is found by simply locally truncating the state using an SVD decomposition. This yields a (locally) optimal truncation, but clearly cannot be used to increase the bond dimension. Note that a globally optimal truncation can be obtained by using the [`SvdCut`](@ref) algorithm in combination with [`approximate`](@ref). Since the output of this method might have a - truncated bonddimension, the new state might not be identical to the input state. + truncated bond dimension, the new state might not be identical to the input state. The truncation is controlled through `trunc`, which dictates how the singular values of the original state are truncated. diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index e520d7b9d..0b692b185 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -108,7 +108,7 @@ included in `ϵ`. In real time with `normalize = false`, `ϵ` is exactly the norm lost to truncation, ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon^2``. -See [Errors and accuracy](@ref) in the manual for what the other error sources are, why the +See [Time evolution accuracy](@ref) in the manual for what the other error sources are, why the per-bond errors combine in a squared manner, and the precise statement and caveats of the norm identity. # Examples From 54d99d3e28dde5e17d7ff900711621e0282b7597 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 14 Aug 2026 11:53:56 +0200 Subject: [PATCH 08/29] typo --- src/algorithms/groundstate/find_groundstate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index e4bd2af0c..bc76686ab 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/src/algorithms/groundstate/find_groundstate.jl @@ -40,7 +40,7 @@ low-bond-dimension initial guess such as a product state. - `environments`: environments corresponding to the converged state - `ϵ::Float64`: final convergence error upon terminating the algorithm, i.e. the quantity compared against the algorithm's `tol`. It measures distance from a variational fixed point and is not - a truncation error. See the manual under [The error convention](@ref)) for more information. + a truncation error. See the manual under [The error convention](@ref) for more information. Which measure it is depends on the algorithm: the sweeping algorithms ([`DMRG`](@ref), [`DMRG2`](@ref), [`VUMPS`](@ref), [`IDMRG`](@ref), [`IDMRG2`](@ref)) report the Galerkin error, whereas [`GradientGrassmann`](@ref) reports From ef5d8493db34fee352017d0f7ccdb3c2c0369931 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 14 Aug 2026 12:01:58 +0200 Subject: [PATCH 09/29] add tests on time evolution errors --- test/algorithms/timestep.jl | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index c503dac45..39a3cd763 100644 --- a/test/algorithms/timestep.jl +++ b/test/algorithms/timestep.jl @@ -263,6 +263,36 @@ end end end +@testset "Truncation error" verbose = true begin + L = 10 + dt = 0.1 + H = force_planar(heisenberg_XXX(Float64, Trivial; spin = 1 // 2, L)) + Random.seed!(7) + ψ₀ = normalize!(complex(FiniteMPS(rand, Float64, L, ℙ^2, ℙ^16))) + + # fixed bond dimension sweep never discards anything, so the reported error is exactly zero + @testset "no truncation" begin + @test last(timestep(ψ₀, H, 0.0, dt, TDVP())) == 0 + @test last(timestep(ψ₀, H, 0.0, dt, BUG())) == 0 + end + + @testset "$(nameof(alg))" for alg in (TDVP2, BUG) + _, _, ϵ_loose = timestep(ψ₀, H, 0.0, dt, alg(; trunc = truncrank(32))) + ψ, _, ϵ_tight = timestep(ψ₀, H, 0.0, dt, alg(; trunc = truncrank(2))) + + @test ϵ_tight > 1.0e-3 # throw away real weight + @test ϵ_loose < ϵ_tight # throw away less weight with a more forgiving truncation + @test norm(ψ) < norm(ψ₀) # discarded weight = norm loss in real time with `normalize = false` + end + + @testset "accumulation over an evolution" begin + alg = TDVP2(; trunc = truncrank(2)) + _, _, ϵ_step = timestep(ψ₀, H, 0.0, dt, alg) + _, _, ϵ_total = time_evolve(ψ₀, H, 0:dt:(4 * dt), alg) + @test ϵ_total >= ϵ_step + end +end + @testset "time_evolve" verbose = true begin t_span = 0:0.1:0.1 algs = [TDVP(), TDVP2(; trunc = truncrank(10)), BUG(; trunc = truncrank(10))] From cebc8472f925b6771cf268a458923ee3f3558708 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 14 Aug 2026 12:18:52 +0200 Subject: [PATCH 10/29] add to changelog --- docs/src/changelog.md | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 1788564bd..793027cb0 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -41,6 +41,16 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Changed +- `timestep`/`timestep!`/`time_evolve`/`time_evolve!` now return `(ψ, envs, ϵ)` instead of + `(ψ, envs)`. `ϵ` is the truncation error, the norm of the singular values discarded by the local + gauge factorisations, summed in squares over the bonds of a step and over the iterations, so that `ϵ²` is the discarded/truncated weight. It is nonzero for `TDVP2`, for `BUG` with + a `trunc`, and for `TDVP` with a bond expansion. Plain single-site `TDVP` never truncates and + returns `0`. Existing `ψ, envs = timestep(...)` call sites keep working! `time_evolve` now + also logs the per-step error at `verbosity ≥ 3` and the total at `verbosity ≥ 2` (previously a + hardcoded `0`), under the algorithm's own name rather than always `"TDVP"`. +- The meaning of every reported error and tolerance is now documented, and the manual has a new + [Errors and accuracy](@ref) section covering ground states, time evolution and excitations + separately. The quantities themselves are unchanged, but what they represent is clarified. - 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, @@ -99,6 +109,16 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Performance +- `TDVP2` now performs its two-site split through the shared `gauge2!` (as two-site DMRG already + did), which removes two sources of waste per local update: + - its right-to-left sweep installed the two sites in the left-to-right order, which made the + lazy orthogonality-view cache re-derive `AR` at the bond from the pre-update tensor, only to + overwrite it on the next install. `gauge2!` installs in sweep order, dropping that redundant + right-orthogonalisation per bond. + - it unconditionally complexified the bond tensor, so for a real-valued state (real Hamiltonian + in imaginary time) every local update allocated a complex copy that the state's own storage + then converted straight back to real. `gauge2!` only complexifies when the state is complex. + ## [0.13.11](https://github.com/QuantumKitHub/MPSKit.jl/compare/v0.13.10...v0.13.11) - 2026-05-04 ### Added From 086f3ce21303c5858060679c7c25660c42d0a330 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 19 Aug 2026 09:56:01 +0200 Subject: [PATCH 11/29] define algorithminfo struct --- src/MPSKit.jl | 2 + src/utility/algorithminfo.jl | 134 +++++++++++++++++++++++++++++++++++ 2 files changed, 136 insertions(+) create mode 100644 src/utility/algorithminfo.jl diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 29077215b..b386d7bce 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -35,6 +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 AlgorithmInfo export TDVP, TDVP2, BUG, WI, WII, TaylorCluster export changebonds, changebonds! export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand @@ -104,6 +105,7 @@ include("utility/allocator.jl") include("utility/logging.jl") using .IterativeLoggers include("utility/iterativesolvers.jl") +include("utility/algorithminfo.jl") include("utility/styles.jl") include("utility/periodicarray.jl") diff --git a/src/utility/algorithminfo.jl b/src/utility/algorithminfo.jl new file mode 100644 index 000000000..1f511effd --- /dev/null +++ b/src/utility/algorithminfo.jl @@ -0,0 +1,134 @@ +""" +$(TYPEDEF) + +Information about how an algorithm arrived at its result, returned as the last value by +[`find_groundstate`](@ref), [`leading_boundary`](@ref), [`approximate`](@ref), [`timestep`](@ref) +and [`time_evolve`](@ref). + +Algorithms in MPSKit produce genuinely different error measures, and a single bare number cannot +carry that distinction. This struct keeps the return signature uniform while naming each quantity, +so that nothing is promised that an algorithm does not actually compute. + +## Convergence + + - `converged`: whether the algorithm reached its stopping criterion, or `nothing` for algorithms + that do not iterate to a fixed point. + - `normres`: the quantity compared against the algorithm's `tol`, or `nothing` when there is none. + Which measure this is depends on the algorithm. + - `numiter`: number of iterations. + +## Truncation + +Both truncation fields are built from the same per-factorisation quantity, namely the 2-norm of the +singular values a single local factorisation discarded, but aggregate it differently, because no +single aggregation answers every question: + + - `ϵ_max`: the largest of them. Independent of system size and iteration count, and on the same + scale as the per-bond threshold set by `trunc`, so this is what to compare against a tolerance + or between runs. It is also what gets logged. + - `ϵ_total`: all of them summed up in squares. This grows with system size and iteration count. + - `numtrunc`: how many local factorisations actually discarded anything. + +These are `0` for an algorithm that never truncates, which does not mean the result is exact. Rather, +it means this particular error channel is absent. See the manual on [Errors and accuracy](@ref) +for what is *not* measured here. +""" +struct AlgorithmInfo{T <: Real} + converged::Union{Bool, Nothing} + normres::Union{T, Nothing} + ϵ_max::T + ϵ_total::T + numtrunc::Int + numiter::Int +end + +""" + AlgorithmInfo(; converged, normres, truncation, numiter) + +Keyword constructor, with every field defaulting to "not produced by this algorithm": no +convergence notion, and nothing truncated. `truncation` accepts a +[`TruncationAccumulator`](@ref), or is left out when the algorithm does not truncate. +""" +function AlgorithmInfo(; + converged = nothing, normres = nothing, truncation = nothing, numiter::Int = 1 + ) + T = _info_scalartype(normres, truncation) + acc = isnothing(truncation) ? TruncationAccumulator(T) : truncation + return AlgorithmInfo{T}( + converged, isnothing(normres) ? nothing : convert(T, normres), + convert(T, acc.ϵ_max), convert(T, sqrt(acc.ϵ_sq)), acc.numtrunc, numiter + ) +end + +_info_scalartype(::Nothing, ::Nothing) = Float64 +_info_scalartype(normres, ::Nothing) = float(typeof(normres)) +_info_scalartype(::Nothing, acc) = _acc_type(acc) +_info_scalartype(normres, acc) = promote_type(float(typeof(normres)), _acc_type(acc)) + +""" + TruncationAccumulator{T} + +Collects the per-factorisation truncation errors of a sweep. Algorithms push errors in as they are +produced with [`push_error!`](@ref) and never decide how they aggregate. +The latter is [`AlgorithmInfo`](@ref)'s job. +""" +mutable struct TruncationAccumulator{T <: Real} + ϵ_max::T + ϵ_sq::T + numtrunc::Int +end +function TruncationAccumulator(::Type{T}) where {T <: Real} + return TruncationAccumulator{T}(zero(T), zero(T), 0) +end +TruncationAccumulator(ψ) = TruncationAccumulator(_truncation_scalartype(ψ)) +_truncation_scalartype(ψ) = real(scalartype(ψ)) +_acc_type(::TruncationAccumulator{T}) where {T} = T + +""" + push_error!(acc::TruncationAccumulator, ϵ) -> acc + +Record the error of a single local factorisation. A factorisation that discarded nothing +(a QR gauge, or a truncation that kept everything) is not counted. +""" +function push_error!(acc::TruncationAccumulator, ϵ) + iszero(ϵ) && return acc + acc.ϵ_max = max(acc.ϵ_max, ϵ) + acc.ϵ_sq += ϵ^2 + acc.numtrunc += 1 + return acc +end + +# combining infos follows the same rules as combining the per-bond errors within one of them: +# worst case for `ϵ_max`, sum of squares for `ϵ_total`, and the later convergence verdict wins +# assumes type stability of the scalars and an ordering of receiving this info +function _combine(a::AlgorithmInfo, b::AlgorithmInfo) + return AlgorithmInfo( + b.converged, b.normres, + max(a.ϵ_max, b.ϵ_max), sqrt(a.ϵ_total^2 + b.ϵ_total^2), + a.numtrunc + b.numtrunc, a.numiter + b.numiter + ) +end + +function Base.show(io::IO, ::MIME"text/plain", info::AlgorithmInfo) + println(io, "AlgorithmInfo:") + if !isnothing(info.converged) + println(io, " converged = ", info.converged, " after ", info.numiter, " iterations") + else + println(io, " ", info.numiter, " iteration", info.numiter == 1 ? "" : "s") + end + isnothing(info.normres) || println(io, " normres = ", info.normres) + if info.numtrunc > 0 + println(io, " ϵ_max = ", info.ϵ_max, "\t(largest single factorization)") + println(io, " ϵ_total = ", info.ϵ_total, "\t(quadrature over ", info.numtrunc, ")") + else + println(io, " no truncation") + end + return nothing +end +function Base.show(io::IO, info::AlgorithmInfo) + return print( + io, "AlgorithmInfo(converged = ", info.converged, ", normres = ", info.normres, + ", ϵ_max = ", info.ϵ_max, ", ϵ_total = ", info.ϵ_total, + ", numtrunc = ", info.numtrunc, ", numiter = ", info.numiter, ")" + ) +end From 190686ceed6d44d6670456e73a2b85f703aa9341 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 19 Aug 2026 09:57:35 +0200 Subject: [PATCH 12/29] update ground state docstrings --- src/algorithms/approximate/approximate.jl | 34 +++++++++---------- .../groundstate/find_groundstate.jl | 18 +++++----- src/algorithms/statmech/leading_boundary.jl | 11 +++--- 3 files changed, 31 insertions(+), 32 deletions(-) diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index f48942ecd..9d450b5a7 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -1,11 +1,11 @@ @doc """ - approximate(ψ₀, (O, ψ), [environments]; kwargs...) -> (ψ, environments, ϵ) - approximate(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, ϵ) - approximate!(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, ϵ) - approximate(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ) - approximate!(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ) - approximate((O, ψ), algorithm) -> (ψ′, ϵ) - approximate!(ψ₀, (O, ψ), algorithm) -> (ψ, ϵ) + approximate(ψ₀, (O, ψ), [environments]; kwargs...) -> (ψ, environments, info) + approximate(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, info) + approximate!(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, info) + approximate(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, info) + approximate!(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, info) + approximate((O, ψ), algorithm) -> (ψ′, info) + approximate!(ψ₀, (O, ψ), algorithm) -> (ψ, info) Compute an approximation to the application of an operator `O` to the state `ψ` in the form of an MPS, using initial guess `ψ₀`. If only a state `ψ` is supplied instead of the `(O, ψ)` pair, @@ -39,16 +39,16 @@ struct itself instead (e.g. `DMRG(; tol, maxiter, verbosity)`). - `ψ`: the approximated state - `environments`: environments corresponding to the result (not returned by `Zipup`, which uses none) -- `ϵ::Float64`: an error measure whose meaning depends on the algorithm: - - for the iterative algorithms (`DMRG`, `DMRG2`, `IDMRG`, `IDMRG2`, `VOMPS`) it is the final - convergence error, i.e. the Galerkin error compared against `tol`. This is the same quantity - [`find_groundstate`](@ref) returns, measuring distance from the variational fixed point. - - for [`Zipup`](@ref) it is instead a truncation error, i.e. the largest 2-norm of the discarded - singular values over all bonds and sweeps. `Zipup` is a single non-iterative sweep, so there is - no convergence measure to report and no `tol` to compare against. +- `info::AlgorithmInfo`: how the algorithm arrived there. Which of its fields are populated depends + on the algorithm: + - the iterative algorithms (`DMRG`, `DMRG2`, `IDMRG`, `IDMRG2`, `VOMPS`) fill `converged`, + `normres` and `numiter`. `normres` is the Galerkin error compared against `tol`, the same + quantity [`find_groundstate`](@ref) reports, measuring distance from the variational fixed point. + - [`Zipup`](@ref) is a single non-iterative sweep, so it has no convergence measure at all: + `converged` and `normres` are `nothing`, and it fills the truncation fields instead. - The two are not comparable, and a small `ϵ` means different things in each case. See the manual on - the `ϵ` convention under [The error convention](@ref). + See [`AlgorithmInfo`](@ref) for the full list and [The error convention](@ref) in the manual for + why a convergence measure and a truncation error are not comparable quantities. # Algorithms @@ -56,7 +56,7 @@ Each algorithm below only supports a subset of the general interface. Check this picking one — in particular, note that **only `DMRG`/`DMRG2` accept a bare state `ψ`**; the infinite algorithms always require an explicit `(O, ψ)` tuple, and **`VOMPS` has no in-place `approximate!`** at all. `Zipup` is a single sweep rather than an iterative optimization, so it uses -no environments and returns `(ψ, ϵ)`; its `ψ₀` is a write destination, not an initial guess, and it +no environments and returns `(ψ, info)`; its `ψ₀` is a write destination, not an initial guess, and it may be omitted. | Algorithm | Scheme | State `ψ₀` | bare `ψ` allowed? | `approximate!` | diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index bc76686ab..9a350d19b 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/src/algorithms/groundstate/find_groundstate.jl @@ -1,6 +1,6 @@ """ - find_groundstate(ψ₀, H, [environments]; kwargs...) -> (ψ, environments, ϵ) - find_groundstate(ψ₀, H, algorithm, [environments]) -> (ψ, environments, ϵ) + find_groundstate(ψ₀, H, [environments]; kwargs...) -> (ψ, environments, info) + find_groundstate(ψ₀, H, algorithm, [environments]) -> (ψ, environments, info) Compute the ground state for Hamiltonian `H` with initial guess `ψ₀`. If no `algorithm` is specified, one is selected automatically from the type of `ψ₀` and the supplied keywords @@ -38,13 +38,11 @@ low-bond-dimension initial guess such as a product state. - `ψ::AbstractMPS`: converged ground state - `environments`: environments corresponding to the converged state -- `ϵ::Float64`: final convergence error upon terminating the algorithm, i.e. the quantity compared - against the algorithm's `tol`. It measures distance from a variational fixed point and is not - a truncation error. See the manual under [The error convention](@ref) for more information. - Which measure it is depends on the algorithm: the sweeping algorithms - ([`DMRG`](@ref), [`DMRG2`](@ref), [`VUMPS`](@ref), [`IDMRG`](@ref), [`IDMRG2`](@ref)) - report the Galerkin error, whereas [`GradientGrassmann`](@ref) reports - the norm of the Riemannian gradient from its optimizer. +- `info::AlgorithmInfo`: how the algorithm terminated. `info.normres` is the quantity compared + against `tol` and `info.converged` says whether it got there. Which measure `normres` is depends + on the algorithm. A truncating algorithm additionally fills `info.ϵ_max`/`info.ϵ_total` + with what its final sweep discarded. See [`AlgorithmInfo`](@ref), + and [The error convention](@ref) in the manual. # Examples @@ -63,7 +61,7 @@ julia> H = FiniteMPOHamiltonian(lattice, ((i, i + 1) => -(X ⊗ X) for i in 1:(L julia> ψ₀ = FiniteMPS(ones(Float64, (ℂ^2)^L)); -julia> ψ, envs, ϵ = find_groundstate(ψ₀, H; verbosity = 0, trunc = truncrank(16)); +julia> ψ, envs, info = find_groundstate(ψ₀, H; verbosity = 0, trunc = truncrank(16)); julia> round(real(expectation_value(ψ, H)); digits = 4) -4.7588 diff --git a/src/algorithms/statmech/leading_boundary.jl b/src/algorithms/statmech/leading_boundary.jl index de1694cdf..8124e4eb8 100644 --- a/src/algorithms/statmech/leading_boundary.jl +++ b/src/algorithms/statmech/leading_boundary.jl @@ -1,6 +1,6 @@ @doc """ - leading_boundary(ψ₀, O, [environments]; kwargs...) -> (ψ, environments, ϵ) - leading_boundary(ψ₀, O, algorithm, environments) -> (ψ, environments, ϵ) + leading_boundary(ψ₀, O, [environments]; kwargs...) -> (ψ, environments, info) + leading_boundary(ψ₀, O, algorithm, environments) -> (ψ, environments, info) Compute the leading boundary MPS for operator `O` with initial guess `ψ`. If not specified, an optimization algorithm will be attempted based on the supplied keywords. @@ -22,9 +22,10 @@ optimization algorithm will be attempted based on the supplied keywords. - `ψ::AbstractMPS`: converged leading boundary MPS - `environments`: environments corresponding to the converged boundary -- `ϵ::Float64`: final convergence error upon terminating the algorithm, i.e. the Galerkin error. - It is not a truncation error; see [`find_groundstate`](@ref) and the manual on the - `ϵ` convention under [The error convention](@ref), and [Ground-state accuracy](@ref). +- `info::AlgorithmInfo`: how the algorithm terminated; `info.normres` is the Galerkin error compared + against `tol`, and `info.converged` whether it got there. It is not a truncation error. See + [`AlgorithmInfo`](@ref), [`find_groundstate`](@ref), and the manual on the `ϵ` convention under + [The error convention](@ref) and [Ground state accuracy](@ref). """ leading_boundary # TODO: alg selector From b9012616f84f0c452f03cd7af8e1ac6a247fb6a0 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 19 Aug 2026 16:50:41 +0200 Subject: [PATCH 13/29] algorithminfo in ground states --- src/algorithms/groundstate/dmrg.jl | 30 +++++++++++++------ .../groundstate/gradient_grassmann.jl | 9 ++++-- src/algorithms/groundstate/idmrg.jl | 5 +++- src/algorithms/groundstate/vumps.jl | 4 +-- 4 files changed, 34 insertions(+), 14 deletions(-) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 940e8e503..105473d11 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -60,8 +60,9 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`approximate`](@ref). """ struct DMRG{A, F, E, G, B} <: Algorithm - "convergence tolerance, compared against the Galerkin error (the tangent-space gradient - norm)" + "convergence tolerance on the Galerkin error (the tangent-space gradient norm). This acts as a + floor: the stopping test is `ϵ ≤ max(tol, maximum(ϵ_trunc))`, which reduces to `ϵ ≤ tol` + when nothing is truncated" tol::Float64 "maximal amount of iterations" @@ -172,8 +173,9 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`approximate`](@ref). """ struct DMRG2{A, G, F, B} <: Algorithm - "convergence tolerance, compared against the Galerkin error (the tangent-space gradient - norm)" + "convergence tolerance on the Galerkin error (the tangent-space gradient norm). This acts as a + floor: the stopping test is `ϵ ≤ max(tol, maximum(ϵ_trunc))`, which reduces to `ϵ ≤ tol` + when nothing is truncated" tol::Float64 "maximal amount of iterations" @@ -258,7 +260,7 @@ _sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1) inner_alg_gauge(alg::Union{DMRG, DMRG2}) = alg_gauge(alg.alg_gauge) """ - find_groundstate!(ψ, H, algorithm, [environments]) -> (ψ, environments, ϵ) + find_groundstate!(ψ, H, algorithm, [environments]) -> (ψ, environments, info) In-place version of [`find_groundstate`](@ref): optimize the finite MPS `ψ` for the Hamiltonian `H`, overwriting the input state instead of working on a copy. @@ -275,8 +277,10 @@ Currently supported for the finite-system algorithms [`DMRG`](@ref) and [`DMRG2` - `ψ::AbstractFiniteMPS`: converged ground state - `environments`: environments corresponding to the converged state -- `ϵ::Float64`: final convergence error upon terminating the algorithm, - i.e. the Galerkin error (see [`find_groundstate`](@ref)) +- `info::AlgorithmInfo`: how the algorithm terminated. `info.normres` is the Galerkin error and + `info.converged` whether it met the stopping test. A truncating gauge also fills + `info.ϵ_max`/`info.ϵ_total` with what the final sweep discarded + (see [`find_groundstate`](@ref) and [`AlgorithmInfo`](@ref)) """ function find_groundstate!( ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs = environments(ψ, H, ψ) @@ -301,10 +305,11 @@ function _find_groundstate_sweep!( ϵ_truncs = zeros(Tr, n) # local truncation error decay_rates = zeros(n) # local observed decay rate of eigensolver fwd, bwd = _sweep_ranges(alg, ψ) + local iter LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ_global, expectation_value(ψ, H, envs)) - for iter in 1:(alg.maxiter) + for outer iter in 1:(alg.maxiter) @timeit timeroutput "sweep" begin # left-to-right for pos in fwd @@ -353,7 +358,14 @@ function _find_groundstate_sweep!( end end end - return ψ, envs, ϵ_global + + acc = TruncationAccumulator(Tr) + foreach(ϵ -> push_error!(acc, ϵ), ϵ_truncs) + info = AlgorithmInfo(; + converged = ϵ_global <= max(alg.tol, maximum(ϵ_truncs)), normres = ϵ_global, + truncation = acc, numiter = iter + ) + return ψ, envs, info end function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2}, envs...; kwargs...) diff --git a/src/algorithms/groundstate/gradient_grassmann.jl b/src/algorithms/groundstate/gradient_grassmann.jl index ca323a5c7..82ec02a63 100644 --- a/src/algorithms/groundstate/gradient_grassmann.jl +++ b/src/algorithms/groundstate/gradient_grassmann.jl @@ -81,7 +81,7 @@ end function find_groundstate( ψ::S, H, alg::GradientGrassmann, envs::P = environments(ψ, H, ψ) - )::Tuple{S, P, Float64} where {S, P} + )::Tuple{S, P, AlgorithmInfo{Float64}} where {S, P} !isa(ψ, FiniteMPS) || dim(ψ.C[end]) == 1 || @warn "This is not fully supported - split the mps up in a sum of mps's and optimize separately" normalize!(ψ) @@ -122,5 +122,10 @@ function find_groundstate( @infov 4 timeroutput end - return x, envs, normgradhistory[end] + normres = normgradhistory[end] + info = AlgorithmInfo(; + converged = normres <= alg.method.gradtol, normres, + numiter = size(normgradhistory, 1) + ) + return x, envs, info end diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index b8732417d..67055a691 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -137,7 +137,10 @@ function _find_groundstate_idmrg(mps, operator, alg::alg_type, envs) where {alg_ alg_gauge = adapt_solver(alg.alg_gauge; iter = it.state.iter, g_global = it.state.ϵ) ψ′ = InfiniteMPS(it.state.mps.AR; alg_gauge.tol, alg_gauge.maxiter) envs = recalculate!(it.state.envs, ψ′, it.state.operator, ψ′) - return ψ′, envs, it.state.ϵ + info = AlgorithmInfo(; + converged = it.state.ϵ <= alg.tol, normres = it.state.ϵ, numiter = it.state.iter + ) + return ψ′, envs, info end end diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index a7c6d5419..5ac67deb5 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -84,12 +84,12 @@ function dominant_eigsolve( if ϵ ≤ alg.tol @infov 4 timeroutput @infov 2 logfinish!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) - return mps, envs, ϵ + return mps, envs, AlgorithmInfo(; converged = true, normres = ϵ, numiter = it.iter) end if it.iter ≥ alg.maxiter @infov 4 timeroutput @warnv 1 logcancel!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) - return mps, envs, ϵ + return mps, envs, AlgorithmInfo(; converged = false, normres = ϵ, numiter = it.iter) end @infov 3 logiter!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) end From e0e695765f064873722c72cbc5d140b2afe636b7 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 19 Aug 2026 16:51:54 +0200 Subject: [PATCH 14/29] algorithminfo in statmech --- src/algorithms/statmech/gradient_grassmann.jl | 7 ++++++- src/algorithms/statmech/idmrg.jl | 4 ++-- src/algorithms/statmech/vomps.jl | 4 ++-- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/algorithms/statmech/gradient_grassmann.jl b/src/algorithms/statmech/gradient_grassmann.jl index 8084cbda2..3a07226bb 100644 --- a/src/algorithms/statmech/gradient_grassmann.jl +++ b/src/algorithms/statmech/gradient_grassmann.jl @@ -19,5 +19,10 @@ function leading_boundary( alg.finalize!, isometrictransport = true ) - return x, envs, normgradhistory[end] + normres = normgradhistory[end] + info = AlgorithmInfo(; + converged = normres <= alg.method.gradtol, normres, + numiter = size(normgradhistory, 1) + ) + return x, envs, info end diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 8e283dd6e..09cd17d21 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -59,7 +59,7 @@ function leading_boundary( ψ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter) recalculate!(envs, ψ, operator, ψ) - return ψ, envs, ϵ + return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, normres = ϵ, numiter = iter) end function leading_boundary( @@ -195,5 +195,5 @@ function leading_boundary( ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter) recalculate!(envs, ψ, operator, ψ) - return ψ, envs, ϵ + return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, normres = ϵ, numiter = iter) end diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl index dfbd81d1b..5733abed1 100644 --- a/src/algorithms/statmech/vomps.jl +++ b/src/algorithms/statmech/vomps.jl @@ -74,11 +74,11 @@ function dominant_eigsolve( for (mps, envs, ϵ) in it if ϵ ≤ alg.tol @infov 2 logfinish!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) - return mps, envs, ϵ + return mps, envs, AlgorithmInfo(; converged = true, normres = ϵ, numiter = it.iter) end if it.iter ≥ alg.maxiter @warnv 1 logcancel!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) - return mps, envs, ϵ + return mps, envs, AlgorithmInfo(; converged = false, normres = ϵ, numiter = it.iter) end @infov 3 logiter!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) end From 936baf78b6d8d61fc1c8847734696b187c55f892 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 19 Aug 2026 17:41:20 +0200 Subject: [PATCH 15/29] return trunc errors in idmrg2 groundstates and stat mech --- src/algorithms/approximate/fvomps.jl | 10 +++++--- src/algorithms/groundstate/idmrg.jl | 37 +++++++++++++++++----------- src/algorithms/statmech/idmrg.jl | 15 +++++++---- 3 files changed, 39 insertions(+), 23 deletions(-) diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl index 47dcfdfef..8dfa5cd19 100644 --- a/src/algorithms/approximate/fvomps.jl +++ b/src/algorithms/approximate/fvomps.jl @@ -1,11 +1,12 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environments(ψ, _environment_args(Oϕ)...)) allocator = default_allocator(ψ, SerialScheduler()) ϵ::Float64 = 2 * alg.tol + local iter log = IterLog("DMRG2") LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ) - for iter in 1:(alg.maxiter) + for outer iter in 1:(alg.maxiter) ϵ = 0.0 for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1] AC2′ = AC2_projection(pos, ψ, Oϕ, envs; alg.backend, allocator) @@ -33,17 +34,18 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environment end end - return ψ, envs, ϵ + return ψ, envs, AlgorithmInfo(; converged = ϵ < alg.tol, normres = ϵ, numiter = iter) end function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments(ψ, _environment_args(Oϕ)...)) allocator = default_allocator(ψ, SerialScheduler()) ϵ::Float64 = 2 * alg.tol + local iter log = IterLog("DMRG") LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ) - for iter in 1:(alg.maxiter) + for outer iter in 1:(alg.maxiter) ϵ = 0.0 for pos in [1:(length(ψ) - 1); length(ψ):-1:2] AC′ = AC_projection(pos, ψ, Oϕ, envs; alg.backend, allocator) @@ -68,5 +70,5 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments end end - return ψ, envs, ϵ + return ψ, envs, AlgorithmInfo(; converged = ϵ < alg.tol, normres = ϵ, numiter = iter) end diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 67055a691..8694f98e5 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -79,16 +79,18 @@ struct IDMRGState{S, O, E, T, A} envs::E iter::Int ϵ::Float64 # TODO: Could be any <:Real + truncation::TruncationAccumulator{Float64} # of the most recent sweep only energy::T timeroutput::TimerOutput allocator::A end function IDMRGState{T}( - mps::S, operator::O, envs::E, iter::Int, ϵ::Float64, energy, + mps::S, operator::O, envs::E, iter::Int, ϵ::Float64, + truncation::TruncationAccumulator{Float64}, energy, timeroutput::TimerOutput, allocator::A, ) where {S, O, E, T, A} return IDMRGState{S, O, E, T, A}( - mps, operator, envs, iter, ϵ, T(energy), timeroutput, allocator + mps, operator, envs, iter, ϵ, truncation, T(energy), timeroutput, allocator ) end @@ -116,7 +118,8 @@ function _find_groundstate_idmrg(mps, operator, alg::alg_type, envs) where {alg_ end end - state = IDMRGState(mps, operator, envs, iter, ϵ, E, timeroutput, allocator) + acc = TruncationAccumulator(Float64) + state = IDMRGState(mps, operator, envs, iter, ϵ, acc, E, timeroutput, allocator) it = IterativeSolver(alg, state) return LoggingExtras.withlevel(; alg.verbosity) do @@ -138,7 +141,8 @@ function _find_groundstate_idmrg(mps, operator, alg::alg_type, envs) where {alg_ ψ′ = InfiniteMPS(it.state.mps.AR; alg_gauge.tol, alg_gauge.maxiter) envs = recalculate!(it.state.envs, ψ′, it.state.operator, ψ′) info = AlgorithmInfo(; - converged = it.state.ϵ <= alg.tol, normres = it.state.ϵ, numiter = it.state.iter + converged = it.state.ϵ <= alg.tol, normres = it.state.ϵ, + truncation = it.state.truncation, numiter = it.state.iter ) return ψ′, envs, info end @@ -148,7 +152,8 @@ function Base.iterate( it::IterativeSolver{alg_type}, state::IDMRGState{<:Any, <:Any, <:Any, T} = it.state ) where {alg_type <: Union{<:IDMRG, <:IDMRG2}, T} timeroutput = state.timeroutput - mps, envs, C_old, E_new = @timeit timeroutput "localupdate" localupdate_step!(it, state) + acc = TruncationAccumulator(Float64) # fresh each sweep, what the state carries is what's last discarded + mps, envs, C_old, E_new = @timeit timeroutput "localupdate" localupdate_step!(it, state, acc) # error criterion C = mps.C[0] @@ -169,14 +174,14 @@ function Base.iterate( # update state it.state = IDMRGState{T}( - mps, state.operator, envs, state.iter + 1, ϵ, E_new, timeroutput, state.allocator, + mps, state.operator, envs, state.iter + 1, ϵ, acc, E_new, timeroutput, state.allocator, ) return (mps, envs, ϵ, ΔE), it.state end function localupdate_step!( - it::IterativeSolver{<:IDMRG}, state + it::IterativeSolver{<:IDMRG}, state, acc ) alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = state.iter, g_global = state.ϵ) return _localupdate_sweep_idmrg!( @@ -186,12 +191,12 @@ function localupdate_step!( end function localupdate_step!( - it::IterativeSolver{<:IDMRG2}, state + it::IterativeSolver{<:IDMRG2}, state, acc ) alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = state.iter, g_global = state.ϵ) return _localupdate_sweep_idmrg2!( state.mps, state.operator, state.envs, alg_eigsolve, - it.trunc, it.alg_svd, state.timeroutput; + it.trunc, it.alg_svd, state.timeroutput, acc; it.backend, state.allocator, ) end @@ -235,7 +240,7 @@ function _localupdate_sweep_idmrg!( end function _localupdate_sweep_idmrg2!( - ψ, H, envs, alg_eigsolve, alg_trunc, alg_svd, timeroutput::TimerOutput; + ψ, H, envs, alg_eigsolve, alg_trunc, alg_svd, timeroutput::TimerOutput, acc::TruncationAccumulator; backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator() ) # @timeit wraps its body in try-finally, which is a new lexical scope: declare locals @@ -249,7 +254,8 @@ function _localupdate_sweep_idmrg2!( _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) end @timeit timeroutput "svd_trunc" begin - al, c, ar = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[pos] = al @@ -272,7 +278,8 @@ function _localupdate_sweep_idmrg2!( _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) end @timeit timeroutput "svd_trunc" begin - al, c, ar = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[end] = al @@ -300,7 +307,8 @@ function _localupdate_sweep_idmrg2!( _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) end @timeit timeroutput "svd_trunc" begin - al, c, ar = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[pos] = al @@ -324,7 +332,8 @@ function _localupdate_sweep_idmrg2!( E, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) end @timeit timeroutput "svd_trunc" begin - al, c, ar = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[end] = al diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 09cd17d21..5bff1e922 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -74,6 +74,7 @@ function leading_boundary( LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ) for outer iter in 1:(alg.maxiter) + acc = TruncationAccumulator(ψ) # fresh each sweep, reported truncation is that of the final sweep alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = ϵ) C_current = ψ.C[:, 0] @@ -84,7 +85,8 @@ function leading_boundary( _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve) for row in 1:size(ψ, 1) - al, c, ar = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[row + 1, site] = al @@ -108,7 +110,8 @@ function leading_boundary( _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve) for row in 1:size(ψ, 1) - al, c, ar = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[row + 1, site] = al @@ -133,7 +136,8 @@ function leading_boundary( _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve) for row in 1:size(ψ, 1) - al, c, ar = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[row + 1, site] = al @@ -155,7 +159,8 @@ function leading_boundary( _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve) for row in 1:size(ψ, 1) - al, c, ar = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[row + 1, end] = al @@ -195,5 +200,5 @@ function leading_boundary( ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter) recalculate!(envs, ψ, operator, ψ) - return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, normres = ϵ, numiter = iter) + return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, normres = ϵ, truncation = acc, numiter = iter) end From 2fa70b53cea5d8db546a4eab7f9db7565ddc2e5f Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 19 Aug 2026 17:54:55 +0200 Subject: [PATCH 16/29] algorithminfo in approximate (with idmrg now returning trunc error) --- src/algorithms/approximate/approximate.jl | 3 +++ src/algorithms/approximate/idmrg.jl | 18 ++++++++----- src/algorithms/approximate/vomps.jl | 4 +-- src/algorithms/approximate/zipup.jl | 31 ++++++++++++----------- 4 files changed, 33 insertions(+), 23 deletions(-) diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index 9d450b5a7..0115afe84 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -44,6 +44,9 @@ struct itself instead (e.g. `DMRG(; tol, maxiter, verbosity)`). - the iterative algorithms (`DMRG`, `DMRG2`, `IDMRG`, `IDMRG2`, `VOMPS`) fill `converged`, `normres` and `numiter`. `normres` is the Galerkin error compared against `tol`, the same quantity [`find_groundstate`](@ref) reports, measuring distance from the variational fixed point. + - the two-site ones (`DMRG2`, `IDMRG2`) which involve a truncated SVD fill the truncation + fields with what their final sweep discarded, i.e. what the returned state is still + throwing away per sweep rather than what the early, unconverged sweeps did. - [`Zipup`](@ref) is a single non-iterative sweep, so it has no convergence measure at all: `converged` and `normres` are `nothing`, and it fills the truncation fields instead. diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl index 7503ac76d..252df7828 100644 --- a/src/algorithms/approximate/idmrg.jl +++ b/src/algorithms/approximate/idmrg.jl @@ -60,7 +60,7 @@ function approximate!( copy!(ψ, ψ′) # ensure output destination is unchanged recalculate!(envs, ψ, toapprox) - return ψ, envs, ϵ + return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, normres = ϵ, numiter = iter) end function approximate!( @@ -77,6 +77,7 @@ function approximate!( LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ) for outer iter in 1:(alg.maxiter) + acc = TruncationAccumulator(ψ) # fresh each sweep, reported truncation is the final sweep's C_current = ψ.C[:, 0] # sweep from left to right @@ -86,7 +87,8 @@ function approximate!( CartesianIndex(row, site), ψ, toapprox, envs; kind = :ACAR, alg.backend, allocator ) - al, c, ar = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[row + 1, site] = al @@ -107,7 +109,8 @@ function approximate!( CartesianIndex(row, size(ψ, 2)), ψ, toapprox, envs; kind = :ALAC, alg.backend, allocator ) - al, c, ar = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[row + 1, end] = al @@ -132,7 +135,8 @@ function approximate!( CartesianIndex(row, site), ψ, toapprox, envs; kind = :ALAC, alg.backend, allocator ) - al, c, ar = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[row + 1, site] = al @@ -152,7 +156,8 @@ function approximate!( CartesianIndex(row, 0), ψ, toapprox, envs; kind = :ACAR, alg.backend, allocator ) - al, c, ar = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) + al, c, ar, ϵ_trunc = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) + push_error!(acc, ϵ_trunc) normalize!(c) ψ.AL[row, end] = al @@ -193,5 +198,6 @@ function approximate!( copy!(ψ, ψ′) # ensure output destination is unchanged recalculate!(envs, ψ, toapprox) - return ψ, envs, ϵ + info = AlgorithmInfo(; converged = ϵ <= alg.tol, normres = ϵ, truncation = acc, numiter = iter) + return ψ, envs, info end diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl index 77e6c8cde..22bdb675f 100644 --- a/src/algorithms/approximate/vomps.jl +++ b/src/algorithms/approximate/vomps.jl @@ -36,11 +36,11 @@ function _approximate_vomps(mps, toapprox, alg::VOMPS, envs) for (mps, envs, ϵ) in it if ϵ ≤ alg.tol @infov 2 logfinish!(log, it.iter, ϵ) - return mps, envs, ϵ + return mps, envs, AlgorithmInfo(; converged = true, normres = ϵ, numiter = it.iter) end if it.iter ≥ alg.maxiter @warnv 1 logcancel!(log, it.iter, ϵ) - return mps, envs, ϵ + return mps, envs, AlgorithmInfo(; converged = false, normres = ϵ, numiter = it.iter) end @infov 3 logiter!(log, it.iter, ϵ) end diff --git a/src/algorithms/approximate/zipup.jl b/src/algorithms/approximate/zipup.jl index c178c9db5..e4c0b8eb5 100644 --- a/src/algorithms/approximate/zipup.jl +++ b/src/algorithms/approximate/zipup.jl @@ -6,13 +6,13 @@ followed by a zip-down sweep in the opposite direction. The MPO and MPS are cont time, and the enlarged virtual bond is truncated immediately. The sweep direction is selected by `left_to_right`. - approximate((O, ϕ), alg::Zipup) -> ψ, ϵ - approximate!(ψ, (O, ϕ), alg::Zipup) -> ψ, ϵ + approximate((O, ϕ), alg::Zipup) -> ψ, info + approximate!(ψ, (O, ϕ), alg::Zipup) -> ψ, info Contrary to the variational algorithms, this algorithm requires no initial guess: the in-place version simply uses `ψ` as the destination of the sweep, overwriting its contents, and may alias `ϕ`. The out-of-place version allocates a destination with the promoted scalar type of `O` and `ϕ`. -Both return the truncation error `ϵ` alongside the approximated state. +Both return an [`AlgorithmInfo`](@ref) alongside the approximated state which contains the truncation information. # Constructors @@ -87,8 +87,8 @@ function approximate(Oϕ::Tuple{Any, <:FiniteMPS}, alg::Zipup) end @doc """ - zip_left_right!(ψ, O, ϕ, alg_zipup, [alg_zipdown]) -> ψ, ϵ - zip_right_left!(ψ, O, ϕ, alg_zipup, [alg_zipdown]) -> ψ, ϵ + zip_left_right!(ψ, O, ϕ, alg_zipup, [alg_zipdown]) -> ψ, info + zip_right_left!(ψ, O, ϕ, alg_zipup, [alg_zipdown]) -> ψ, info Contract the MPO `O` with the MPS `ϕ` in a single sweep, truncating the enlarged virtual bond at every site with `alg_zipup`, and write the result into `ψ`. `zip_left_right!` zips up from left to right, @@ -96,8 +96,9 @@ site with `alg_zipup`, and write the result into `ψ`. `zip_left_right!` zips up opposite direction imposes a final truncation with `alg_zipdown` in a locally gauged basis, leaving the gauge center of `ψ` at the far end. The destination may alias `ϕ`. -Also returns the truncation error `ϵ`, the largest 2-norm of the discarded singular values over all -bonds and both sweeps, i.e. the worst single bond. +Also returns an [`AlgorithmInfo`](@ref) describing the truncation. Being a single sweep +rather than an iterative optimisation, there is no convergence measure, +so `converged` and `normres` are `nothing`. """ zip_left_right! @doc (@doc zip_left_right!) zip_right_left! @@ -115,14 +116,14 @@ function zip_left_right!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown A = storagetype(eltype(ψ)) Fₗ = fuser(A, left_virtualspace(Aϕs[1]), left_virtualspace(O, 1)) - ϵ = zero(real(scalartype(ψ))) + acc = TruncationAccumulator(ψ) # zip up from left to right, leaving the gauge center on the last site for i in 1:(N - 1) Aᶻ = _fuse_mpo_mps_left(O[i], Aϕs[i], Fₗ) AL, Fₗ, ϵᵢ = left_gauge(Aᶻ, alg_zipup) # right factor doubles as the next left fuser ψ.ALs[i] = AL - ϵ = max(ϵ, ϵᵢ) + push_error!(acc, ϵᵢ) end Fᵣ = fuser(A, right_virtualspace(Aϕs[N]), right_virtualspace(O, N)) ψ.ACs[N] = _fuse_mpo_mps(O[N], Aϕs[N], Fₗ, Fᵣ) @@ -131,11 +132,11 @@ function zip_left_right!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown if !isnothing(alg_zipdown) for i in N:-1:2 ψ, ϵᵢ = right_gauge!(ψ, i, ψ.AC[i], alg_zipdown) - ϵ = max(ϵ, ϵᵢ) + push_error!(acc, ϵᵢ) end end - return ψ, ϵ + return ψ, AlgorithmInfo(; truncation = acc) end function zip_right_left!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown = nothing) @@ -149,14 +150,14 @@ function zip_right_left!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown # replaces them on the next site Vᵣ = right_virtualspace(Aϕs[N]) ⊗ right_virtualspace(O, N) Fᵣ = isomorphism(A, Vᵣ, fuse(Vᵣ)) - ϵ = zero(real(scalartype(ψ))) + acc = TruncationAccumulator(ψ) # zip up from right to left, leaving the gauge center on the first site for i in N:-1:2 Aᶻ = _fuse_mpo_mps_right(O[i], Aϕs[i], Fᵣ) Fᵣ, AR, ϵᵢ = _right_gauge_zip(Aᶻ, alg_zipup) # left factor doubles as the next right fuser ψ.ARs[i] = AR - ϵ = max(ϵ, ϵᵢ) + push_error!(acc, ϵᵢ) end # the carry is oriented such that it can simply be composed with the last local tensor Fₗ = fuser(A, left_virtualspace(Aϕs[1]), left_virtualspace(O, 1)) @@ -166,11 +167,11 @@ function zip_right_left!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown if !isnothing(alg_zipdown) for i in 1:(N - 1) ψ, ϵᵢ = left_gauge!(ψ, i, ψ.AC[i], alg_zipdown) - ϵ = max(ϵ, ϵᵢ) + push_error!(acc, ϵᵢ) end end - return ψ, ϵ + return ψ, AlgorithmInfo(; truncation = acc) end # `right_gauge` for the tensors of a right-to-left zip-up sweep: these are already partitioned across From e12b799d6e6862dd227bae1671eb32c2a468c573 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 19 Aug 2026 18:07:12 +0200 Subject: [PATCH 17/29] algorithminfo in time evolution --- src/algorithms/timestep/bug.jl | 18 ++++---- src/algorithms/timestep/tdvp.jl | 32 +++++++------- src/algorithms/timestep/time_evolve.jl | 61 +++++++++++++------------- 3 files changed, 54 insertions(+), 57 deletions(-) diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl index 4f53916c4..5c3ecc39a 100644 --- a/src/algorithms/timestep/bug.jl +++ b/src/algorithms/timestep/bug.jl @@ -25,13 +25,14 @@ To restore a maximal dimension of `D`, apply [`changebonds`](@ref) with an [`Svd !!! note By default the state is not renormalized, as the (loss of) norm accumulates useful information. In real time the squared norm lost is exactly the weight discarded by the bond cuts, - ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon^2``, with `ϵ` the truncation error - returned by [`timestep`](@ref). In imaginary time the norm also carries the physical decay of - the weight and no longer isolates the truncation. `ϵ` is exactly zero when not truncating. + ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon_{\\text{total}}^2``, with + `ϵ_total` from the [`AlgorithmInfo`](@ref) returned by [`timestep`](@ref). In imaginary time + the norm also carries the physical decay of the weight and no longer isolates the truncation. + Both reported errors are exactly zero when not truncating. Pass `normalize = true` to `timestep`/`time_evolve` to renormalize after every half-sweep instead. !!! warning - `ϵ` counts only the cuts in step (i). Step (iii) augments the basis without truncating, so within + The reported errors count only the cuts in step (i). Step (iii) augments the basis without truncating, so within a sweep the bond dimension temporarily exceeds `trunc` and the norm identity above applies to the state returned at the end of the step. @@ -180,8 +181,7 @@ function timestep!( # the sweep is serial, so a single allocator serves all local updates allocator = default_allocator(ψ, SerialScheduler()) - # discarded weight of the step, accumulated in squares over the bond cuts - ϵ² = zero(real(scalartype(ψ))) + acc = TruncationAccumulator(ψ) # left→right half-sweep (root = last site): `t → t + dt / 2` @timeit timeroutput "half-sweep" for site in 1:L @@ -189,7 +189,7 @@ function timestep!( site, Val(:right), ψ, H, alg, envs, t, h, allocator; imaginary_evolution, normalize, timeroutput ) - ϵ² += ϵ^2 + push_error!(acc, ϵ) end # right→left half-sweep (root = first site): `t + dt / 2 → t + dt` @@ -198,10 +198,10 @@ function timestep!( site, Val(:left), ψ, H, alg, envs, t + h, h, allocator; imaginary_evolution, normalize, timeroutput ) - ϵ² += ϵ^2 + push_error!(acc, ϵ) end - return ψ, envs, sqrt(ϵ²) + return ψ, envs, AlgorithmInfo(; truncation = acc) end # copying version diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 532528517..79c5d8681 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -14,13 +14,14 @@ state-preserving, as required for a consistent time evolution. By default the norm is not preserved: neither the bond expansion nor the truncation renormalizes, so the state norm keeps useful information. In real time this is exact, namely the squared norm drops by precisely the truncated ("discarded") weight, - ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon^2``, - with `ϵ` the truncation error returned by [`timestep`](@ref). In imaginary time the norm also carries - the physical decay of the weight, so it no longer isolates the truncation. Without `trunc` nothing is - discarded at all and the norm is conserved exactly in real time. + ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon_{\\text{total}}^2``, + with `ϵ_total` from the [`AlgorithmInfo`](@ref) returned by [`timestep`](@ref). In imaginary + time the norm also carries the physical decay of the weight, so it no longer isolates the + truncation. Without `trunc` nothing is discarded at all and the norm is conserved exactly in + real 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. + like a ground state search. This is independent of `imaginary_evolution`. CBE is only available for finite MPS. # Fields @@ -136,7 +137,7 @@ function _timestep_infinite( recalculate!(envs, ψ′, H) # infinite one-site TDVP has a fixed bond dimension and never truncates so nothing is discarded # the gauge-fixing residual is controlled by `tolgauge`, not reported here - return ψ′, envs, zero(real(scalartype(ψ′))) + return ψ′, envs, AlgorithmInfo() end function timestep!( @@ -155,9 +156,7 @@ function _timestep_finite!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs, allocator; imaginary_evolution::Bool, normalize::Bool ) - # discarded weight of the step, accumulated in squares over the local gauges - # stays exactly zero for QR gauging - ϵ² = zero(real(scalartype(ψ))) + acc = TruncationAccumulator(ψ) # sweep left to right for i in 1:(length(ψ) - 1) @@ -173,7 +172,7 @@ function _timestep_finite!( # 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) - ϵ² += ϵ^2 + push_error!(acc, ϵ) # 4. evolve the bond tensor backward Hc = C_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator) @@ -203,7 +202,7 @@ function _timestep_finite!( # 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) - ϵ² += ϵ^2 + push_error!(acc, ϵ) # 4. evolve the bond tensor backward Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator) @@ -220,7 +219,7 @@ function _timestep_finite!( imaginary_evolution ) - return ψ, envs, sqrt(ϵ²) + return ψ, envs, AlgorithmInfo(; truncation = acc) end """ @@ -283,8 +282,7 @@ function _timestep2_finite!( # the two-site center always has to be split back up, so the gauge is always a truncated SVD alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trunc) - # discarded weight of the step, accumulated in squares over the local gauges - ϵ² = zero(real(scalartype(ψ))) + acc = TruncationAccumulator(ψ) # sweep left to right for i in 1:(length(ψ) - 1) @@ -294,7 +292,7 @@ function _timestep2_finite!( # the norm of the discarded singular values is the truncation error _, ϵ = gauge2!(ψ, i, Val(:right), ac2′, alg_gauge; normalize) - ϵ² += ϵ^2 + push_error!(acc, ϵ) if i != (length(ψ) - 1) Hac = AC_hamiltonian(i + 1, ψ, H, ψ, envs; alg.backend, allocator) @@ -312,7 +310,7 @@ function _timestep2_finite!( ac2′ = integrate(Hac2, ac2, t + dt / 2, dt / 2, alg.integrator; imaginary_evolution) _, ϵ = gauge2!(ψ, i - 1, Val(:left), ac2′, alg_gauge; normalize) - ϵ² += ϵ^2 + push_error!(acc, ϵ) if i != 2 Hac = AC_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator) @@ -323,7 +321,7 @@ function _timestep2_finite!( end end - return ψ, envs, sqrt(ϵ²) + return ψ, envs, AlgorithmInfo(; truncation = acc) end # copying version diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 0b692b185..4c4bd0b68 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -1,6 +1,6 @@ """ - time_evolve(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ, envs, ϵ) - time_evolve!(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ₀, envs, ϵ) + time_evolve(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ, envs, info) + time_evolve!(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ₀, envs, info) Time-evolve the initial state `ψ₀` with Hamiltonian `H` over a given time span by stepping through each of the time points obtained by iterating t_span. @@ -26,13 +26,13 @@ through each of the time points obtained by iterating t_span. - `ψ`: the time-stepped state - `envs`: the updated environment manager -- `ϵ`: the truncation error accumulated over the whole evolution, i.e. the per-step errors of - [`timestep`](@ref), ``\\epsilon = \\sqrt{\\sum_{\\text{steps}} \\epsilon_i^2}``. - In real time with `normalize = false` this is exactly the norm lost to truncation over the whole - evolution, ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon^2``. This is zero for - algorithms that never truncate. See [`timestep`](@ref) for what it does and does not measure. +- `info::AlgorithmInfo`: the truncation performed over the whole evolution. + See [`AlgorithmInfo`](@ref) and [Time evolution accuracy](@ref) in the manual + for the difference and when to use which reported error measure, + and [`timestep`](@ref) for what neither measures. -The per-step error is logged at `verbosity ≥ 3` and the accumulated total at `verbosity ≥ 2`. +`ϵ_max` is logged per step at `verbosity ≥ 3` and for the whole evolution at `verbosity ≥ 2`. +The size-independent measure is used here for the same reason as the ground state algorithms. """ function time_evolve end, function time_evolve! end @@ -43,30 +43,31 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time verbosity::Int = 0, imaginary_evolution::Bool = false, normalize::Bool = false ) log = IterLog(string(nameof(typeof(alg)))) - # discarded weight over the whole evolution, accumulated in squares over the steps - ϵ² = zero(real(scalartype(ψ))) + info = AlgorithmInfo(; truncation = TruncationAccumulator(ψ), numiter = 0) LoggingExtras.withlevel(; verbosity) do @infov 2 loginit!(log, 0.0, first(t_span)) for iter in 1:(length(t_span) - 1) t = t_span[iter] dt = t_span[iter + 1] - t - ψ, envs, ϵ = $timestep(ψ, H, t, dt, alg, envs; imaginary_evolution, normalize) + ψ, envs, info_step = $timestep( + ψ, H, t, dt, alg, envs; imaginary_evolution, normalize + ) ψ, envs = alg.finalize(t, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)} - ϵ² += ϵ^2 + info = _combine(info, info_step) - # per-step truncation error; the running total is what is reported at the end - @infov 3 logiter!(log, iter, convert(Float64, ϵ), t) + # log the size-independent error measure + @infov 3 logiter!(log, iter, convert(Float64, info_step.ϵ_max), t) end - @infov 2 logfinish!(log, length(t_span), convert(Float64, sqrt(ϵ²)), t_span[end]) + @infov 2 logfinish!(log, length(t_span), convert(Float64, info.ϵ_max), t_span[end]) end - return ψ, envs, sqrt(ϵ²) + return ψ, envs, info end end """ - timestep(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ, envs, ϵ) - timestep!(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ₀, envs, ϵ) + timestep(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ, envs, info) + timestep!(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ₀, envs, info) Time-step the state `ψ₀` with Hamiltonian `H` over a given time step `dt` at time `t`, solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``. @@ -92,24 +93,22 @@ solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``. - `ψ`: the time-stepped state - `envs`: the updated environment manager -- `ϵ`: the truncation error of the step (see below) +- `info::AlgorithmInfo`: what the step truncated (see below) # Truncation error -`ϵ` is the truncation error of the step, i.e. the 2-norm of the sum of squared singular values discarded -by the local gauge factorizations. In other words, `ϵ²` is the truncated ("discarded") weight. +A step performs many local factorisations, each discarding some weight. Rather than collapse those +into one number, `info` reports both aggregations under names that say what they are: +`info.ϵ_max` is the largest single one (size-independent, comparable against `trunc` and across +runs) and `info.ϵ_total` sums them in squares. -It is nonzero only for algorithms that truncate such as [`TDVP2`](@ref), [`BUG`](@ref) with a `trunc`, and -[`TDVP`](@ref) with a bond expansion, while it is exactly `0` for one-site [`TDVP`](@ref), -which runs at fixed bond dimension. A zero `ϵ` does not mean the step was exact, but that this -particular error channel is absent. In particular, the projection and time-discretisation errors are never -included in `ϵ`. +Both are non-zero only for algorithms that truncate ([`TDVP2`](@ref), [`BUG`](@ref) with a +`trunc`, and [`TDVP`](@ref) with a bond expansion), and are exactly `0` for one-site +[`TDVP`](@ref), which runs at fixed bond dimension. A zero here does not mean the step was exact, +but that this particular error channel is absent. -In real time with `normalize = false`, `ϵ` is exactly the norm lost to truncation, -``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon^2``. - -See [Time evolution accuracy](@ref) in the manual for what the other error sources are, why the -per-bond errors combine in a squared manner, and the precise statement and caveats of the norm identity. +See [`AlgorithmInfo`](@ref) for the fields, and [Time evolution accuracy](@ref) in the manual +for the other error sources. # Examples From f75b86fde2a07de08dcd98338f65802b131d0346 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 19 Aug 2026 18:34:39 +0200 Subject: [PATCH 18/29] update and expand on tests --- test/algorithms/approximate.jl | 23 ++++---- test/algorithms/groundstate.jl | 102 ++++++++++++++++++--------------- test/algorithms/timestep.jl | 42 ++++++++++---- test/misc/multifusion.jl | 16 +++--- 4 files changed, 106 insertions(+), 77 deletions(-) diff --git a/test/algorithms/approximate.jl b/test/algorithms/approximate.jl index fdf71df46..6e8d49aa3 100644 --- a/test/algorithms/approximate.jl +++ b/test/algorithms/approximate.jl @@ -116,9 +116,10 @@ end # both sweep directions, with and without the zip-down pass for left_to_right in (true, false), trunc′ in (trunc, (notrunc(), trunc)) - got, ϵ = approximate((O, ψ), Zipup(; trunc = trunc′, left_to_right)) + got, info = approximate((O, ψ), Zipup(; trunc = trunc′, left_to_right)) @test norm(ref - got) / norm(ref) < 1.0e-10 - @test ϵ < 1.0e-10 + @test info.ϵ_max < 1.0e-10 + @test isnothing(info.converged) && isnothing(info.normres) end @test norm(ψ - ψ_copy) < 1.0e-12 @@ -133,11 +134,11 @@ end zipup_trunc = truncrank(2Dcut) & truncerror(; rtol = rtol / 10) ref_tr = changebonds(O * ψ, SvdCut(; trunc = final_trunc); normalize = false) - got_one_sweep, ϵ_one_sweep = approximate((O, ψ), Zipup(; trunc = final_trunc, left_to_right)) + got_one_sweep, info_one_sweep = approximate((O, ψ), Zipup(; trunc = final_trunc, left_to_right)) got_two_sweep, _ = approximate( (O, ψ), Zipup(; trunc = (zipup_trunc, final_trunc), left_to_right) ) - @test ϵ_one_sweep > 0 + @test info_one_sweep.ϵ_max > 0 err_one_sweep = norm(ref_tr - got_one_sweep) / norm(ref_tr) err_two_sweep = norm(ref_tr - got_two_sweep) / norm(ref_tr) @@ -150,26 +151,26 @@ end (pspace, Dspace, Dcut) in zipup_spacelist, left_to_right in (true, false) O, ψ = _random_mpo_mps(pspace, Dspace) alg = Zipup(; trunc = (truncrank(2Dcut), truncrank(Dcut)), left_to_right) - ref, ϵ_ref = approximate((O, ψ), alg) + ref, info_ref = approximate((O, ψ), alg) # empty destination dst = similar(ψ, ComplexF64) - got, ϵ = approximate!(dst, (O, ψ), alg) + got, info = approximate!(dst, (O, ψ), alg) @test got === dst @test norm(ref - got) / norm(ref) < 1.0e-12 - @test ϵ ≈ ϵ_ref + @test info.ϵ_max ≈ info_ref.ϵ_max # a destination with unrelated contents is overwritten entirely dst = FiniteMPS(rand, ComplexF64, length(ψ), pspace, oneunit(Dspace) ⊕ Dspace ⊕ Dspace) - got, ϵ = approximate!(dst, (O, ψ), alg) + got, info = approximate!(dst, (O, ψ), alg) @test norm(ref - got) / norm(ref) < 1.0e-12 - @test ϵ ≈ ϵ_ref + @test info.ϵ_max ≈ info_ref.ϵ_max # the input may serve as its own destination - got, ϵ = approximate!(ψ, (O, ψ), alg) + got, info = approximate!(ψ, (O, ψ), alg) @test got === ψ @test norm(ref - got) / norm(ref) < 1.0e-12 - @test ϵ ≈ ϵ_ref + @test info.ϵ_max ≈ info_ref.ϵ_max end @testset "Zip-up with non-trivial boundary spaces $(spacetype(pspace))" for (pspace, Dspace, _) in zipup_spacelist diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl index bbb6f12e5..b1579e8d2 100644 --- a/test/algorithms/groundstate.jl +++ b/test/algorithms/groundstate.jl @@ -27,19 +27,20 @@ verbosity_conv = 1 v₀ = variance(ψ₀, H) # test logging - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2) ) - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10), envs ) v = variance(ψ, H) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 + @test info.numtrunc == 0 # the algorithm object carries no scratch space of its own - the sweep's allocator is # obtained per solve - so re-using one across solves has to reproduce the answer @@ -54,19 +55,23 @@ verbosity_conv = 1 v₀ = variance(ψ₀, H) trunc = truncrank(floor(Int, D * 1.5)) # test logging - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H, DMRG2(; verbosity = verbosity_full, maxiter = 2, trunc) ) - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ, H, DMRG2(; verbosity = verbosity_conv, maxiter = 10, trunc), envs ) v = variance(ψ, H) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 + + @test info.numtrunc > 0 + @test info.ϵ_max > 0 + @test info.ϵ_max <= info.ϵ_total <= sqrt(info.numtrunc) * info.ϵ_max end @testset "CBEDMRG" begin @@ -77,17 +82,17 @@ verbosity_conv = 1 trunc = truncrank(D) # test logging - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_expand = expand, trunc) ) - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_expand = expand, trunc), envs ) v = variance(ψ, H) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 # the bond should have grown to the truncation target @@ -105,17 +110,17 @@ verbosity_conv = 1 trunc = truncrank(D) # test logging - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_expand = expand, trunc) ) - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 15, alg_expand = expand, trunc), envs ) v = variance(ψ, H) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 # the bond should have grown to the truncation target @@ -131,17 +136,17 @@ verbosity_conv = 1 trunc = truncrank(D) # test logging - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_gauge, trunc) ) - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_gauge, trunc), envs ) v = variance(ψ, H) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 # the bond should have grown to the truncation target @@ -155,13 +160,13 @@ verbosity_conv = 1 Random.seed!(1234) ψ_bad = bad_initial_state(H_heis, L_heis) - ψ_stuck, envs_stuck, δ_stuck = find_groundstate( + ψ_stuck, envs_stuck, info_stuck = find_groundstate( ψ_bad, H_heis, DMRG(; verbosity = verbosity_conv, maxiter = 30) ) E_stuck = real(expectation_value(ψ_stuck, H_heis, envs_stuck)) alg_gauge = DMRG3S(0.1, ExponentialDecay(0.8)) - ψ_escape, envs_escape, δ_escape = find_groundstate( + ψ_escape, envs_escape, info_escape = find_groundstate( ψ_bad, H_heis, DMRG(; verbosity = verbosity_conv, maxiter = 30, alg_gauge, trunc = truncrank(20), @@ -179,17 +184,17 @@ verbosity_conv = 1 v₀ = variance(ψ₀, H) # test logging - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H, GradientGrassmann(; verbosity = verbosity_full, maxiter = 2) ) - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ, H, GradientGrassmann(; verbosity = verbosity_conv, maxiter = 50), envs ) v = variance(ψ, H) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 @test v < v₀ && v < 1.0e-2 end end @@ -214,7 +219,7 @@ end ψ₀ = repeat(InfiniteMPS(ℙ^2, ℙ^D), unit_cell_size) H = repeat(H_ref, unit_cell_size) - ψ′, envs, δ = with_scheduler(scheduler) do + ψ′, envs, info = with_scheduler(scheduler) do # test logging ψ₁, = find_groundstate( ψ₀, H, VUMPS(; tol, verbosity = verbosity_full, maxiter = 2) @@ -224,7 +229,7 @@ end v = variance(ψ′, H, envs) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 end @@ -234,15 +239,16 @@ end H = repeat(H_ref, unit_cell_size) # test logging - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ, H, IDMRG(; tol, verbosity = verbosity_full, maxiter = 2) ) - ψ, envs, δ = find_groundstate(ψ, H, IDMRG(; tol, verbosity = verbosity_conv)) + ψ, envs, info = find_groundstate(ψ, H, IDMRG(; tol, verbosity = verbosity_conv)) v = variance(ψ, H, envs) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 + @test info.numtrunc == 0 @test v < v₀ @test v < 1.0e-2 end @@ -254,19 +260,23 @@ end trunc = trunctol(; atol = 1.0e-8) # test logging - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ, H, IDMRG2(; tol, verbosity = verbosity_full, maxiter = 2, trunc) ) - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ, H, IDMRG2(; tol, verbosity = verbosity_conv, trunc) ) v = variance(ψ, H, envs) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 + + @test info.numtrunc > 0 + @test info.ϵ_max > 0 + @test info.ϵ_max <= info.ϵ_total <= sqrt(info.numtrunc) * info.ϵ_max end # the gradient is computed concurrently over the unit cell, so the scheduler decides its @@ -277,7 +287,7 @@ end ψ₀ = repeat(InfiniteMPS(ℙ^2, ℙ^D), unit_cell_size) H = repeat(H_ref, unit_cell_size) - ψ′, envs, δ = with_scheduler(scheduler) do + ψ′, envs, info = with_scheduler(scheduler) do # test logging ψ₁, = find_groundstate( ψ₀, H, GradientGrassmann(; tol, verbosity = verbosity_full, maxiter = 2) @@ -289,7 +299,7 @@ end v = variance(ψ′, H, envs) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 end @@ -300,12 +310,12 @@ end alg = VUMPS(; tol = 100 * tol, verbosity = verbosity_conv, maxiter = 10) & GradientGrassmann(; tol, verbosity = verbosity_conv, maxiter = 50) - ψ, envs, δ = find_groundstate(ψ, H, alg) + ψ, envs, info = find_groundstate(ψ, H, alg) v = variance(ψ, H, envs) # test using low variance - @test sum(δ) ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 end @@ -335,13 +345,13 @@ end @testset "DMRG" begin # test logging passes - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H_lazy, DMRG(; tol, verbosity = verbosity_full, maxiter = 1) ) # compare states alg = DMRG(; tol, verbosity = verbosity_conv) - ψ, envs, δ = find_groundstate(ψ, H_lazy, alg) + ψ, envs, info = find_groundstate(ψ, H_lazy, alg) @test abs(dot(ψ₀, ψ)) ≈ 1 atol = atol end @@ -349,28 +359,28 @@ end @testset "DMRG2" begin # test logging passes trunc = truncrank(floor(Int, D * 1.5)) - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H_lazy, DMRG2(; tol, verbosity = verbosity_full, maxiter = 1, trunc) ) # compare states alg = DMRG2(; tol, verbosity = verbosity_conv, trunc) ψ, = find_groundstate(ψ₀, H, alg) - ψ_lazy, envs, δ = find_groundstate(ψ₀, H_lazy, alg) + ψ_lazy, envs, info = find_groundstate(ψ₀, H_lazy, alg) @test abs(dot(ψ₀, ψ_lazy)) ≈ 1 atol = atol end @testset "GradientGrassmann" begin # test logging passes - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H_lazy, GradientGrassmann(; tol, verbosity = verbosity_full, maxiter = 2) ) # compare states alg = GradientGrassmann(; tol, verbosity = verbosity_conv) ψ, = find_groundstate(ψ₀, H, alg) - ψ_lazy, envs, δ = find_groundstate(ψ₀, H_lazy, alg) + ψ_lazy, envs, info = find_groundstate(ψ₀, H_lazy, alg) @test abs(dot(ψ₀, ψ_lazy)) ≈ 1 atol = atol end @@ -397,26 +407,26 @@ end @testset "VUMPS" begin # test logging passes - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H_lazy, VUMPS(; tol, verbosity = verbosity_full, maxiter = 2) ) # compare states alg = VUMPS(; tol, verbosity = verbosity_conv) - ψ, envs, δ = find_groundstate(ψ, H_lazy, alg) + ψ, envs, info = find_groundstate(ψ, H_lazy, alg) @test abs(dot(ψ₀, ψ)) ≈ 1 atol = atol end @testset "IDMRG" begin # test logging passes - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H_lazy, IDMRG(; tol, verbosity = verbosity_full, maxiter = 2) ) # compare states alg = IDMRG(; tol, verbosity = verbosity_conv, maxiter = 300) - ψ, envs, δ = find_groundstate(ψ, H_lazy, alg) + ψ, envs, info = find_groundstate(ψ, H_lazy, alg) @test abs(dot(ψ₀, ψ)) ≈ 1 atol = atol end @@ -428,26 +438,26 @@ end trunc = truncrank(floor(Int, D * 1.5)) # test logging passes - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀′, H_lazy′, IDMRG2(; tol, verbosity = verbosity_full, maxiter = 2, trunc) ) # compare states alg = IDMRG2(; tol, verbosity = verbosity_conv, trunc) - ψ, envs, δ = find_groundstate(ψ, H_lazy′, alg) + ψ, envs, info = find_groundstate(ψ, H_lazy′, alg) @test abs(dot(ψ₀′, ψ)) ≈ 1 atol = atol end @testset "GradientGrassmann" begin # test logging passes - ψ, envs, δ = find_groundstate( + ψ, envs, info = find_groundstate( ψ₀, H_lazy, GradientGrassmann(; tol, verbosity = verbosity_full, maxiter = 2) ) # compare states alg = GradientGrassmann(; tol, verbosity = verbosity_conv) - ψ, envs, δ = find_groundstate(ψ₀, H_lazy, alg) + ψ, envs, info = find_groundstate(ψ₀, H_lazy, alg) @test abs(dot(ψ₀, ψ)) ≈ 1 atol = atol end diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index 39a3cd763..0350b8a61 100644 --- a/test/algorithms/timestep.jl +++ b/test/algorithms/timestep.jl @@ -272,24 +272,42 @@ end # fixed bond dimension sweep never discards anything, so the reported error is exactly zero @testset "no truncation" begin - @test last(timestep(ψ₀, H, 0.0, dt, TDVP())) == 0 - @test last(timestep(ψ₀, H, 0.0, dt, BUG())) == 0 + for alg in (TDVP(), BUG()) + info = last(timestep(ψ₀, H, 0.0, dt, alg)) + @test info isa MPSKit.AlgorithmInfo + @test info.ϵ_max == 0 + @test info.ϵ_total == 0 + @test info.numtrunc == 0 + @test info.numsteps == 1 + end end @testset "$(nameof(alg))" for alg in (TDVP2, BUG) - _, _, ϵ_loose = timestep(ψ₀, H, 0.0, dt, alg(; trunc = truncrank(32))) - ψ, _, ϵ_tight = timestep(ψ₀, H, 0.0, dt, alg(; trunc = truncrank(2))) - - @test ϵ_tight > 1.0e-3 # throw away real weight - @test ϵ_loose < ϵ_tight # throw away less weight with a more forgiving truncation - @test norm(ψ) < norm(ψ₀) # discarded weight = norm loss in real time with `normalize = false` + _, _, loose = timestep(ψ₀, H, 0.0, dt, alg(; trunc = truncrank(32))) + ψ, _, tight = timestep(ψ₀, H, 0.0, dt, alg(; trunc = truncrank(2))) + + @test tight.ϵ_total > 1.0e-3 # throw away real weight + @test tight.ϵ_max > 1.0e-4 + @test loose.ϵ_max < tight.ϵ_max # throw away less weight with a more forgiving truncation + @test tight.numtrunc > 0 + @test tight.ϵ_total >= tight.ϵ_max + @test tight.ϵ_total <= sqrt(tight.numtrunc) * tight.ϵ_max + # `ϵ_total` = norm loss in real time with `normalize = false` + @test norm(ψ)^2 ≈ norm(ψ₀)^2 - tight.ϵ_total^2 atol = 1.0e-12 end - @testset "accumulation over an evolution" begin + @testset "aggregation over an evolution" begin alg = TDVP2(; trunc = truncrank(2)) - _, _, ϵ_step = timestep(ψ₀, H, 0.0, dt, alg) - _, _, ϵ_total = time_evolve(ψ₀, H, 0:dt:(4 * dt), alg) - @test ϵ_total >= ϵ_step + nsteps = 4 + _, _, step = timestep(ψ₀, H, 0.0, dt, alg) + ψ, _, total = time_evolve(ψ₀, H, 0:dt:(nsteps * dt), alg) + + @test total.numsteps == nsteps + @test total.numtrunc >= step.numtrunc + @test total.ϵ_total >= step.ϵ_total + @test norm(ψ)^2 ≈ norm(ψ₀)^2 - total.ϵ_total^2 atol = 1.0e-12 + @test total.ϵ_max >= step.ϵ_max + @test total.ϵ_max <= total.ϵ_total end end diff --git a/test/misc/multifusion.jl b/test/misc/multifusion.jl index 3096062a2..b2df5e81a 100644 --- a/test/misc/multifusion.jl +++ b/test/misc/multifusion.jl @@ -77,16 +77,16 @@ module TestMultifusion V = Vect[I](M => 8) init = FiniteMPS(L, PD, V; left = Vect[I](M => 1), right = Vect[I](M => 1)) v₀ = variance(init, H) - ψ, envs, δ = find_groundstate(init, H, DMRG()) + ψ, envs, info = find_groundstate(init, H, DMRG()) v = variance(ψ, H) E = expectation_value(ψ, H, envs) - ψ2, envs2, δ2 = find_groundstate(init, H, DMRG2(; trunc = trunctol(; atol = 1.0e-6))) + ψ2, envs2, info2 = find_groundstate(init, H, DMRG2(; trunc = trunctol(; atol = 1.0e-6))) v2 = variance(ψ2, H) E2 = expectation_value(ψ2, H, envs2) - @test δ ≈ 0 atol = 1.0e-3 - @test δ2 ≈ 0 atol = 1.0e-3 + @test info.normres ≈ 0 atol = 1.0e-3 + @test info2.normres ≈ 0 atol = 1.0e-3 @test v < v₀ && v2 < v₀ @test isapprox(E, E2; atol = 1.0e-6) @@ -109,21 +109,21 @@ module TestMultifusion init = InfiniteMPS([PD, PD], [V, V]) v₀ = variance(init, H) tol = 1.0e-10 - ψ, envs, δ = find_groundstate(init, H, IDMRG(; tol = tol, maxiter = 400)) + ψ, envs, info = find_groundstate(init, H, IDMRG(; tol = tol, maxiter = 400)) E = expectation_value(ψ, H, envs) v = variance(ψ, H) - ψ2, envs2, δ2 = find_groundstate(init, H, IDMRG2(; tol = tol, trunc = trunctol(; atol = 1.0e-6), maxiter = 400)) + ψ2, envs2, info2 = find_groundstate(init, H, IDMRG2(; tol = tol, trunc = trunctol(; atol = 1.0e-6), maxiter = 400)) E2 = expectation_value(ψ2, H, envs2) v2 = variance(ψ2, H) - ψ3, envs3, δ3 = find_groundstate(init, H, VUMPS(; tol = tol, maxiter = 400)) + ψ3, envs3, info3 = find_groundstate(init, H, VUMPS(; tol = tol, maxiter = 400)) E3 = expectation_value(ψ3, H, envs3) v3 = variance(ψ3, H) @test isapprox(E, E2; atol = 1.0e-6) @test isapprox(E, E3; atol = 1.0e-6) - for delta in [δ, δ2, δ3] + for delta in [info.normres, info2.normres, info3.normres] @test delta ≈ 0 atol = 1.0e-3 end for var in [v, v2, v3] From 22390424a36c919cec3802a0e9c7846c06fdec2f Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Thu, 20 Aug 2026 08:28:23 +0200 Subject: [PATCH 19/29] rewrite docs addressing comments + explaining algorithminfo --- docs/src/man/algorithms.md | 180 +++++++++++++++++++++++++---------- src/utility/algorithminfo.jl | 13 ++- 2 files changed, 137 insertions(+), 56 deletions(-) diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 60a76ae2e..2543c5596 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -7,7 +7,7 @@ DocTestSetup = :(using MPSKit, TensorKit, MPSKitModels) Here is a collection of the algorithms that have been added to MPSKit.jl. If a particular algorithm is missing, feel free to let us know via an issue, or contribute via a PR. -## Groundstates +## Ground states One of the most prominent use-cases of MPS is to obtain the ground state of a given (quasi-) one-dimensional quantum Hamiltonian. In MPSKit.jl, this can be achieved through `find_groundstate`: @@ -16,7 +16,7 @@ In MPSKit.jl, this can be achieved through `find_groundstate`: find_groundstate ``` -The returned error measures convergence to a variational fixed point, which is not the same as accuracy; see [Ground-state accuracy](@ref). +The returned error measures convergence to a variational fixed point, which is not the same as accuracy; see [Ground state accuracy](@ref). There are a variety of algorithms that have been developed over the years, and many of them have been implemented in MPSKit. Keep in mind that some of them are exclusive to finite or infinite systems, while others may work for both. @@ -35,7 +35,7 @@ Here, we enumerate some of their properties in hopes of pointing you in the righ ### DMRG -Probably the most widely used algorithm for optimizing groundstates with MPS is [`DMRG`](@ref) and its variants. +Probably the most widely used algorithm for optimizing ground states with MPS is [`DMRG`](@ref) and its variants. This algorithm sweeps through the system, optimizing a single site or pair of sites while keeping all others fixed. Since this local problem can be solved efficiently, the global optimal state follows by alternating through the system. However, because of the single-site nature of this algorithm, this can never alter the bond dimension of the state, such that there is no way of dynamically increasing the precision. @@ -277,88 +277,159 @@ isapprox(Es[1] - E₀, 2(g - 1); rtol=1e-2) # infinite analytical result ## Errors and accuracy -Most algorithms in MPSKit report an error alongside their result, and the manual pages above refer to it as `ϵ` throughout. -That single name covers genuinely different quantities, and the differences matter. -This section clarifies the differences, and explains what they do and, just as important, what they don't measure. +The algorithms that solve for a state, particularly [`find_groundstate`](@ref), [`leading_boundary`](@ref), [`approximate`](@ref), [`timestep`](@ref) and [`time_evolve`](@ref), return an [`AlgorithmInfo`](@ref) as their last value, describing how they arrived at their result. +[`excitations`](@ref) and [`changebonds`](@ref) report nothing. +What limits their accuracy is covered below all the same. +A single bare number could not do this job, because the quantities involved are genuinely different. +A convergence measure and a truncation error answer different questions, are not comparable, and not every algorithm produces both. +The struct therefore names each field and leaves the ones an algorithm does not compute as `nothing`, rather than promising a meaning that isn't there. + +```@docs; canonical=false +AlgorithmInfo +``` + +The rest of this section explains what those quantities are, and - equally important - what they do not measure. ### The error convention -Where `ϵ` is a *truncation* error, it is the 2-norm of the discarded singular values of a single factorisation, so that `ϵ²` is the discarded weight and the squared norm of the factorised tensor drops by exactly `ϵ²`. -The name "discarded weight" refers to it truly representing a probability: the ``\sigma_\alpha^2`` are the eigenvalues of the reduced density matrix across the cut, i.e. the statistical weights of the Schmidt states, summing to 1 for a normalised state. -So `ϵ²` is the probability weight thrown away and `ϵ` is the corresponding amplitude, which is why squares appear wherever these errors are combined. -Note that `ϵ` is absolute rather than relative: it is not divided by the norm of the state, which under a non-renormalising algorithm drifts away from 1 precisely as truncation accumulates. +**In theory:** A truncation error measures how much a factorisation changed the tensor it acted on: replacing ``A`` by its rank-restricted approximation ``\tilde{A}`` costs + +```math +\epsilon = \lVert A - \tilde{A} \rVert . +``` -What differs between algorithms is how the per-bond values are *aggregated*. -For this reason, values between algorithms are not directly comparable. +Because a truncated SVD keeps the largest singular values, the discarded part is orthogonal to the kept part, and ``\epsilon`` is exactly the 2-norm of the discarded singular values ([Schollwöck](@cite schollwoeck2011)). + +**In practice:** That is precisely what MPSKit computes: every factorisation reports the 2-norm of what it discarded, and nothing more. In particular, the value is absolute, i.e. it is not normalised by ``\lVert A \rVert``. + +**Where the two are often conflated:** In DMRG the same quantity is usually called the "discarded weight" and is read as a probability: when ``A`` is the bond tensor of a normalised state, the ``\sigma_\alpha^2`` are eigenvalues of the reduced density matrix across the cut which sum to 1 for a normalised state, so ``\epsilon^2`` is the probability of the discarded subspace. +Even though viewing it as a probability is intuitively useful, it is not what makes the value well-defined, and it does not always apply here. +In particular, under non-renormalising algorithms the state norm drifts away from 1 precisely as truncation accumulates. +The probability intuition holds only when the factorised object is normalised, which is not always the case. +The definition ``\epsilon = \lVert A - \tilde{A} \rVert`` holds in every case, though. + +What differs between algorithms is how these per-factorisation values are summed up (*aggregated*) into the numbers they report, which is the subject of the next subsection. !!! warning - Not every `ϵ` is a truncation error. - The `ϵ` returned by [`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms is a *convergence* measure, with no truncation interpretation at all. - The two are unrelated quantities that happen to share a name; see the two sections below. + A convergence measure and a truncation error are unrelated quantities. + [`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms fill in the `normres` field, which has no truncation interpretation at all, while the truncating algorithms fill in `ϵ_max`/`ϵ_total`, which say nothing about convergence. + An algorithm that does both fills both, and they should not be compared with each other. + +#### Aggregating truncation errors + +This applies to every truncating algorithm. +A sweep of [`DMRG2`](@ref), [`IDMRG2`](@ref), [`TDVP2`](@ref), [`BUG`](@ref) or [`Zipup`](@ref) performs many local factorisations, each with its own ``\epsilon_k``, and there is no single number that answers every question one might ask of them. +Rather than pick one and hope the caller wants that one, [`AlgorithmInfo`](@ref) carries both aggregations under names that say what they are. + +`ϵ_max`, the largest single ``\epsilon_k``, is a worst case. +The point of it is that it is still a per-factorisation quantity: it is one of the ``\epsilon_k``, not a combination of them, so it does not grow just because the chain is longer or more sweeps or steps were taken. +That is what makes it the one to watch over the course of a run, or to compare between runs at different sizes. +An increase in `ϵ_max` means individual bonds are being cut harder, whereas an increase in `ϵ_total` may only mean there were more bonds to cut. +This is the same reason the ground state algorithms report a maximum over local gradient norms rather than a total. +It is also why `ϵ_max` is the one [`time_evolve`](@ref) logs; the ground state algorithms log their `normres` instead, since for them convergence rather than truncation is what the sweep is driving. -### Ground-state accuracy +It is also the field that a `trunc` setting most directly controls, though how directly depends on the strategy: -[`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms return the quantity their `tol` is compared against. -For the sweeping algorithms ([`DMRG`](@ref), [`DMRG2`](@ref), [`VUMPS`](@ref), [`IDMRG`](@ref), [`IDMRG2`](@ref)) this is the Galerkin error: the norm of the local gradient projected orthogonally to the current state. -It vanishes exactly at a variational fixed point. -[`GradientGrassmann`](@ref) instead reports the norm of the Riemannian gradient from its optimizer. -Both vanish at a fixed point and both are gradient norms. -However, since they are taken in different metrics, their magnitudes are not directly comparable. -In particular, a `tol` tuned for one is not a `tol` tuned for the other. +* [`truncerror`](@extref MatrixAlgebraKit.truncerror) bounds the discarded weight of each factorisation, which is exactly ``\epsilon_k``, so `ϵ_max` should come out at or below the tolerance you set. +* [`trunctol`](@extref MatrixAlgebraKit.trunctol) bounds each individual singular value instead. Discarding ``k`` of them leaves ``\epsilon_k \le \sqrt{k}\,\texttt{atol}``, so `ϵ_max` lands near the tolerance but is not bounded by it. +* [`truncrank`](@extref MatrixAlgebraKit.truncrank) fixes the rank and says nothing about magnitudes at all. Here, `ϵ_max` is not something you set but something you read off. It is thus the consequence of that choice of bond dimension. + +`ϵ_total`, all of them combined in quadrature, + +```math +\epsilon_{\text{total}} = \sqrt{\textstyle\sum_k \epsilon_k^2} , +``` + +is the one that adds up to something. +Note that it is the *squares* that are summed, and the root taken at the end, because ``\epsilon^2`` is the quantity that accumulates exactly: each local truncation is an orthogonal projection, so it removes exactly ``\epsilon_k^2`` from the squared norm of the tensor it acts on. +Summing the squares therefore tracks a running "cost" rather than merely a worst case. +The price is that it is extensive: it grows like ``\sqrt{N}`` in the number of truncations, so it is not comparable between different system sizes or sweep counts. + +Whether that running cost is also the error of the *final state* depends on what the algorithm does between truncations, so it is not a property of the aggregation itself. +It does hold for real-time evolution, which is worked out in [The norm as a record of truncation](@ref). +A variational sweep, by contrast, renormalises as it goes, so there `ϵ_total` is a diagnostic of how hard the truncation is working rather than a norm deficit. + +Neither field is the distance to the untruncated solution. +Bounding that gives the linear sum ``\lVert \psi_{\text{untruncated}} - \psi \rVert \le \sum_k \epsilon_k``, which is a third quantity again, and always the largest of the three. + +### Ground state accuracy + +[`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms report the quantity their `tol` is compared against as `normres`, together with a `converged` flag. +**In theory:** Convergence is measured by the (norm of the) variational gradient: the component of ``H \lvert \psi \rangle`` that points away from the current state but still lies in the tangent space of the variational manifold. +It vanishes exactly at a variational fixed point, and its norm is what both families of algorithms report; for the sweeping algorithms that norm is known as the Galerkin error. + +**In practice:** The sweeping algorithms ([`DMRG`](@ref), [`DMRG2`](@ref), [`VUMPS`](@ref), [`IDMRG`](@ref), [`IDMRG2`](@ref)) report the Galerkin error, computed per site as the norm of the local update projected onto the orthogonal complement of the current tensor, and then aggregated as the maximum over sites. +[`GradientGrassmann`](@ref) instead reports the gradient norm supplied by its optimiser, taken over the whole state at once in the preconditioned Grassmann metric. + +**Why they differ:** These are the same underlying gradient, not different physical quantities. +They differ in the metric it is measured in (the Grassmann gradient is preconditioned) and in how the per-site contributions are combined (a maximum versus a single global norm). +The maximum is a deliberate practical choice, as it keeps the reported number independent of system size. +The consequence of the mismatch is that a `tol` tuned for one is not a `tol` tuned for the other. In other words, convergence is only defined relative to the manifold you are optimising over. -A single-site algorithm at a fixed bond dimension can drive its `ϵ` to machine precision and still be far from the true ground state, because the error that remains is the bond dimension itself, which no amount of further sweeping can address. -A small `ϵ` certifies a fixed point, not an accurate state. +A single-site algorithm at a fixed bond dimension can drive its `normres` to machine precision and still be far from the true ground state, because the error that remains is the bond dimension itself, which no amount of further sweeping can address. +A small `normres` certifies a fixed point, not an accurate state. Growing the bond dimension is the job of the two-site algorithms ([`DMRG2`](@ref), [`IDMRG2`](@ref)) or of a bond expansion ([`DMRG`](@ref) with an `alg_expand`, or an expanding `alg_gauge` such as [`DMRG3S`](@ref)); see also [`changebonds`](@ref). Once an algorithm does truncate, the two error notions interact. -The Galerkin error cannot fall below the level set by the weight being discarded each sweep, so a truncating scheme converges once `ϵ` reaches the truncation error rather than the (unreachable) bare `tol`. +The Galerkin error cannot fall below the level set by the weight being discarded each sweep, so a truncating scheme converges once `normres` reaches the truncation error rather than the (unreachable) bare `tol`. [`DMRG`](@ref)/[`DMRG2`](@ref) account for this: their stopping test is `ϵ ≤ max(tol, maximum(ϵ_trunc))`, which reduces to the plain `ϵ ≤ tol` when nothing is truncated. -Neither measure is an error bar on an observable. -For that, the standard route is the energy variance ``\langle H^2 \rangle - \langle H \rangle^2`` and extrapolation of observables towards zero variance ([Hubig et al.](@cite hubig2018)). +Neither measure is an error bar on an observable, and no cheap substitute for one exists. +The energy variance ``\langle H^2 \rangle - \langle H \rangle^2`` is an independent and more demanding measure. +Note what it actually quantifies, namely how far the state is from being an *exact eigenstate*, which is not the same thing as the error on some other observable. +The usual practical route is to compute a series of states at increasing bond dimension and extrapolate observables against the variance towards zero ([Hubig et al.](@cite hubig2018)). +This is an empirical extrapolation rather than a bound. ### Time evolution accuracy -Unlike a ground-state search, a time evolution has no convergence criterion to run to. +Unlike a ground state search, a time evolution has no convergence criterion to run to. There is no fixed point, and the error is made at every step. Three sources behave differently and only one of them is reported. * **Truncation error.** Whenever a bond is cut back down, the discarded singular values are lost from the state. - [`timestep`](@ref) and [`time_evolve`](@ref) return this as their third value `ϵ`, the norm of the discarded component, so that `ϵ²` is the discarded weight. + [`timestep`](@ref) and [`time_evolve`](@ref) report this through the `ϵ_max`/`ϵ_total` fields of the [`AlgorithmInfo`](@ref) they return. It is the error you control through the algorithm's `trunc`, and the only one that is free to compute, since the truncating SVD produces it anyway. It is non-zero for [`TDVP2`](@ref), for [`BUG`](@ref) with a `trunc`, and for [`TDVP`](@ref) with a bond expansion. - Plain single-site [`TDVP`](@ref) runs at fixed bond dimension and returns exactly `0`. + Plain single-site [`TDVP`](@ref) runs at fixed bond dimension and reports exactly `0`. * **Projection error.** Single-site [`TDVP`](@ref) confines the evolution to the tangent space of a fixed-bond-dimension manifold, ``\lVert (1 - P_{T_\psi}) H \psi \rVert``. The component of the exact evolution pointing off that manifold is simply dropped, and this happens even with no truncation and exact local solves. - It is not reported as measuring it costs an extra effective-Hamiltonian application per site. + It is not reported, since measuring it costs an extra effective-Hamiltonian application per site. This is what a bond expansion (CBE) exists to reduce ([Li et al.](@cite li2024)). -* **Time-discretization error.** - The projector splitting is globally ``O(dt^2)`` for the symmetric back-and-forth sweep ([Lubich et al.](@cite lubich2015), [Paeckel et al.](@cite paeckel2019)), so it is controlled by `dt` alone. +* **Splitting error.** + The evolution is not applied in one piece. + Rather, it is split into local terms that are integrated in sequence, which do not commute. + This is a Trotter-type error. + Note, however, that what is split differs per method: [`TDVP`](@ref) splits the tangent-space projector into site and bond terms, while the MPO methods ([`WI`](@ref), [`WII`](@ref), [`TaylorCluster`](@ref)) split the Hamiltonian in the more familiar sense. + For TDVP's symmetric back-and-forth sweep the result is globally ``O(dt^2)`` ([Lubich et al.](@cite lubich2015), [Paeckel et al.](@cite paeckel2019)), so it is controlled by `dt` alone. This can only be estimated by comparing one step of `dt` against two of `dt / 2`. -A trustworthy run needs all three under control, not just a small `ϵ`. -In practice: pick `dt` from a convergence check, pick `trunc` from the accumulated `ϵ`, and use a bond-adaptive scheme ([`TDVP2`](@ref), [`BUG`](@ref), or [`TDVP`](@ref) with `alg_expand`) whenever entanglement grows during the evolution, since a fixed bond dimension silently converts entanglement growth into projection error. +A trustworthy run needs all three under control, not just a small truncation error. +In practice: pick `dt` from a convergence check, pick `trunc` from the reported truncation error, and use a bond-adaptive scheme ([`TDVP2`](@ref), [`BUG`](@ref), or [`TDVP`](@ref) with `alg_expand`) whenever entanglement grows during the evolution, since a fixed bond dimension silently converts entanglement growth into projection error. + +Note that these do not all shrink together, so there is a sweet spot in `dt` rather than a "smaller is always better". +Decreasing `dt` reduces the splitting error, but it also means more steps to reach the same final time, and every step truncates again. +The accumulated truncation error grows with the number of steps, as does the compute time. +The projection error does not improve at all, since it represents a rate at which the exact solution leaves the manifold. +Shrinking `dt` in this case only samples this rate more finely. +The practical consequence is that below some `dt` the total error stops improving and eventually gets worse, and the remedy at that point is a larger bond dimension rather than a smaller step. -#### Summing local errors in squares and relation to the norm +#### The norm as a record of truncation -``\epsilon^2`` is the quantity that adds exactly. -Each local truncation is an orthogonal projection, so it removes exactly ``\epsilon_k^2`` from the squared norm, and the substeps between truncations preserve the norm. -Summing the squares therefore tracks a conserved "cost". -The alternative is to sum ``\epsilon_k``: ``\lVert \psi_{\text{untruncated}} - \psi \rVert \le \sum_k \epsilon_k``, a distance to the untruncated solution, which is a different and always larger quantity. +How the per-factorisation errors are aggregated into `ϵ_max`/`ϵ_total` is described under [Aggregating truncation errors](@ref); what follows is specific to time evolution. By default none of the time evolution algorithms renormalize (`normalize = false`), which is deliberate. In real time the local exponentials are unitary, so truncation is the only thing that changes the norm and it becomes a running record of what truncation has cost, ```math -\lVert \psi \rVert^2 = \lVert \psi_0 \rVert^2 - \epsilon^2 . +\lVert \psi \rVert^2 = \lVert \psi_0 \rVert^2 - \epsilon_{\text{total}}^2 . ``` -The reported `ϵ` is the norm deficit, and this composes across steps. +The reported `ϵ_total` is the norm deficit, and this composes across steps. This follows from the following two facts put together, one per half of a local update. 1) An SVD truncation is an orthogonal projection onto the kept Schmidt vectors and is 2-norm optimal at that rank ([Schollwöck](@cite schollwoeck2011)), so the kept and discarded parts are orthogonal. By Pythagoras the squared norm drops by exactly the discarded weight, the usual way of quantifying truncation during a time evolution ([Paeckel et al.](@cite paeckel2019)). @@ -372,10 +443,12 @@ These two hold for the local updates of a step, so composing over all steps give It is also specific to real time with `normalize = false`: * **Imaginary time** evolves with the non-unitary ``\exp(-H dt)``, which rescales the state on its own. - The norm then moves for two independent reasons, namely the physical decay of the weight and the truncation loss. One cannot separate them from each other. `ϵ` still counts only the truncation. + The norm then moves for two independent reasons, namely the physical decay of the weight and the truncation loss. + One cannot separate them from each other. + The reported errors still count only the truncation. * **`normalize = true`** renormalizes at every local update, destroying the identity by construction. - This is usually what you want for imaginary-time evolution used as a ground-state or thermal-state search. - `ϵ` is still reported and is unaffected. + This is usually what you want for imaginary-time evolution used as a ground state or thermal-state search. + The reported errors are unaffected. * **No truncation at all** (plain single-site [`TDVP`](@ref), or [`BUG`](@ref) with a QR gauge) gives ``\epsilon = 0``, and in real time the norm is then conserved exactly. * An **`InfiniteMPS`** is regauged to norm 1 per site structurally, so its norm carries no such information and `normalize` has no effect. @@ -384,19 +457,22 @@ It is also specific to real time with `normalize = false`: [`excitations`](@ref) returns only `(energies, states)`: there is no error term, and none of the sources below is reported back to you. They are worth knowing about, because the dominant one is usually not the one the algorithm is working on. -* **Inherited ground-state error.** +* **Inherited ground state error.** Every method here builds on a ground state you supply and treats it as exact. Its error propagates straight into the gap, and since a gap is a difference of two large energies, it is typically the limiting factor. - A well-converged ground state (in the sense of `ϵ` and bond dimension) is necessary for a meaningful excitation calculation. + A well-converged ground state (in the sense of `normres` and bond dimension) is necessary for a meaningful excitation calculation. * **Ansatz limitation.** [`QuasiparticleAnsatz`](@ref) varies over the single-quasiparticle tangent space on top of a fixed ground state. It is variational within that space and well suited to isolated quasiparticle branches, but multi-particle continua are not representable in it, so results there are not to be trusted. - For infinite systems the momentum superposition itself is exact, so momentum is a good quantum number and no error enters through it. + How well it does on an isolated branch is controlled by the gaps around the targeted level. + The ansatz approximates an eigenvalue that is separated from the rest of the spectrum in its momentum sector with an error bounded exponentially in the size of the local operator's support, at a rate set by the gap below *and* above that eigenvalue ([Haegeman et al.](@cite haegeman2013)). + A branch that is nearly degenerate with the ground state, or that sits just below a continuum, therefore converges much more slowly than an isolated one. * **Eigensolver convergence.** - The local eigenvalue problem is solved with KrylovKit, and a run that fails to converge `num` states emits a warning on the residual when the verbosity is set high enough. - This residual is neither returned nor thrown, so it is worth not running with warnings suppressed. + The eigenvalue problem is solved with KrylovKit, and a run that fails to converge `num` states emits a warning carrying the residual when the verbosity is set high enough. + That residual is neither returned nor thrown, so it is worth not running with warnings suppressed. + Its convergence is governed by the same spectral structure as the ansatz above: levels that are well separated converge quickly, while nearly degenerate ones converge slowly and are the ones most likely to come back unconverged. * **Penalty-based orthogonality** ([`FiniteExcited`](@ref)). Higher states are found by minimising ``H + \lambda \sum_i |\psi_i\rangle\langle\psi_i|`` against the previously converged states, with ``\lambda`` the `weight` field. @@ -419,7 +495,7 @@ state. changebonds ``` -All of these are controlled by a `trunc`, and the weight they discard is measured the same way as explained in the `ϵ` convention under [The error convention](@ref). +All of these are controlled by a `trunc`, and the weight they discard is measured the same way as described under [The error convention](@ref). `changebonds` does not report it, since every algorithm has its own interpretation of the discarded singular values. There are several different algorithms implemented, each having their own advantages and diff --git a/src/utility/algorithminfo.jl b/src/utility/algorithminfo.jl index 1f511effd..f0473da8a 100644 --- a/src/utility/algorithminfo.jl +++ b/src/utility/algorithminfo.jl @@ -23,12 +23,17 @@ Both truncation fields are built from the same per-factorisation quantity, namel singular values a single local factorisation discarded, but aggregate it differently, because no single aggregation answers every question: - - `ϵ_max`: the largest of them. Independent of system size and iteration count, and on the same - scale as the per-bond threshold set by `trunc`, so this is what to compare against a tolerance - or between runs. It is also what gets logged. - - `ϵ_total`: all of them summed up in squares. This grows with system size and iteration count. + - `ϵ_max`: the largest of them. It is still a per-factorisation quantity rather than a combination + of them, so it does not grow with system size or iteration count, which is what makes it + comparable between runs. It is also the field a `trunc` setting most directly controls, though + how directly depends on the strategy. + - `ϵ_total`: all of them combined in quadrature, ``\\sqrt{\\sum_k \\epsilon_k^2}``. This grows + with system size and iteration count, so unlike `ϵ_max` it is not comparable between runs. - `numtrunc`: how many local factorisations actually discarded anything. +See [Aggregating truncation errors](@ref) for how the two relate to a `trunc` setting, and for the +per-strategy caveats. + These are `0` for an algorithm that never truncates, which does not mean the result is exact. Rather, it means this particular error channel is absent. See the manual on [Errors and accuracy](@ref) for what is *not* measured here. From e64b642818912e1f34b5b163a53b560119fbc909 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Thu, 20 Aug 2026 10:28:56 +0200 Subject: [PATCH 20/29] have examples use info --- examples/classic2d/1.hard-hexagon/main.jl | 2 +- examples/quantum1d/1.ising-cft/main.jl | 2 +- examples/quantum1d/2.haldane/main.jl | 6 +++--- examples/quantum1d/3.ising-dqpt/main.jl | 10 +++++----- examples/quantum1d/4.xxz-heisenberg/main.jl | 12 +++++++----- examples/quantum1d/5.haldane-spt/main.jl | 2 +- 6 files changed, 18 insertions(+), 16 deletions(-) diff --git a/examples/classic2d/1.hard-hexagon/main.jl b/examples/classic2d/1.hard-hexagon/main.jl index c7b3531b1..5ee8e0a5b 100644 --- a/examples/classic2d/1.hard-hexagon/main.jl +++ b/examples/classic2d/1.hard-hexagon/main.jl @@ -44,7 +44,7 @@ Additionally, we can compute the entanglement entropy as well as the correlation D = 10 V = virtual_space(D) ψ₀ = InfiniteMPS([P], [V]) -ψ, envs, = leading_boundary( +ψ, envs, info = leading_boundary( ψ₀, mpo, VUMPS(; verbosity = 0, alg_eigsolve = MPSKit.Defaults.alg_eigsolve(; ishermitian = false)) ) # use non-hermitian eigensolver diff --git a/examples/quantum1d/1.ising-cft/main.jl b/examples/quantum1d/1.ising-cft/main.jl index a3f45887e..c2434cfc1 100644 --- a/examples/quantum1d/1.ising-cft/main.jl +++ b/examples/quantum1d/1.ising-cft/main.jl @@ -122,7 +122,7 @@ can reach higher system sizes. L_mps = 20 H_mps = periodic_boundary_conditions(transverse_field_ising(), L_mps) D = 64 -ψ, envs, δ = find_groundstate(FiniteMPS(L_mps, ℂ^2, ℂ^D), H_mps, DMRG()); +ψ, envs, info = find_groundstate(FiniteMPS(L_mps, ℂ^2, ℂ^D), H_mps, DMRG()); md""" Excitations on top of the ground state can be found through the use of the quasiparticle diff --git a/examples/quantum1d/2.haldane/main.jl b/examples/quantum1d/2.haldane/main.jl index 98ac72a0e..4cb689e66 100644 --- a/examples/quantum1d/2.haldane/main.jl +++ b/examples/quantum1d/2.haldane/main.jl @@ -42,7 +42,7 @@ H = heisenberg_XXX(symmetry, chain; J, spin) physical_space = SU2Space(1 => 1) virtual_space = SU2Space(0 => 12, 1 => 12, 2 => 5, 3 => 3) ψ₀ = FiniteMPS(L, physical_space, virtual_space) -ψ, envs, delta = find_groundstate(ψ₀, H, DMRG(; verbosity = 0)) +ψ, envs, info = find_groundstate(ψ₀, H, DMRG(; verbosity = 0)) E₀ = real(expectation_value(ψ, H)) En_1, st_1 = excitations(H, QuasiparticleAnsatz(), ψ, envs; sector = SU2Irrep(1)) En_2, st_2 = excitations(H, QuasiparticleAnsatz(), ψ, envs; sector = SU2Irrep(2)) @@ -72,7 +72,7 @@ Ls = 12:4:30 @info "computing L = $L" ψ₀ = FiniteMPS(L, physical_space, virtual_space) H = heisenberg_XXX(symmetry, FiniteChain(L); J, spin) - ψ, envs, delta = find_groundstate(ψ₀, H, DMRG(; verbosity = 0)) + ψ, envs, info = find_groundstate(ψ₀, H, DMRG(; verbosity = 0)) En_1, st_1 = excitations(H, QuasiparticleAnsatz(), ψ, envs; sector = SU2Irrep(1)) En_2, st_2 = excitations(H, QuasiparticleAnsatz(), ψ, envs; sector = SU2Irrep(2)) return real(En_2[1] - En_1[1]) @@ -103,7 +103,7 @@ chain = InfiniteChain(1) H = heisenberg_XXX(symmetry, chain; J, spin) virtual_space_inf = Rep[SU₂](1 // 2 => 16, 3 // 2 => 16, 5 // 2 => 8, 7 // 2 => 4) ψ₀_inf = InfiniteMPS([physical_space], [virtual_space_inf]) -ψ_inf, envs_inf, delta_inf = find_groundstate(ψ₀_inf, H; verbosity = 0) +ψ_inf, envs_inf, info_inf = find_groundstate(ψ₀_inf, H; verbosity = 0) kspace = range(0, π, 16) Es, _ = excitations(H, QuasiparticleAnsatz(), kspace, ψ_inf, envs_inf; sector = SU2Irrep(1)) diff --git a/examples/quantum1d/3.ising-dqpt/main.jl b/examples/quantum1d/3.ising-dqpt/main.jl index 50f7a3a23..8c9219526 100644 --- a/examples/quantum1d/3.ising-dqpt/main.jl +++ b/examples/quantum1d/3.ising-dqpt/main.jl @@ -33,7 +33,7 @@ H₀ = transverse_field_ising(FiniteChain(L); g = -0.5) md""" ## Finite MPS quenching -We can define a helper function that measures the loschmith echo +We can define a helper function that measures the Loschmidt echo: """ echo(ψ₀::FiniteMPS, ψₜ::FiniteMPS) = -2 * log(abs(dot(ψ₀, ψₜ))) / length(ψ₀) @@ -46,7 +46,7 @@ We will initially use a two-site TDVP scheme to dynamically increase the bond di H₁ = transverse_field_ising(FiniteChain(L); g = -2.0) ψₜ = deepcopy(ψ₀) dt = 0.01 -ψₜ, envs = timestep(ψₜ, H₁, 0, dt, TDVP2(; trunc = truncrank(20))); +ψₜ, envs, info = timestep(ψₜ, H₁, 0, dt, TDVP2(; trunc = truncrank(20))); md""" "envs" is a kind of cache object that keeps track of all environments in `ψ`. It is often advantageous to re-use the environment, so that MPSKit doesn't need to recalculate everything. @@ -68,7 +68,7 @@ function finite_sim(L; dt = 0.05, finaltime = 5.0) for t in times[2:end] alg = t > 3 * dt ? TDVP() : TDVP2(; trunc = truncrank(50)) - ψₜ, envs = timestep(ψₜ, H₁, 0, dt, alg, envs) + ψₜ, envs, info = timestep(ψₜ, H₁, 0, dt, alg, envs) push!(echos, echo(ψₜ, ψ₀)) end @@ -114,7 +114,7 @@ H₁ = transverse_field_ising(; g = -2.0) # a single timestep is easy dt = 0.01 -ψₜ, envs = timestep(ψₜ, H₁, 0, dt, TDVP(), envs); +ψₜ, envs, info = timestep(ψₜ, H₁, 0, dt, TDVP(), envs); md""" With performance in mind we should once again try to re-use these "envs" cache objects. @@ -135,7 +135,7 @@ function infinite_sim(dt = 0.05, finaltime = 5.0) if t < 50dt # if t is sufficiently small, we increase the bond dimension ψₜ, envs = changebonds(ψₜ, H₁, OptimalExpand(; trunc = truncrank(1)), envs) end - ψₜ, envs = timestep(ψₜ, H₁, 0, dt, TDVP(), envs) + ψₜ, envs, info = timestep(ψₜ, H₁, 0, dt, TDVP(), envs) push!(echos, echo(ψₜ, ψ₀)) end diff --git a/examples/quantum1d/4.xxz-heisenberg/main.jl b/examples/quantum1d/4.xxz-heisenberg/main.jl index 36ce04d77..53d7c9e00 100644 --- a/examples/quantum1d/4.xxz-heisenberg/main.jl +++ b/examples/quantum1d/4.xxz-heisenberg/main.jl @@ -31,7 +31,8 @@ md""" The ground state can then be found by calling `find_groundstate`. """ -groundstate, cache, delta = find_groundstate(state, H, VUMPS()); +groundstate, cache, info = find_groundstate(state, H, VUMPS()); +info md""" As you can see, VUMPS struggles to converge. @@ -39,7 +40,8 @@ On its own, that is already quite curious. Maybe we can do better using another algorithm, such as gradient descent. """ -groundstate, cache, delta = find_groundstate(state, H, GradientGrassmann(; maxiter = 20)); +groundstate, cache, info = find_groundstate(state, H, GradientGrassmann(; maxiter = 20)); +info md""" Convergence is quite slow and even fails after sufficiently many iterations. @@ -71,7 +73,7 @@ Alternatively, the Hamiltonian can be constructed directly on a two-site unit ce ## H2 = repeat(H, 2); -- copies the one-site version H2 = heisenberg_XXX(ComplexF64, Trivial, InfiniteChain(2); spin = 1 // 2) -groundstate, envs, delta = find_groundstate( +groundstate, envs, info = find_groundstate( state, H2, VUMPS(; maxiter = 100, tol = 1.0e-12) ); @@ -80,7 +82,7 @@ We get convergence, but it takes an enormous amount of iterations. The reason behind this becomes more obvious at higher bond dimensions: """ -groundstate, envs, delta = find_groundstate( +groundstate, envs, info = find_groundstate( state, H2, IDMRG2(; trunc = truncrank(50), maxiter = 20, tol = 1.0e-12) ); entanglementplot(groundstate) @@ -125,4 +127,4 @@ Even though the bond dimension is higher than in the example without symmetry, c println(dim(V1)) println(dim(V2)) -groundstate, cache, delta = find_groundstate(state, H2, VUMPS(; maxiter = 400, tol = 1.0e-12)); +groundstate, cache, info = find_groundstate(state, H2, VUMPS(; maxiter = 400, tol = 1.0e-12)); diff --git a/examples/quantum1d/5.haldane-spt/main.jl b/examples/quantum1d/5.haldane-spt/main.jl index 937228e6d..2fd4b7294 100644 --- a/examples/quantum1d/5.haldane-spt/main.jl +++ b/examples/quantum1d/5.haldane-spt/main.jl @@ -80,7 +80,7 @@ non-injective MPS. ℋ = SU2Space(1 => 1) V_wrong = SU2Space(0 => 8, 1 // 2 => 8, 1 => 3, 3 // 2 => 3) ψ = InfiniteMPS(ℋ, V_wrong) -ψ, environments, δ = find_groundstate(ψ, H, VUMPS(; maxiter = 10)) +ψ, environments, info = find_groundstate(ψ, H, VUMPS(; maxiter = 10)) sectors = SU2Irrep[0, 1 // 2, 1, 3 // 2] transferplot(ψ; sectors, title = "Transfer matrix spectrum", legend = :outertop) From 808d257d4c03daa6d19ab1137d3851cbd67f4403 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Thu, 20 Aug 2026 17:12:16 +0200 Subject: [PATCH 21/29] scoping is very hard --- src/algorithms/approximate/idmrg.jl | 2 +- src/algorithms/statmech/idmrg.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl index 252df7828..f3b2a6d01 100644 --- a/src/algorithms/approximate/idmrg.jl +++ b/src/algorithms/approximate/idmrg.jl @@ -72,7 +72,7 @@ function approximate!( ϵ::Float64 = 2 * alg.tol log = IterLog("IDMRG2") O, ϕ = toapprox - local iter + local iter, acc LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ) diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 5bff1e922..cab17d0fb 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -69,7 +69,7 @@ function leading_boundary( size(ψ, 2) < 2 && throw(ArgumentError("unit cell should be >= 2")) ϵ::Float64 = 2 * alg.tol log = IterLog("IDMRG2") - local iter + local iter, acc LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ) From 57f2e9aee68017cb899e2830b40208cdfe098c07 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Thu, 20 Aug 2026 17:12:32 +0200 Subject: [PATCH 22/29] brain lag --- test/algorithms/timestep.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/algorithms/timestep.jl b/test/algorithms/timestep.jl index 0350b8a61..dd925a923 100644 --- a/test/algorithms/timestep.jl +++ b/test/algorithms/timestep.jl @@ -278,7 +278,7 @@ end @test info.ϵ_max == 0 @test info.ϵ_total == 0 @test info.numtrunc == 0 - @test info.numsteps == 1 + @test info.numiter == 1 end end @@ -302,7 +302,7 @@ end _, _, step = timestep(ψ₀, H, 0.0, dt, alg) ψ, _, total = time_evolve(ψ₀, H, 0:dt:(nsteps * dt), alg) - @test total.numsteps == nsteps + @test total.numiter == nsteps @test total.numtrunc >= step.numtrunc @test total.ϵ_total >= step.ϵ_total @test norm(ψ)^2 ≈ norm(ψ₀)^2 - total.ϵ_total^2 atol = 1.0e-12 From fb0287beeb10b713b2d785f381061359c3afc806 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Thu, 20 Aug 2026 17:23:07 +0200 Subject: [PATCH 23/29] update changelog --- docs/src/changelog.md | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 793027cb0..47b3bce3f 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -27,7 +27,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec a single sweep, optionally followed by a sweep in the opposite direction that imposes the final 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 `(ψ, ϵ)`. + `ψ` is a write target rather than an initial guess and may alias `ϕ`; they return `(ψ, info)`. - `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 @@ -41,16 +41,18 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Changed -- `timestep`/`timestep!`/`time_evolve`/`time_evolve!` now return `(ψ, envs, ϵ)` instead of - `(ψ, envs)`. `ϵ` is the truncation error, the norm of the singular values discarded by the local - gauge factorisations, summed in squares over the bonds of a step and over the iterations, so that `ϵ²` is the discarded/truncated weight. It is nonzero for `TDVP2`, for `BUG` with - a `trunc`, and for `TDVP` with a bond expansion. Plain single-site `TDVP` never truncates and - returns `0`. Existing `ψ, envs = timestep(...)` call sites keep working! `time_evolve` now - also logs the per-step error at `verbosity ≥ 3` and the total at `verbosity ≥ 2` (previously a - hardcoded `0`), under the algorithm's own name rather than always `"TDVP"`. +- The following algorithms now return an `AlgorithmInfo` in place of a bare error or nothing: `find_groundstate`, + `find_groundstate!`, `leading_boundary`, `approximate` and `approximate!` return + `(ψ, envs, info)` instead of `(ψ, envs, ϵ)`, `Zipup` returns `(ψ, info)`, and + `timestep`/`timestep!`/`time_evolve`/`time_evolve!` gain the same third value where they + previously returned none. This was motivated by the fact that a single number could not + carry what these algorithms actually produce. To migrate, replace `ϵ` with `info.normres` + for convergence measures and `info.ϵ_max` for truncation errors.([#512](https://github.com/QuantumKitHub/MPSKit.jl/pull/512)) - The meaning of every reported error and tolerance is now documented, and the manual has a new [Errors and accuracy](@ref) section covering ground states, time evolution and excitations - separately. The quantities themselves are unchanged, but what they represent is clarified. + separately. Each is written as what the quantity is in principle, what MPSKit actually computes, + and why the two differ where they do. Aside from the time-evolution return value, the + quantities themselves are unchanged.([#512](https://github.com/QuantumKitHub/MPSKit.jl/pull/512)) - 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, From d3419b44a3722c189f7be6a499180ba476113e12 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 21 Aug 2026 08:28:17 +0200 Subject: [PATCH 24/29] missed one --- test/algorithms/approximate.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/algorithms/approximate.jl b/test/algorithms/approximate.jl index c95175b49..b9b58e25f 100644 --- a/test/algorithms/approximate.jl +++ b/test/algorithms/approximate.jl @@ -140,11 +140,11 @@ end # both sweep directions, with and without the zip-down pass for left_to_right in (true, false), trunc′ in (trunc, (notrunc(), trunc)) - got_s, ϵ = approximate((expH, ψ), Zipup(; trunc = trunc′, left_to_right)) + got_s, info = approximate((expH, ψ), Zipup(; trunc = trunc′, left_to_right)) normalize!(got_s) @test norm(ref_s - got_s) < 0.002 @test norm(ψ - got_s) > 0.002 - @test ϵ < 1.0e-10 + @test info.ϵ_max < 1.0e-10 end end From 0ea04708d0c09783846fd47fc2ce41fa249b17dc Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Mon, 24 Aug 2026 12:09:50 +0200 Subject: [PATCH 25/29] fix sloppiness on reporting galerkin error where not the case [skip ci] --- docs/src/man/algorithms.md | 15 +++++++++++---- src/algorithms/approximate/approximate.jl | 10 +++++++--- src/algorithms/groundstate/idmrg.jl | 10 ++++++---- src/algorithms/statmech/leading_boundary.jl | 18 +++++++++++++----- 4 files changed, 37 insertions(+), 16 deletions(-) diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 2543c5596..306824e2c 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -356,15 +356,22 @@ Bounding that gives the linear sum ``\lVert \psi_{\text{untruncated}} - \psi \rV [`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms report the quantity their `tol` is compared against as `normres`, together with a `converged` flag. **In theory:** Convergence is measured by the (norm of the) variational gradient: the component of ``H \lvert \psi \rangle`` that points away from the current state but still lies in the tangent space of the variational manifold. -It vanishes exactly at a variational fixed point, and its norm is what both families of algorithms report; for the sweeping algorithms that norm is known as the Galerkin error. +It vanishes exactly at a variational fixed point, and its norm is what most of these algorithms report. +For the sweeping algorithms that norm is known as the Galerkin error. -**In practice:** The sweeping algorithms ([`DMRG`](@ref), [`DMRG2`](@ref), [`VUMPS`](@ref), [`IDMRG`](@ref), [`IDMRG2`](@ref)) report the Galerkin error, computed per site as the norm of the local update projected onto the orthogonal complement of the current tensor, and then aggregated as the maximum over sites. +**In practice:** The sweeping algorithms ([`DMRG`](@ref), [`DMRG2`](@ref), [`VUMPS`](@ref), [`VOMPS`](@ref)) report the Galerkin error, computed per site as the norm of the local update projected onto the orthogonal complement of the current tensor, and then aggregated as the maximum over sites. [`GradientGrassmann`](@ref) instead reports the gradient norm supplied by its optimiser, taken over the whole state at once in the preconditioned Grassmann metric. -**Why they differ:** These are the same underlying gradient, not different physical quantities. +[`IDMRG`](@ref) and [`IDMRG2`](@ref) report neither, and this is easy to miss because their `tol` sits alongside the others. +Their `normres` is a *fixed-point residual*: the change in the center bond tensor from one sweep to the next, ``\lVert C - C_{\text{prev}} \rVert``, with both tensors projected onto their common space when the bond space changed. +It measures how much a sweep still moves the state, which is weaker than measuring how far the state is from stationarity, since an algorithm crawling through a slow region reports a small value for the same reason a converged one does. +For [`IDMRG2`](@ref) the projection onto the common space also means that a change in bond dimension between sweeps is projected out of the measure rather than counted in it. + +**Why they differ:** The Galerkin and Grassmann measures are the same underlying gradient, not different physical quantities. They differ in the metric it is measured in (the Grassmann gradient is preconditioned) and in how the per-site contributions are combined (a maximum versus a single global norm). The maximum is a deliberate practical choice, as it keeps the reported number independent of system size. -The consequence of the mismatch is that a `tol` tuned for one is not a `tol` tuned for the other. +The IDMRG residual is not that gradient at all: it certifies that the sweeps have stopped moving, not that there is nowhere left to move to. +The consequence of these mismatches is that a `tol` tuned for one algorithm is not a `tol` tuned for another, and this is worth keeping in mind when swapping algorithms at a fixed `tol`. In other words, convergence is only defined relative to the manifold you are optimising over. A single-site algorithm at a fixed bond dimension can drive its `normres` to machine precision and still be far from the true ground state, because the error that remains is the bond dimension itself, which no amount of further sweeping can address. diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index 0115afe84..e0d9dfa0d 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -29,7 +29,8 @@ algorithm for you based on the type of `ψ₀` (`DMRG`/`DMRG2` for a finite MPS, `IDMRG2` for an infinite MPS) and only accepts the `(O, ψ)` tuple form of `toapprox`. Once you pass an explicit `algorithm`, keywords are no longer accepted here — configure the algorithm struct itself instead (e.g. `DMRG(; tol, maxiter, verbosity)`). -- `tol::Float64`: convergence tolerance, compared against the Galerkin error (see Returns below) +- `tol::Float64`: convergence tolerance, compared against `info.normres` (see Returns below). Which + quantity that is depends on the algorithm - `maxiter::Int`: maximum amount of iterations - `verbosity::Int`: display progress information - `trunc`: if supplied, a truncated two-site sweep (`DMRG2`/`IDMRG2`) is prepended to @@ -42,8 +43,11 @@ struct itself instead (e.g. `DMRG(; tol, maxiter, verbosity)`). - `info::AlgorithmInfo`: how the algorithm arrived there. Which of its fields are populated depends on the algorithm: - the iterative algorithms (`DMRG`, `DMRG2`, `IDMRG`, `IDMRG2`, `VOMPS`) fill `converged`, - `normres` and `numiter`. `normres` is the Galerkin error compared against `tol`, the same - quantity [`find_groundstate`](@ref) reports, measuring distance from the variational fixed point. + `normres` and `numiter`. `normres` is the quantity compared against `tol`, the same one + [`find_groundstate`](@ref) reports for that algorithm. For `DMRG`, `DMRG2` and `VOMPS` it is the + Galerkin error, measuring distance from the variational fixed point. For `IDMRG` and `IDMRG2` it + is instead a fixed-point residual, the change in the center bond tensor over a sweep, which says + the sweeps have stopped moving rather than that the state is variationally stationary. - the two-site ones (`DMRG2`, `IDMRG2`) which involve a truncated SVD fill the truncation fields with what their final sweep discarded, i.e. what the returned state is still throwing away per sweep rather than what the early, unconverged sweeps did. diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 8694f98e5..acb4176f3 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -12,8 +12,9 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_boundary`](@ref), and [`approximate`](@ref). """ @kwdef struct IDMRG{A, B} <: Algorithm - "convergence tolerance, compared against the Galerkin error (the tangent-space gradient - norm)" + "convergence tolerance, compared against the change in the center bond tensor over a sweep. + This is a fixed-point residual measuring how much a sweep still moves the state, and is not + equivalent to the Galerkin error that [`DMRG`](@ref) and [`VUMPS`](@ref) report" tol::Float64 = Defaults.tol "maximal amount of iterations" @@ -46,8 +47,9 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_boundary`](@ref), and [`approximate`](@ref). """ @kwdef struct IDMRG2{A, S, B} <: Algorithm - "convergence tolerance, compared against the Galerkin error (the tangent-space gradient - norm)" + "convergence tolerance, compared against the change in the center bond tensor over a sweep. + This is a fixed-point residual measuring how much a sweep still moves the state, and is not + equivalent to the Galerkin error that [`DMRG2`](@ref) reports" tol::Float64 = Defaults.tol "maximal amount of iterations" diff --git a/src/algorithms/statmech/leading_boundary.jl b/src/algorithms/statmech/leading_boundary.jl index 8124e4eb8..28061d67f 100644 --- a/src/algorithms/statmech/leading_boundary.jl +++ b/src/algorithms/statmech/leading_boundary.jl @@ -14,7 +14,8 @@ optimization algorithm will be attempted based on the supplied keywords. # Keyword Arguments -- `tol::Float64`: convergence tolerance, compared against the Galerkin error (see Returns below) +- `tol::Float64`: convergence tolerance, compared against `info.normres` (see Returns below). Which + quantity that is depends on the algorithm - `maxiter::Int`: maximum amount of iterations - `verbosity::Int`: display progress information @@ -22,10 +23,17 @@ optimization algorithm will be attempted based on the supplied keywords. - `ψ::AbstractMPS`: converged leading boundary MPS - `environments`: environments corresponding to the converged boundary -- `info::AlgorithmInfo`: how the algorithm terminated; `info.normres` is the Galerkin error compared - against `tol`, and `info.converged` whether it got there. It is not a truncation error. See - [`AlgorithmInfo`](@ref), [`find_groundstate`](@ref), and the manual on the `ϵ` convention under - [The error convention](@ref) and [Ground state accuracy](@ref). +- `info::AlgorithmInfo`: how the algorithm terminated; `info.normres` is the quantity compared + against `tol`, and `info.converged` whether it got there. It is not a truncation error. + Which quantity `normres` is depends on the algorithm: + - [`VUMPS`](@ref) and [`VOMPS`](@ref) report the Galerkin error, the maximum over sites of the + local update projected onto the orthogonal complement of the current tensor. + - [`GradientGrassmann`](@ref) reports the norm of the Riemannian gradient from its optimiser. + - [`IDMRG`](@ref) and [`IDMRG2`](@ref) report a fixed-point residual instead, the change in the + center bond tensor over a sweep. For `Multiline` methods this is extensive in the number of rows. + + See [`AlgorithmInfo`](@ref), [`find_groundstate`](@ref), and the manual on the `ϵ` convention + under [The error convention](@ref) and [Ground state accuracy](@ref). """ leading_boundary # TODO: alg selector From 4deddef092c957a6caec207391f6ec2719a618bf Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Mon, 24 Aug 2026 12:16:28 +0200 Subject: [PATCH 26/29] account properly for optimkit's returned history --- src/algorithms/groundstate/gradient_grassmann.jl | 4 ++-- src/algorithms/statmech/gradient_grassmann.jl | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/algorithms/groundstate/gradient_grassmann.jl b/src/algorithms/groundstate/gradient_grassmann.jl index 82ec02a63..70776b559 100644 --- a/src/algorithms/groundstate/gradient_grassmann.jl +++ b/src/algorithms/groundstate/gradient_grassmann.jl @@ -122,10 +122,10 @@ function find_groundstate( @infov 4 timeroutput end - normres = normgradhistory[end] + normres = normgradhistory[end, 2] # full history returned as [fhistory normgradhistory] info = AlgorithmInfo(; converged = normres <= alg.method.gradtol, normres, - numiter = size(normgradhistory, 1) + numiter = size(normgradhistory, 1) - 1 # history starts with initial point before first iteration ) return x, envs, info end diff --git a/src/algorithms/statmech/gradient_grassmann.jl b/src/algorithms/statmech/gradient_grassmann.jl index 3a07226bb..00ef1346d 100644 --- a/src/algorithms/statmech/gradient_grassmann.jl +++ b/src/algorithms/statmech/gradient_grassmann.jl @@ -19,10 +19,11 @@ function leading_boundary( alg.finalize!, isometrictransport = true ) - normres = normgradhistory[end] + + normres = normgradhistory[end, 2] # full history returned as [fhistory normgradhistory] info = AlgorithmInfo(; converged = normres <= alg.method.gradtol, normres, - numiter = size(normgradhistory, 1) + numiter = size(normgradhistory, 1) - 1 # history starts with initial point before first iteration ) return x, envs, info end From f396fbfbe22634f6d4f540946dbabc484df44a89 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 28 Aug 2026 10:40:46 +0200 Subject: [PATCH 27/29] update `AlgorithmInfo` to carry a `Dict` instead of fixed fields --- README.md | 2 +- docs/src/changelog.md | 6 +- docs/src/man/algorithms.md | 68 +++-- src/MPSKit.jl | 2 +- src/algorithms/approximate/approximate.jl | 24 +- src/algorithms/approximate/fvomps.jl | 4 +- src/algorithms/approximate/idmrg.jl | 4 +- src/algorithms/approximate/vomps.jl | 4 +- src/algorithms/approximate/zipup.jl | 6 +- src/algorithms/groundstate/dmrg.jl | 26 +- .../groundstate/find_groundstate.jl | 13 +- .../groundstate/gradient_grassmann.jl | 8 +- src/algorithms/groundstate/idmrg.jl | 19 +- src/algorithms/groundstate/vumps.jl | 6 +- src/algorithms/statmech/gradient_grassmann.jl | 2 +- src/algorithms/statmech/idmrg.jl | 4 +- src/algorithms/statmech/leading_boundary.jl | 21 +- src/algorithms/statmech/vomps.jl | 6 +- src/algorithms/timestep/bug.jl | 2 +- src/algorithms/timestep/tdvp.jl | 5 +- src/algorithms/timestep/time_evolve.jl | 26 +- src/utility/algorithminfo.jl | 278 +++++++++++++----- test/algorithms/approximate.jl | 2 +- test/algorithms/groundstate.jl | 24 +- test/misc/multifusion.jl | 6 +- 25 files changed, 383 insertions(+), 185 deletions(-) diff --git a/README.md b/README.md index 0543e08b8..4aeb0ded1 100644 --- a/README.md +++ b/README.md @@ -40,7 +40,7 @@ g_values = 0.1:0.1:2 M = @showprogress map(g_values) do g H = transverse_field_ising(; g=g) - groundstate, environment, δ = find_groundstate(init_state, H, VUMPS(; verbosity=0)) + groundstate, environment, info = find_groundstate(init_state, H, VUMPS(; verbosity=0)) return abs(expectation_value(groundstate, 1 => σᶻ())) end diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 7a63473d7..692bc0699 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -46,8 +46,10 @@ When releasing a new version, move the "Unreleased" changes to a new version sec `(ψ, envs, info)` instead of `(ψ, envs, ϵ)`, `Zipup` returns `(ψ, info)`, and `timestep`/`timestep!`/`time_evolve`/`time_evolve!` gain the same third value where they previously returned none. This was motivated by the fact that a single number could not - carry what these algorithms actually produce. To migrate, replace `ϵ` with `info.normres` - for convergence measures and `info.ϵ_max` for truncation errors.([#512](https://github.com/QuantumKitHub/MPSKit.jl/pull/512)) + carry what these algorithms actually produce. To migrate, replace `ϵ` with + `convergence_measure(info)` for convergence measures and `info.max_truncation_error` or `info.ϵ_max` + for truncation errors. See the updated docs or `AlgorithmInfo`'s docstring for more information. + ([#512](https://github.com/QuantumKitHub/MPSKit.jl/pull/512)) - The meaning of every reported error and tolerance is now documented, and the manual has a new [Errors and accuracy](@ref) section covering ground states, time evolution and excitations separately. Each is written as what the quantity is in principle, what MPSKit actually computes, diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 306824e2c..d4f05e7a3 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -282,7 +282,13 @@ The algorithms that solve for a state, particularly [`find_groundstate`](@ref), What limits their accuracy is covered below all the same. A single bare number could not do this job, because the quantities involved are genuinely different. A convergence measure and a truncation error answer different questions, are not comparable, and not every algorithm produces both. -The struct therefore names each field and leaves the ones an algorithm does not compute as `nothing`, rather than promising a meaning that isn't there. +Conversely, two algorithms that both report "a convergence measure" do not necessarily report the same measure. +[`AlgorithmInfo`](@ref) therefore carries a set of named entries, and each algorithm fills in only the ones it actually computes. +Entries are keyed by `Symbol`, and read either as a property (`info.max_truncation_error`) or by indexing (`info[:max_truncation_error]`). +The usual dictionary interface is supported, in particular `haskey`, `keys`, `values` and `get`. +Nothing is promised that was never measured: asking for an entry an algorithm did not report is an error naming what it did report, rather than a value that was never computed, and the name of the entry says which quantity it is rather than leaving that to prose. +Display the returned object (or call `keys(info)`) to see what a given algorithm reported. +The docstrings of every algorithm report these as well. ```@docs; canonical=false AlgorithmInfo @@ -312,7 +318,7 @@ What differs between algorithms is how these per-factorisation values are summed !!! warning A convergence measure and a truncation error are unrelated quantities. - [`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms fill in the `normres` field, which has no truncation interpretation at all, while the truncating algorithms fill in `ϵ_max`/`ϵ_total`, which say nothing about convergence. + [`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms fill in a convergence entry (`galerkin`, `gradientnorm`, `bondresidual` or `localchange` depending on the algorithm), which has no truncation interpretation at all, while the truncating algorithms fill in `max_truncation_error`/`total_truncation_error`, which say nothing about convergence. An algorithm that does both fills both, and they should not be compared with each other. #### Aggregating truncation errors @@ -321,20 +327,20 @@ This applies to every truncating algorithm. A sweep of [`DMRG2`](@ref), [`IDMRG2`](@ref), [`TDVP2`](@ref), [`BUG`](@ref) or [`Zipup`](@ref) performs many local factorisations, each with its own ``\epsilon_k``, and there is no single number that answers every question one might ask of them. Rather than pick one and hope the caller wants that one, [`AlgorithmInfo`](@ref) carries both aggregations under names that say what they are. -`ϵ_max`, the largest single ``\epsilon_k``, is a worst case. +`max_truncation_error`, the largest single ``\epsilon_k``, is a worst case. The point of it is that it is still a per-factorisation quantity: it is one of the ``\epsilon_k``, not a combination of them, so it does not grow just because the chain is longer or more sweeps or steps were taken. That is what makes it the one to watch over the course of a run, or to compare between runs at different sizes. -An increase in `ϵ_max` means individual bonds are being cut harder, whereas an increase in `ϵ_total` may only mean there were more bonds to cut. +An increase in `max_truncation_error` means individual bonds are being cut harder, whereas an increase in `total_truncation_error` may only mean there were more bonds to cut. This is the same reason the ground state algorithms report a maximum over local gradient norms rather than a total. -It is also why `ϵ_max` is the one [`time_evolve`](@ref) logs; the ground state algorithms log their `normres` instead, since for them convergence rather than truncation is what the sweep is driving. +It is also why `max_truncation_error` is the one [`time_evolve`](@ref) logs; the ground state algorithms log their convergence measure instead, since for them convergence rather than truncation is what the sweep is driving. It is also the field that a `trunc` setting most directly controls, though how directly depends on the strategy: -* [`truncerror`](@extref MatrixAlgebraKit.truncerror) bounds the discarded weight of each factorisation, which is exactly ``\epsilon_k``, so `ϵ_max` should come out at or below the tolerance you set. -* [`trunctol`](@extref MatrixAlgebraKit.trunctol) bounds each individual singular value instead. Discarding ``k`` of them leaves ``\epsilon_k \le \sqrt{k}\,\texttt{atol}``, so `ϵ_max` lands near the tolerance but is not bounded by it. -* [`truncrank`](@extref MatrixAlgebraKit.truncrank) fixes the rank and says nothing about magnitudes at all. Here, `ϵ_max` is not something you set but something you read off. It is thus the consequence of that choice of bond dimension. +* [`truncerror`](@extref MatrixAlgebraKit.truncerror) bounds the discarded weight of each factorisation, which is exactly ``\epsilon_k``, so `max_truncation_error` should come out at or below the tolerance you set. +* [`trunctol`](@extref MatrixAlgebraKit.trunctol) bounds each individual singular value instead. Discarding ``k`` of them leaves ``\epsilon_k \le \sqrt{k}\,\texttt{atol}``, so `max_truncation_error` lands near the tolerance but is not bounded by it. +* [`truncrank`](@extref MatrixAlgebraKit.truncrank) fixes the rank and says nothing about magnitudes at all. Here, `max_truncation_error` is not something you set but something you read off. It is thus the consequence of that choice of bond dimension. -`ϵ_total`, all of them combined in quadrature, +`total_truncation_error`, all of them combined in quadrature, ```math \epsilon_{\text{total}} = \sqrt{\textstyle\sum_k \epsilon_k^2} , @@ -347,14 +353,38 @@ The price is that it is extensive: it grows like ``\sqrt{N}`` in the number of t Whether that running cost is also the error of the *final state* depends on what the algorithm does between truncations, so it is not a property of the aggregation itself. It does hold for real-time evolution, which is worked out in [The norm as a record of truncation](@ref). -A variational sweep, by contrast, renormalises as it goes, so there `ϵ_total` is a diagnostic of how hard the truncation is working rather than a norm deficit. +A variational sweep, by contrast, renormalises as it goes, so there `total_truncation_error` is a diagnostic of how hard the truncation is working rather than a norm deficit. Neither field is the distance to the untruncated solution. Bounding that gives the linear sum ``\lVert \psi_{\text{untruncated}} - \psi \rVert \le \sum_k \epsilon_k``, which is a third quantity again, and always the largest of the three. +##### Which factorisations are recorded + +Both aggregations and the `numtrunc` count are built from the set of per-factorisation errors the algorithm records. +This set is not the same for every algorithm. + +[`IDMRG`](@ref), [`IDMRG2`](@ref), [`TDVP2`](@ref), [`BUG`](@ref) and [`Zipup`](@ref) record every factorisation as it happens. +Thus, `numtrunc` is a count of factorisations, and a sweep that visits a bond twice contributes two entries. + +[`DMRG`](@ref) and [`DMRG2`](@ref) keep one slot per update position instead, overwritten as the sweep passes over it, and record the slots once at the end. +In these cases, `numtrunc` is the number of positions whose most recent cut discarded something, and it can never exceed `length(psi)` (single-site) or `length(psi) - 1` (two-site), no matter how many sweeps ran. +A reported `numtrunc` less than the number of bonds means some bonds last discarded a non-zero weight. +These are typically the bonds nearest the two ends. + +The difference here is made deliberately. +What a returned variational state still throws away at a given bond is the last cut made there, not the sum of every cut ever made there. +This means that recording the latest value per position is the more meaningful quantity for a sweeping ground-state algorithm. +Thus, `numtrunc` counts different things in the two families, and `total_truncation_error` correspondingly sums a different number of terms, so neither is directly comparable between, say, [`DMRG2`](@ref) and [`IDMRG2`](@ref). +Read `numtrunc` more as "how many recorded errors were non-zero", and less as a tally of SVD calls. + +The two error entries also read under the short aliases `ϵ_max` and `ϵ_total`. +They are only ever stored under the descriptive names, so `keys` and displaying/showing them returns one name per quantity. + ### Ground state accuracy -[`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms report the quantity their `tol` is compared against as `normres`, together with a `converged` flag. +[`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms report the quantity their `tol` is compared against, together with a `converged` flag. +Because these are not the same quantity from one algorithm to the next, each is stored under a key that names it (`galerkin`, `gradientnorm`, `bondresidual` or `localchange`). +[`convergence_measure`](@ref) returns whichever of them is present, for code that only wants the number. **In theory:** Convergence is measured by the (norm of the) variational gradient: the component of ``H \lvert \psi \rangle`` that points away from the current state but still lies in the tangent space of the variational manifold. It vanishes exactly at a variational fixed point, and its norm is what most of these algorithms report. For the sweeping algorithms that norm is known as the Galerkin error. @@ -363,7 +393,7 @@ For the sweeping algorithms that norm is known as the Galerkin error. [`GradientGrassmann`](@ref) instead reports the gradient norm supplied by its optimiser, taken over the whole state at once in the preconditioned Grassmann metric. [`IDMRG`](@ref) and [`IDMRG2`](@ref) report neither, and this is easy to miss because their `tol` sits alongside the others. -Their `normres` is a *fixed-point residual*: the change in the center bond tensor from one sweep to the next, ``\lVert C - C_{\text{prev}} \rVert``, with both tensors projected onto their common space when the bond space changed. +Their `bondresidual` is a *fixed-point residual*: the change in the center bond tensor from one sweep to the next, ``\lVert C - C_{\text{prev}} \rVert``, with both tensors projected onto their common space when the bond space changed. It measures how much a sweep still moves the state, which is weaker than measuring how far the state is from stationarity, since an algorithm crawling through a slow region reports a small value for the same reason a converged one does. For [`IDMRG2`](@ref) the projection onto the common space also means that a change in bond dimension between sweeps is projected out of the measure rather than counted in it. @@ -374,12 +404,12 @@ The IDMRG residual is not that gradient at all: it certifies that the sweeps hav The consequence of these mismatches is that a `tol` tuned for one algorithm is not a `tol` tuned for another, and this is worth keeping in mind when swapping algorithms at a fixed `tol`. In other words, convergence is only defined relative to the manifold you are optimising over. -A single-site algorithm at a fixed bond dimension can drive its `normres` to machine precision and still be far from the true ground state, because the error that remains is the bond dimension itself, which no amount of further sweeping can address. -A small `normres` certifies a fixed point, not an accurate state. +A single-site algorithm at a fixed bond dimension can drive its convergence measure to machine precision and still be far from the true ground state, because the error that remains is the bond dimension itself, which no amount of further sweeping can address. +A small convergence measure certifies a fixed point, not an accurate state. Growing the bond dimension is the job of the two-site algorithms ([`DMRG2`](@ref), [`IDMRG2`](@ref)) or of a bond expansion ([`DMRG`](@ref) with an `alg_expand`, or an expanding `alg_gauge` such as [`DMRG3S`](@ref)); see also [`changebonds`](@ref). Once an algorithm does truncate, the two error notions interact. -The Galerkin error cannot fall below the level set by the weight being discarded each sweep, so a truncating scheme converges once `normres` reaches the truncation error rather than the (unreachable) bare `tol`. +The Galerkin error cannot fall below the level set by the weight being discarded each sweep, so a truncating scheme converges once `galerkin` reaches the truncation error rather than the (unreachable) bare `tol`. [`DMRG`](@ref)/[`DMRG2`](@ref) account for this: their stopping test is `ϵ ≤ max(tol, maximum(ϵ_trunc))`, which reduces to the plain `ϵ ≤ tol` when nothing is truncated. Neither measure is an error bar on an observable, and no cheap substitute for one exists. @@ -396,7 +426,7 @@ Three sources behave differently and only one of them is reported. * **Truncation error.** Whenever a bond is cut back down, the discarded singular values are lost from the state. - [`timestep`](@ref) and [`time_evolve`](@ref) report this through the `ϵ_max`/`ϵ_total` fields of the [`AlgorithmInfo`](@ref) they return. + [`timestep`](@ref) and [`time_evolve`](@ref) report this through the `max_truncation_error`/`total_truncation_error` entries of the [`AlgorithmInfo`](@ref) they return. It is the error you control through the algorithm's `trunc`, and the only one that is free to compute, since the truncating SVD produces it anyway. It is non-zero for [`TDVP2`](@ref), for [`BUG`](@ref) with a `trunc`, and for [`TDVP`](@ref) with a bond expansion. Plain single-site [`TDVP`](@ref) runs at fixed bond dimension and reports exactly `0`. @@ -427,7 +457,7 @@ The practical consequence is that below some `dt` the total error stops improvin #### The norm as a record of truncation -How the per-factorisation errors are aggregated into `ϵ_max`/`ϵ_total` is described under [Aggregating truncation errors](@ref); what follows is specific to time evolution. +How the per-factorisation errors are aggregated into `max_truncation_error`/`total_truncation_error` is described under [Aggregating truncation errors](@ref); what follows is specific to time evolution. By default none of the time evolution algorithms renormalize (`normalize = false`), which is deliberate. In real time the local exponentials are unitary, so truncation is the only thing that changes the norm and it becomes a running record of what truncation has cost, @@ -436,7 +466,7 @@ In real time the local exponentials are unitary, so truncation is the only thing \lVert \psi \rVert^2 = \lVert \psi_0 \rVert^2 - \epsilon_{\text{total}}^2 . ``` -The reported `ϵ_total` is the norm deficit, and this composes across steps. +The reported `total_truncation_error` is the norm deficit, and this composes across steps. This follows from the following two facts put together, one per half of a local update. 1) An SVD truncation is an orthogonal projection onto the kept Schmidt vectors and is 2-norm optimal at that rank ([Schollwöck](@cite schollwoeck2011)), so the kept and discarded parts are orthogonal. By Pythagoras the squared norm drops by exactly the discarded weight, the usual way of quantifying truncation during a time evolution ([Paeckel et al.](@cite paeckel2019)). @@ -467,7 +497,7 @@ They are worth knowing about, because the dominant one is usually not the one th * **Inherited ground state error.** Every method here builds on a ground state you supply and treats it as exact. Its error propagates straight into the gap, and since a gap is a difference of two large energies, it is typically the limiting factor. - A well-converged ground state (in the sense of `normres` and bond dimension) is necessary for a meaningful excitation calculation. + A well-converged ground state (in the sense of its convergence measure and bond dimension) is necessary for a meaningful excitation calculation. * **Ansatz limitation.** [`QuasiparticleAnsatz`](@ref) varies over the single-quasiparticle tangent space on top of a fixed ground state. diff --git a/src/MPSKit.jl b/src/MPSKit.jl index b386d7bce..8813e1fed 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 AlgorithmInfo +export AlgorithmInfo, convergence_measure export TDVP, TDVP2, BUG, WI, WII, TaylorCluster export changebonds, changebonds! export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index e0d9dfa0d..d788e6dff 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -29,8 +29,8 @@ algorithm for you based on the type of `ψ₀` (`DMRG`/`DMRG2` for a finite MPS, `IDMRG2` for an infinite MPS) and only accepts the `(O, ψ)` tuple form of `toapprox`. Once you pass an explicit `algorithm`, keywords are no longer accepted here — configure the algorithm struct itself instead (e.g. `DMRG(; tol, maxiter, verbosity)`). -- `tol::Float64`: convergence tolerance, compared against `info.normres` (see Returns below). Which - quantity that is depends on the algorithm +- `tol::Float64`: convergence tolerance, compared against the convergence entry of the returned + `info` (see Returns below). Which quantity that is depends on the algorithm - `maxiter::Int`: maximum amount of iterations - `verbosity::Int`: display progress information - `trunc`: if supplied, a truncated two-site sweep (`DMRG2`/`IDMRG2`) is prepended to @@ -42,17 +42,23 @@ struct itself instead (e.g. `DMRG(; tol, maxiter, verbosity)`). - `environments`: environments corresponding to the result (not returned by `Zipup`, which uses none) - `info::AlgorithmInfo`: how the algorithm arrived there. Which of its fields are populated depends on the algorithm: - - the iterative algorithms (`DMRG`, `DMRG2`, `IDMRG`, `IDMRG2`, `VOMPS`) fill `converged`, - `normres` and `numiter`. `normres` is the quantity compared against `tol`, the same one - [`find_groundstate`](@ref) reports for that algorithm. For `DMRG`, `DMRG2` and `VOMPS` it is the - Galerkin error, measuring distance from the variational fixed point. For `IDMRG` and `IDMRG2` it - is instead a fixed-point residual, the change in the center bond tensor over a sweep, which says - the sweeps have stopped moving rather than that the state is variationally stationary. + - the iterative algorithms (`DMRG`, `DMRG2`, `IDMRG`, `IDMRG2`, `VOMPS`) fill `converged` and + `numiter`, plus the quantity compared against `tol` under a key naming which measure it is: + - `VOMPS` reports `galerkin`, the Galerkin error, measuring distance from the variational + fixed point. + - `IDMRG` and `IDMRG2` report `bondresidual`, the change in the center bond tensor over a + sweep, which says the sweeps have stopped moving rather than that the state is stationary. + - `DMRG` and `DMRG2` report `localchange`, the largest relative change of a local tensor over + a sweep. Note this is not the Galerkin error they report in [`find_groundstate`](@ref). + + [`convergence_measure`](@ref) returns whichever of these is present, for code that only wants + the number. - the two-site ones (`DMRG2`, `IDMRG2`) which involve a truncated SVD fill the truncation fields with what their final sweep discarded, i.e. what the returned state is still throwing away per sweep rather than what the early, unconverged sweeps did. - [`Zipup`](@ref) is a single non-iterative sweep, so it has no convergence measure at all: - `converged` and `normres` are `nothing`, and it fills the truncation fields instead. + it reports no `converged` entry and none of the convergence entries, and fills the truncation + entries instead. See [`AlgorithmInfo`](@ref) for the full list and [The error convention](@ref) in the manual for why a convergence measure and a truncation error are not comparable quantities. diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl index 8dfa5cd19..5ca2a7045 100644 --- a/src/algorithms/approximate/fvomps.jl +++ b/src/algorithms/approximate/fvomps.jl @@ -34,7 +34,7 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environment end end - return ψ, envs, AlgorithmInfo(; converged = ϵ < alg.tol, normres = ϵ, numiter = iter) + return ψ, envs, AlgorithmInfo(; converged = ϵ < alg.tol, localchange = ϵ, numiter = iter) end function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments(ψ, _environment_args(Oϕ)...)) @@ -70,5 +70,5 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments end end - return ψ, envs, AlgorithmInfo(; converged = ϵ < alg.tol, normres = ϵ, numiter = iter) + return ψ, envs, AlgorithmInfo(; converged = ϵ < alg.tol, localchange = ϵ, numiter = iter) end diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl index f3b2a6d01..3f050b549 100644 --- a/src/algorithms/approximate/idmrg.jl +++ b/src/algorithms/approximate/idmrg.jl @@ -60,7 +60,7 @@ function approximate!( copy!(ψ, ψ′) # ensure output destination is unchanged recalculate!(envs, ψ, toapprox) - return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, normres = ϵ, numiter = iter) + return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, numiter = iter) end function approximate!( @@ -198,6 +198,6 @@ function approximate!( copy!(ψ, ψ′) # ensure output destination is unchanged recalculate!(envs, ψ, toapprox) - info = AlgorithmInfo(; converged = ϵ <= alg.tol, normres = ϵ, truncation = acc, numiter = iter) + info = AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, truncation = acc, numiter = iter) return ψ, envs, info end diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl index 22bdb675f..2acdbb00b 100644 --- a/src/algorithms/approximate/vomps.jl +++ b/src/algorithms/approximate/vomps.jl @@ -36,11 +36,11 @@ function _approximate_vomps(mps, toapprox, alg::VOMPS, envs) for (mps, envs, ϵ) in it if ϵ ≤ alg.tol @infov 2 logfinish!(log, it.iter, ϵ) - return mps, envs, AlgorithmInfo(; converged = true, normres = ϵ, numiter = it.iter) + return mps, envs, AlgorithmInfo(; converged = true, galerkin = ϵ, numiter = it.iter) end if it.iter ≥ alg.maxiter @warnv 1 logcancel!(log, it.iter, ϵ) - return mps, envs, AlgorithmInfo(; converged = false, normres = ϵ, numiter = it.iter) + return mps, envs, AlgorithmInfo(; converged = false, galerkin = ϵ, numiter = it.iter) end @infov 3 logiter!(log, it.iter, ϵ) end diff --git a/src/algorithms/approximate/zipup.jl b/src/algorithms/approximate/zipup.jl index e4c0b8eb5..6e5f4bc57 100644 --- a/src/algorithms/approximate/zipup.jl +++ b/src/algorithms/approximate/zipup.jl @@ -96,9 +96,9 @@ site with `alg_zipup`, and write the result into `ψ`. `zip_left_right!` zips up opposite direction imposes a final truncation with `alg_zipdown` in a locally gauged basis, leaving the gauge center of `ψ` at the far end. The destination may alias `ϕ`. -Also returns an [`AlgorithmInfo`](@ref) describing the truncation. Being a single sweep -rather than an iterative optimisation, there is no convergence measure, -so `converged` and `normres` are `nothing`. +Also returns an [`AlgorithmInfo`](@ref) describing the truncation. Being a single sweep +rather than an iterative optimisation, there is no convergence measure, so it reports neither +`converged` nor any convergence entry; [`convergence_measure`](@ref) returns `nothing` for it. """ zip_left_right! @doc (@doc zip_left_right!) zip_right_left! diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 105473d11..05543e32e 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -60,9 +60,9 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`approximate`](@ref). """ struct DMRG{A, F, E, G, B} <: Algorithm - "convergence tolerance on the Galerkin error (the tangent-space gradient norm). This acts as a - floor: the stopping test is `ϵ ≤ max(tol, maximum(ϵ_trunc))`, which reduces to `ϵ ≤ tol` - when nothing is truncated" + "convergence tolerance on the Galerkin error (the tangent-space gradient norm), reported as the + `galerkin` entry of the returned [`AlgorithmInfo`](@ref). This acts as a floor: the stopping + test is `ϵ ≤ max(tol, maximum(ϵ_trunc))`, which reduces to `ϵ ≤ tol` when nothing is truncated" tol::Float64 "maximal amount of iterations" @@ -173,9 +173,9 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`approximate`](@ref). """ struct DMRG2{A, G, F, B} <: Algorithm - "convergence tolerance on the Galerkin error (the tangent-space gradient norm). This acts as a - floor: the stopping test is `ϵ ≤ max(tol, maximum(ϵ_trunc))`, which reduces to `ϵ ≤ tol` - when nothing is truncated" + "convergence tolerance on the Galerkin error (the tangent-space gradient norm), reported as the + `galerkin` entry of the returned [`AlgorithmInfo`](@ref). This acts as a floor: the stopping + test is `ϵ ≤ max(tol, maximum(ϵ_trunc))`, which reduces to `ϵ ≤ tol` when nothing is truncated" tol::Float64 "maximal amount of iterations" @@ -277,10 +277,14 @@ Currently supported for the finite-system algorithms [`DMRG`](@ref) and [`DMRG2` - `ψ::AbstractFiniteMPS`: converged ground state - `environments`: environments corresponding to the converged state -- `info::AlgorithmInfo`: how the algorithm terminated. `info.normres` is the Galerkin error and - `info.converged` whether it met the stopping test. A truncating gauge also fills - `info.ϵ_max`/`info.ϵ_total` with what the final sweep discarded - (see [`find_groundstate`](@ref) and [`AlgorithmInfo`](@ref)) +- `info::AlgorithmInfo`: how the algorithm terminated. `info.galerkin` is the Galerkin error (also + reachable as [`convergence_measure`](@ref)) and `info.converged` whether it met the + stopping test. A truncating gauge also fills `info.max_truncation_error`/`info.total_truncation_error` + with what the final sweep discarded (see [`find_groundstate`](@ref) and [`AlgorithmInfo`](@ref)). + + These are built from one recorded value per update position, overwritten as the sweep passes + over it, so what is reported at each position is the most recent cut there rather than every + cut made during the sweep. See the manual on [Aggregating truncation errors](@ref). """ function find_groundstate!( ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs = environments(ψ, H, ψ) @@ -362,7 +366,7 @@ function _find_groundstate_sweep!( acc = TruncationAccumulator(Tr) foreach(ϵ -> push_error!(acc, ϵ), ϵ_truncs) info = AlgorithmInfo(; - converged = ϵ_global <= max(alg.tol, maximum(ϵ_truncs)), normres = ϵ_global, + converged = ϵ_global <= max(alg.tol, maximum(ϵ_truncs)), galerkin = ϵ_global, truncation = acc, numiter = iter ) return ψ, envs, info diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index 9a350d19b..67dc3e8ab 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/src/algorithms/groundstate/find_groundstate.jl @@ -38,11 +38,14 @@ low-bond-dimension initial guess such as a product state. - `ψ::AbstractMPS`: converged ground state - `environments`: environments corresponding to the converged state -- `info::AlgorithmInfo`: how the algorithm terminated. `info.normres` is the quantity compared - against `tol` and `info.converged` says whether it got there. Which measure `normres` is depends - on the algorithm. A truncating algorithm additionally fills `info.ϵ_max`/`info.ϵ_total` - with what its final sweep discarded. See [`AlgorithmInfo`](@ref), - and [The error convention](@ref) in the manual. +- `info::AlgorithmInfo`: how the algorithm terminated. `info.converged` says whether it met its + stopping criterion. The quantity that was compared against `tol` is stored under a key naming + which measure it is: `galerkin` for [`DMRG`](@ref), [`DMRG2`](@ref) and [`VUMPS`](@ref), + `gradientnorm` for [`GradientGrassmann`](@ref), and `bondresidual` for [`IDMRG`](@ref) and + [`IDMRG2`](@ref). [`convergence_measure`](@ref) returns whichever of these is present, + for code that only wants the number. A truncating algorithm additionally fills + `info.max_truncation_error`/`info.total_truncation_error` with what its final sweep discarded. + See [`AlgorithmInfo`](@ref) for the full vocabulary, and [The error convention](@ref) in the manual. # Examples diff --git a/src/algorithms/groundstate/gradient_grassmann.jl b/src/algorithms/groundstate/gradient_grassmann.jl index 70776b559..09f554444 100644 --- a/src/algorithms/groundstate/gradient_grassmann.jl +++ b/src/algorithms/groundstate/gradient_grassmann.jl @@ -15,8 +15,8 @@ with a preconditioner to induce the metric from the Hilbert space inner product. algorithm to construct - `finalize!`: finalizer algorithm - `tol = Defaults.tol`: convergence tolerance, compared against the norm of the Riemannian - (Grassmann) gradient reported by the optimizer. This is also the `ϵ` returned by - [`find_groundstate`](@ref) for this algorithm. + (Grassmann) gradient reported by the optimizer, which [`find_groundstate`](@ref) returns as the + `gradientnorm` entry of its [`AlgorithmInfo`](@ref). - `maxiter = Defaults.maxiter`: maximum amount of iterations - `verbosity = Defaults.verbosity - 1`: level of information display - `hasconverged = OptimKit.DefaultHasConverged(tol)`: convergence criterium @@ -81,7 +81,7 @@ end function find_groundstate( ψ::S, H, alg::GradientGrassmann, envs::P = environments(ψ, H, ψ) - )::Tuple{S, P, AlgorithmInfo{Float64}} where {S, P} + )::Tuple{S, P, AlgorithmInfo} where {S, P} !isa(ψ, FiniteMPS) || dim(ψ.C[end]) == 1 || @warn "This is not fully supported - split the mps up in a sum of mps's and optimize separately" normalize!(ψ) @@ -124,7 +124,7 @@ function find_groundstate( normres = normgradhistory[end, 2] # full history returned as [fhistory normgradhistory] info = AlgorithmInfo(; - converged = normres <= alg.method.gradtol, normres, + converged = normres <= alg.method.gradtol, gradientnorm = normres, numiter = size(normgradhistory, 1) - 1 # history starts with initial point before first iteration ) return x, envs, info diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index acb4176f3..82cd5d8e2 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -12,9 +12,10 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_boundary`](@ref), and [`approximate`](@ref). """ @kwdef struct IDMRG{A, B} <: Algorithm - "convergence tolerance, compared against the change in the center bond tensor over a sweep. - This is a fixed-point residual measuring how much a sweep still moves the state, and is not - equivalent to the Galerkin error that [`DMRG`](@ref) and [`VUMPS`](@ref) report" + "convergence tolerance, compared against the change in the center bond tensor over a sweep, + reported as the `bondresidual` entry of the returned [`AlgorithmInfo`](@ref). This is a + fixed-point residual measuring how much a sweep still moves the state, and is not equivalent to + the Galerkin error that [`DMRG`](@ref) and [`VUMPS`](@ref) report" tol::Float64 = Defaults.tol "maximal amount of iterations" @@ -47,9 +48,10 @@ $(TYPEDFIELDS) Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_boundary`](@ref), and [`approximate`](@ref). """ @kwdef struct IDMRG2{A, S, B} <: Algorithm - "convergence tolerance, compared against the change in the center bond tensor over a sweep. - This is a fixed-point residual measuring how much a sweep still moves the state, and is not - equivalent to the Galerkin error that [`DMRG2`](@ref) reports" + "convergence tolerance, compared against the change in the center bond tensor over a sweep, + reported as the `bondresidual` entry of the returned [`AlgorithmInfo`](@ref). This is a + fixed-point residual measuring how much a sweep still moves the state, and is not equivalent to + the Galerkin error that [`DMRG2`](@ref) reports" tol::Float64 = Defaults.tol "maximal amount of iterations" @@ -143,8 +145,9 @@ function _find_groundstate_idmrg(mps, operator, alg::alg_type, envs) where {alg_ ψ′ = InfiniteMPS(it.state.mps.AR; alg_gauge.tol, alg_gauge.maxiter) envs = recalculate!(it.state.envs, ψ′, it.state.operator, ψ′) info = AlgorithmInfo(; - converged = it.state.ϵ <= alg.tol, normres = it.state.ϵ, - truncation = it.state.truncation, numiter = it.state.iter + converged = it.state.ϵ <= alg.tol, bondresidual = it.state.ϵ, + truncation = alg isa IDMRG2 ? it.state.truncation : nothing, + numiter = it.state.iter ) return ψ′, envs, info end diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index 5ac67deb5..fe5481a2b 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -18,7 +18,7 @@ Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`leading_bou """ @kwdef struct VUMPS{F, B} <: Algorithm "convergence tolerance, compared against the Galerkin error (the tangent-space gradient - norm)" + norm), reported as the `galerkin` entry of the returned [`AlgorithmInfo`](@ref)" tol::Float64 = Defaults.tol "maximal amount of iterations" @@ -84,12 +84,12 @@ function dominant_eigsolve( if ϵ ≤ alg.tol @infov 4 timeroutput @infov 2 logfinish!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) - return mps, envs, AlgorithmInfo(; converged = true, normres = ϵ, numiter = it.iter) + return mps, envs, AlgorithmInfo(; converged = true, galerkin = ϵ, numiter = it.iter) end if it.iter ≥ alg.maxiter @infov 4 timeroutput @warnv 1 logcancel!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) - return mps, envs, AlgorithmInfo(; converged = false, normres = ϵ, numiter = it.iter) + return mps, envs, AlgorithmInfo(; converged = false, galerkin = ϵ, numiter = it.iter) end @infov 3 logiter!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) end diff --git a/src/algorithms/statmech/gradient_grassmann.jl b/src/algorithms/statmech/gradient_grassmann.jl index 00ef1346d..67e1422d7 100644 --- a/src/algorithms/statmech/gradient_grassmann.jl +++ b/src/algorithms/statmech/gradient_grassmann.jl @@ -22,7 +22,7 @@ function leading_boundary( normres = normgradhistory[end, 2] # full history returned as [fhistory normgradhistory] info = AlgorithmInfo(; - converged = normres <= alg.method.gradtol, normres, + converged = normres <= alg.method.gradtol, gradientnorm = normres, numiter = size(normgradhistory, 1) - 1 # history starts with initial point before first iteration ) return x, envs, info diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 000fcb81a..6bcc920a7 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -59,7 +59,7 @@ function leading_boundary( ψ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter) recalculate!(envs, ψ, operator, ψ) - return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, normres = ϵ, numiter = iter) + return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, numiter = iter) end function leading_boundary( @@ -201,5 +201,5 @@ function leading_boundary( ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter) recalculate!(envs, ψ, operator, ψ) - return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, normres = ϵ, truncation = acc, numiter = iter) + return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, truncation = acc, numiter = iter) end diff --git a/src/algorithms/statmech/leading_boundary.jl b/src/algorithms/statmech/leading_boundary.jl index 28061d67f..6960bfa3f 100644 --- a/src/algorithms/statmech/leading_boundary.jl +++ b/src/algorithms/statmech/leading_boundary.jl @@ -14,8 +14,8 @@ optimization algorithm will be attempted based on the supplied keywords. # Keyword Arguments -- `tol::Float64`: convergence tolerance, compared against `info.normres` (see Returns below). Which - quantity that is depends on the algorithm +- `tol::Float64`: convergence tolerance, compared against the convergence entry of the returned + `info` (see Returns below). Which quantity that is depends on the algorithm - `maxiter::Int`: maximum amount of iterations - `verbosity::Int`: display progress information @@ -23,14 +23,17 @@ optimization algorithm will be attempted based on the supplied keywords. - `ψ::AbstractMPS`: converged leading boundary MPS - `environments`: environments corresponding to the converged boundary -- `info::AlgorithmInfo`: how the algorithm terminated; `info.normres` is the quantity compared - against `tol`, and `info.converged` whether it got there. It is not a truncation error. - Which quantity `normres` is depends on the algorithm: - - [`VUMPS`](@ref) and [`VOMPS`](@ref) report the Galerkin error, the maximum over sites of the +- `info::AlgorithmInfo`: how the algorithm terminated; `info.converged` says whether it got there. + The quantity compared against `tol` is stored under a key naming which measure it is, and is + never a truncation error: + - [`VUMPS`](@ref) and [`VOMPS`](@ref) report `galerkin`, the maximum over sites of the local update projected onto the orthogonal complement of the current tensor. - - [`GradientGrassmann`](@ref) reports the norm of the Riemannian gradient from its optimiser. - - [`IDMRG`](@ref) and [`IDMRG2`](@ref) report a fixed-point residual instead, the change in the - center bond tensor over a sweep. For `Multiline` methods this is extensive in the number of rows. + - [`GradientGrassmann`](@ref) reports `gradientnorm` from its optimiser. + - [`IDMRG`](@ref) and [`IDMRG2`](@ref) report `bondresidual`, the change in the center bond + tensor over a sweep. For `Multiline` methods this is extensive in the number of rows. + + [`convergence_measure`](@ref) returns whichever is present, for code that only wants the + number. See [`AlgorithmInfo`](@ref), [`find_groundstate`](@ref), and the manual on the `ϵ` convention under [The error convention](@ref) and [Ground state accuracy](@ref). diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl index 5733abed1..626d17b05 100644 --- a/src/algorithms/statmech/vomps.jl +++ b/src/algorithms/statmech/vomps.jl @@ -19,7 +19,7 @@ Used as the `algorithm` argument of [`leading_boundary`](@ref) and [`approximate """ @kwdef struct VOMPS{F, B} <: Algorithm "convergence tolerance, compared against the Galerkin error (the tangent-space gradient - norm)" + norm), reported as the `galerkin` entry of the returned [`AlgorithmInfo`](@ref)" tol::Float64 = Defaults.tol "maximal amount of iterations" @@ -74,11 +74,11 @@ function dominant_eigsolve( for (mps, envs, ϵ) in it if ϵ ≤ alg.tol @infov 2 logfinish!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) - return mps, envs, AlgorithmInfo(; converged = true, normres = ϵ, numiter = it.iter) + return mps, envs, AlgorithmInfo(; converged = true, galerkin = ϵ, numiter = it.iter) end if it.iter ≥ alg.maxiter @warnv 1 logcancel!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) - return mps, envs, AlgorithmInfo(; converged = false, normres = ϵ, numiter = it.iter) + return mps, envs, AlgorithmInfo(; converged = false, galerkin = ϵ, numiter = it.iter) end @infov 3 logiter!(log, it.iter, ϵ, expectation_value(mps, operator, envs)) end diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl index 5c3ecc39a..c15320b4e 100644 --- a/src/algorithms/timestep/bug.jl +++ b/src/algorithms/timestep/bug.jl @@ -26,7 +26,7 @@ To restore a maximal dimension of `D`, apply [`changebonds`](@ref) with an [`Svd By default the state is not renormalized, as the (loss of) norm accumulates useful information. In real time the squared norm lost is exactly the weight discarded by the bond cuts, ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon_{\\text{total}}^2``, with - `ϵ_total` from the [`AlgorithmInfo`](@ref) returned by [`timestep`](@ref). In imaginary time + `total_truncation_error` from the [`AlgorithmInfo`](@ref) returned by [`timestep`](@ref). In imaginary time the norm also carries the physical decay of the weight and no longer isolates the truncation. Both reported errors are exactly zero when not truncating. Pass `normalize = true` to `timestep`/`time_evolve` to renormalize after every half-sweep instead. diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 79c5d8681..83ccbdb1f 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -15,7 +15,7 @@ state-preserving, as required for a consistent time evolution. so the state norm keeps useful information. In real time this is exact, namely the squared norm drops by precisely the truncated ("discarded") weight, ``\\lVert \\psi \\rVert^2 = \\lVert \\psi_0 \\rVert^2 - \\epsilon_{\\text{total}}^2``, - with `ϵ_total` from the [`AlgorithmInfo`](@ref) returned by [`timestep`](@ref). In imaginary + with `total_truncation_error` from the [`AlgorithmInfo`](@ref) returned by [`timestep`](@ref). In imaginary time the norm also carries the physical decay of the weight, so it no longer isolates the truncation. Without `trunc` nothing is discarded at all and the norm is conserved exactly in real time. @@ -135,7 +135,8 @@ function _timestep_infinite( end recalculate!(envs, ψ′, H) - # infinite one-site TDVP has a fixed bond dimension and never truncates so nothing is discarded + # infinite one-site TDVP runs at fixed bond dimension and never truncates, so it doesn't + # report truncation entries (rather than report zeros that would look like a measurement) # the gauge-fixing residual is controlled by `tolgauge`, not reported here return ψ′, envs, AlgorithmInfo() end diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 4c4bd0b68..6483952dd 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -26,12 +26,14 @@ through each of the time points obtained by iterating t_span. - `ψ`: the time-stepped state - `envs`: the updated environment manager -- `info::AlgorithmInfo`: the truncation performed over the whole evolution. +- `info::AlgorithmInfo`: the truncation performed over the whole evolution, accumulated from the + individual steps, together with `numiter`, the number of steps taken. An algorithm that never + truncates reports no truncation entries at all. See [`AlgorithmInfo`](@ref) and [Time evolution accuracy](@ref) in the manual for the difference and when to use which reported error measure, and [`timestep`](@ref) for what neither measures. -`ϵ_max` is logged per step at `verbosity ≥ 3` and for the whole evolution at `verbosity ≥ 2`. +`max_truncation_error` is logged per step at `verbosity ≥ 3` and for the whole evolution at `verbosity ≥ 2`. The size-independent measure is used here for the same reason as the ground state algorithms. """ function time_evolve end, function time_evolve! end @@ -43,7 +45,7 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time verbosity::Int = 0, imaginary_evolution::Bool = false, normalize::Bool = false ) log = IterLog(string(nameof(typeof(alg)))) - info = AlgorithmInfo(; truncation = TruncationAccumulator(ψ), numiter = 0) + info = AlgorithmInfo(; numiter = 0) LoggingExtras.withlevel(; verbosity) do @infov 2 loginit!(log, 0.0, first(t_span)) for iter in 1:(length(t_span) - 1) @@ -57,9 +59,10 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time info = _combine(info, info_step) # log the size-independent error measure - @infov 3 logiter!(log, iter, convert(Float64, info_step.ϵ_max), t) + # for a non-truncating algorithm, log zero + @infov 3 logiter!(log, iter, convert(Float64, get(info_step, :max_truncation_error, 0.0)), t) end - @infov 2 logfinish!(log, length(t_span), convert(Float64, info.ϵ_max), t_span[end]) + @infov 2 logfinish!(log, length(t_span), convert(Float64, get(info, :max_truncation_error, 0.0)), t_span[end]) end return ψ, envs, info end @@ -99,15 +102,16 @@ solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``. A step performs many local factorisations, each discarding some weight. Rather than collapse those into one number, `info` reports both aggregations under names that say what they are: -`info.ϵ_max` is the largest single one (size-independent, comparable against `trunc` and across -runs) and `info.ϵ_total` sums them in squares. +`info.max_truncation_error` is the largest single one (size-independent, comparable against `trunc` and across +runs) and `info.total_truncation_error` sums them in squares. Both are non-zero only for algorithms that truncate ([`TDVP2`](@ref), [`BUG`](@ref) with a -`trunc`, and [`TDVP`](@ref) with a bond expansion), and are exactly `0` for one-site -[`TDVP`](@ref), which runs at fixed bond dimension. A zero here does not mean the step was exact, -but that this particular error channel is absent. +`trunc`, and [`TDVP`](@ref) with a bond expansion). A finite-system step that truncates but +happened to discard nothing reports them as exactly `0`, whereas infinite one-site [`TDVP`](@ref) +never truncates and reports no truncation entries at all. Neither case means the step was exact, +but rather that this particular source of error is either absent or idle. -See [`AlgorithmInfo`](@ref) for the fields, and [Time evolution accuracy](@ref) in the manual +See [`AlgorithmInfo`](@ref) for the entries, and [Time evolution accuracy](@ref) in the manual for the other error sources. # Examples diff --git a/src/utility/algorithminfo.jl b/src/utility/algorithminfo.jl index f0473da8a..7021768cf 100644 --- a/src/utility/algorithminfo.jl +++ b/src/utility/algorithminfo.jl @@ -5,70 +5,179 @@ Information about how an algorithm arrived at its result, returned as the last v [`find_groundstate`](@ref), [`leading_boundary`](@ref), [`approximate`](@ref), [`timestep`](@ref) and [`time_evolve`](@ref). -Algorithms in MPSKit produce genuinely different error measures, and a single bare number cannot -carry that distinction. This struct keeps the return signature uniform while naming each quantity, -so that nothing is promised that an algorithm does not actually compute. +Algorithms in MPSKit produce genuinely different measures, and not every algorithm even has access +to the same information. To avoid reporting two different quantities under one name, the information is +carried in a `Dict{Symbol, Any}` that each algorithm fills with only the entries it actually +computes. -## Convergence +Entries are read as properties (`info.galerkin`), by indexing (`info[:galerkin]`), or through the +usual dictionary interface (`keys`, `haskey`, `get`, `pairs`, `length`). Asking for an entry the +algorithm never reported is an error that names what it did report, rather than a silent +`nothing`. Displaying the object (or calling `keys(info)`) shows what a given algorithm actually produced. - - `converged`: whether the algorithm reached its stopping criterion, or `nothing` for algorithms - that do not iterate to a fixed point. - - `normres`: the quantity compared against the algorithm's `tol`, or `nothing` when there is none. - Which measure this is depends on the algorithm. - - `numiter`: number of iterations. +## The vocabulary -## Truncation +The keys below are the ones currently in use. Each algorithm's own docstring states which of them +it reports. Nothing prevents an algorithm from adding its own. -Both truncation fields are built from the same per-factorisation quantity, namely the 2-norm of the -singular values a single local factorisation discarded, but aggregate it differently, because no -single aggregation answers every question: +### Convergence - - `ϵ_max`: the largest of them. It is still a per-factorisation quantity rather than a combination - of them, so it does not grow with system size or iteration count, which is what makes it - comparable between runs. It is also the field a `trunc` setting most directly controls, though - how directly depends on the strategy. - - `ϵ_total`: all of them combined in quadrature, ``\\sqrt{\\sum_k \\epsilon_k^2}``. This grows - with system size and iteration count, so unlike `ϵ_max` it is not comparable between runs. - - `numtrunc`: how many local factorisations actually discarded anything. + - `converged::Bool`: whether the algorithm met its stopping criterion. + - `numiter::Int`: number of iterations (sweeps or steps). + +The quantity that was compared against the algorithm's `tol` is stored under a name that says +which measure it is: + + - `galerkin`: the Galerkin error, i.e. the maximum over sites of the local update projected onto + the orthogonal complement of the current tensor. Reported by [`DMRG`](@ref), [`DMRG2`](@ref), + [`VUMPS`](@ref) and [`VOMPS`](@ref) when solving for a state. + - `gradientnorm`: the norm of the Riemannian (Grassmann) gradient, as supplied by the optimiser. + Reported by [`GradientGrassmann`](@ref). + - `bondresidual`: the change in the center bond tensor over a sweep. This is a fixed-point + residual: it says the sweeps have stopped moving, which is weaker than saying the state is + variationally stationary. Reported by [`IDMRG`](@ref) and [`IDMRG2`](@ref). + - `localchange`: the largest relative change of a local tensor over a sweep. Reported by + [`DMRG`](@ref) and [`DMRG2`](@ref) inside [`approximate`](@ref). + +[`convergence_measure`](@ref) returns whichever of these is present, for code that only wants +"the number that was compared against `tol`" without caring which one it is. + +### Truncation + +Both truncation entries are built from the same per-factorisation quantity, namely the 2-norm of +the singular values a single local factorisation discarded, but aggregate it differently, because +no single aggregation answers every question: + + - `max_truncation_error`: the largest of them. It is still a per-factorisation quantity rather + than a combination of them, so it does not grow with system size or iteration count, which is + what makes it comparable between runs. It is also the entry a `trunc` setting most directly + controls, though how directly depends on the strategy. + - `total_truncation_error`: all of them combined in quadrature, + ``\\sqrt{\\sum_k \\epsilon_k^2}``. This grows with system size and iteration count, so unlike + `max_truncation_error` it is not comparable between runs. + - `numtrunc`: how many of the recorded errors were non-zero, i.e. how many actually discarded + anything. + +The two error entries also read under the short aliases `ϵ_max` and `ϵ_total` (`info.ϵ_max`, +`info[:ϵ_total]`, `haskey(info, :ϵ_max)`). They are only ever stored under the descriptive names, +so `keys` and displaying/showing them returns one name per quantity. + +Which factorisations get recorded is not the same for every algorithm, and this is worth knowing +before comparing `numtrunc` (or `total_truncation_error`) between them: + + - [`IDMRG`](@ref), [`IDMRG2`](@ref), [`TDVP2`](@ref), [`BUG`](@ref) and [`Zipup`](@ref) record + every factorisation as it happens, so `numtrunc` is a count of factorisations. A sweep that + visits a bond twice contributes twice. + - [`DMRG`](@ref) and [`DMRG2`](@ref) instead keep one slot per update position, overwritten as + the sweep passes, and record those slots once at the end. `numtrunc` is therefore the number of + positions whose most recent cut discarded something . This is never more than `length(ψ)` for + [`DMRG`](@ref) or `length(ψ) - 1` for [`DMRG2`](@ref), however many sweeps ran and however many + SVDs each performed. + +The sweeping choice is deliberate: what the returned state still throws away at a bond is the last +cut made there, not the sum of every cut ever made there. It does mean `numtrunc` counts different +things in the two families, so read it as "how many recorded errors were non-zero" rather than as a +tally of SVD calls. See [Aggregating truncation errors](@ref) for how the two relate to a `trunc` setting, and for the per-strategy caveats. -These are `0` for an algorithm that never truncates, which does not mean the result is exact. Rather, -it means this particular error channel is absent. See the manual on [Errors and accuracy](@ref) -for what is *not* measured here. +An algorithm that truncates reports all three even on a run where it happened to discard nothing, +so `max_truncation_error == 0` means "truncated, but cut nothing away", whereas the entries being +absent altogether means the algorithm never truncates. Neither says the result is exact. See the +manual on [Errors and accuracy](@ref) for what is *not* measured here. """ -struct AlgorithmInfo{T <: Real} - converged::Union{Bool, Nothing} - normres::Union{T, Nothing} - ϵ_max::T - ϵ_total::T - numtrunc::Int - numiter::Int +struct AlgorithmInfo + data::Dict{Symbol, Any} end """ - AlgorithmInfo(; converged, normres, truncation, numiter) + AlgorithmInfo(; truncation = nothing, kwargs...) + +Build an [`AlgorithmInfo`](@ref) from the entries an algorithm actually produced. Every keyword +becomes an entry. -Keyword constructor, with every field defaulting to "not produced by this algorithm": no -convergence notion, and nothing truncated. `truncation` accepts a -[`TruncationAccumulator`](@ref), or is left out when the algorithm does not truncate. +A keyword whose value is `nothing` is omitted rather than stored. This is how an algorithm +reports a quantity it computes only on some branches: write `galerkin = measured ? g : nothing` +to leave the entry out where there is nothing to report, instead of assembling a different keyword +set per branch. A missing entry is an error to read, so pass `nothing` only where absence is +the meaning you intend. + +`truncation` is special: it accepts a [`TruncationAccumulator`](@ref) and expands into the +`max_truncation_error`/`total_truncation_error`/`numtrunc` entries. Leave it out for an algorithm +that does not truncate. +`numiter` defaults to `1` for the single-shot algorithms. """ -function AlgorithmInfo(; - converged = nothing, normres = nothing, truncation = nothing, numiter::Int = 1 - ) - T = _info_scalartype(normres, truncation) - acc = isnothing(truncation) ? TruncationAccumulator(T) : truncation - return AlgorithmInfo{T}( - converged, isnothing(normres) ? nothing : convert(T, normres), - convert(T, acc.ϵ_max), convert(T, sqrt(acc.ϵ_sq)), acc.numtrunc, numiter - ) +function AlgorithmInfo(; truncation = nothing, kwargs...) + data = Dict{Symbol, Any}() + for (key, value) in kwargs + isnothing(value) || (data[_canonical_key(key)] = value) + end + get!(data, :numiter, 1) + if !isnothing(truncation) + data[:max_truncation_error] = truncation.ϵ_max + data[:total_truncation_error] = sqrt(truncation.ϵ_sq) + data[:numtrunc] = truncation.numtrunc + end + return AlgorithmInfo(data) end -_info_scalartype(::Nothing, ::Nothing) = Float64 -_info_scalartype(normres, ::Nothing) = float(typeof(normres)) -_info_scalartype(::Nothing, acc) = _acc_type(acc) -_info_scalartype(normres, acc) = promote_type(float(typeof(normres)), _acc_type(acc)) +# the entries holding "the number that was compared against `tol`" +const convergence_keys = (:galerkin, :gradientnorm, :bondresidual, :localchange) + +""" + convergence_measure(info::AlgorithmInfo) + +The quantity that was compared against the algorithm's `tol`, whichever of +`$(join(convergence_keys, "`/`"))` the algorithm reported, or `nothing` for an algorithm that does +not iterate towards a fixed point and reports none of them. + +Use this when you only want the number, and read the specific entry when the kind of measure +matters, since these are not comparable with one another. +""" +function convergence_measure(info::AlgorithmInfo) + data = getfield(info, :data) + for key in convergence_keys + haskey(data, key) && return data[key] + end + return nothing +end + +# short aliases for the two truncation entries +# entries remain stored under the descriptive name +const _key_aliases = Dict{Symbol, Symbol}( + :ϵ_max => :max_truncation_error, + :ϵ_total => :total_truncation_error, +) +_canonical_key(key::Symbol) = get(_key_aliases, key, key) + +# dictionary interface +Base.getindex(info::AlgorithmInfo, key::Symbol) = getfield(info, :data)[_canonical_key(key)] +Base.haskey(info::AlgorithmInfo, key::Symbol) = haskey(getfield(info, :data), _canonical_key(key)) +function Base.get(info::AlgorithmInfo, key::Symbol, default) + return get(getfield(info, :data), _canonical_key(key), default) +end +Base.keys(info::AlgorithmInfo) = keys(getfield(info, :data)) +Base.values(info::AlgorithmInfo) = values(getfield(info, :data)) +Base.pairs(info::AlgorithmInfo) = pairs(getfield(info, :data)) +Base.length(info::AlgorithmInfo) = length(getfield(info, :data)) + +# property sugar: asking for an entry the algorithm never reported is an +# error naming what it did report, rather than a silent `nothing` +Base.propertynames(info::AlgorithmInfo) = Tuple(sort!(collect(keys(getfield(info, :data))))) +function Base.getproperty(info::AlgorithmInfo, key::Symbol) + key === :data && return getfield(info, :data) + data = getfield(info, :data) + canonical = _canonical_key(key) + haskey(data, canonical) && return data[canonical] + return _no_entry_error(info, canonical) +end + +@noinline function _no_entry_error(info::AlgorithmInfo, key::Symbol) + reported = join(propertynames(info), ", ") + msg = "this AlgorithmInfo has no entry `$key`; this algorithm reports $reported." + throw(ArgumentError(msg)) +end """ TruncationAccumulator{T} @@ -104,36 +213,69 @@ function push_error!(acc::TruncationAccumulator, ϵ) end # combining infos follows the same rules as combining the per-bond errors within one of them: -# worst case for `ϵ_max`, sum of squares for `ϵ_total`, and the later convergence verdict wins +# worst case for `ϵ_max`, sum of squares for `ϵ_total`, later convergence verdict wins, +# and counts are summed # assumes type stability of the scalars and an ordering of receiving this info +_later_wins(_, later) = later +const _combine_rules = Dict{Symbol, Any}( + :max_truncation_error => max, + :total_truncation_error => (a, b) -> sqrt(a^2 + b^2), + :numtrunc => +, + :numiter => +, +) + function _combine(a::AlgorithmInfo, b::AlgorithmInfo) - return AlgorithmInfo( - b.converged, b.normres, - max(a.ϵ_max, b.ϵ_max), sqrt(a.ϵ_total^2 + b.ϵ_total^2), - a.numtrunc + b.numtrunc, a.numiter + b.numiter - ) + dataₐ = copy(getfield(a, :data)) + for (key, value) in getfield(b, :data) + dataₐ[key] = haskey(dataₐ, key) ? + get(_combine_rules, key, _later_wins)(dataₐ[key], value) : value + end + return AlgorithmInfo(dataₐ) end +# custom show +# entries are displayed in a fixed order, with anything outside the known vocabulary listed last +const _show_order = (convergence_keys..., :max_truncation_error, :total_truncation_error) +const _show_handled = (_show_order..., :converged, :numiter, :numtrunc) + function Base.show(io::IO, ::MIME"text/plain", info::AlgorithmInfo) + data = getfield(info, :data) println(io, "AlgorithmInfo:") - if !isnothing(info.converged) - println(io, " converged = ", info.converged, " after ", info.numiter, " iterations") - else - println(io, " ", info.numiter, " iteration", info.numiter == 1 ? "" : "s") + numiter = get(data, :numiter, nothing) + tab_space = " " + + if haskey(data, :converged) + println( + io, tab_space, rpad("converged", 22), " = ", data[:converged], + " after ", numiter, " iterations" + ) + elseif !isnothing(numiter) + println(io, tab_space, numiter, " iteration", numiter == 1 ? "" : "s") end - isnothing(info.normres) || println(io, " normres = ", info.normres) - if info.numtrunc > 0 - println(io, " ϵ_max = ", info.ϵ_max, "\t(largest single factorization)") - println(io, " ϵ_total = ", info.ϵ_total, "\t(quadrature over ", info.numtrunc, ")") - else - println(io, " no truncation") + + for key in _show_order + haskey(data, key) || continue + suffix = if key === :max_truncation_error + "\t(largest single factorisation)" + elseif key === :total_truncation_error + "\t(quadrature over $(get(data, :numtrunc, 0)) truncations)" + else + "" + end + println(io, tab_space, rpad(string(key), 22), " = ", data[key], suffix) + end + haskey(data, :numtrunc) && data[:numtrunc] == 0 && println(io, " no truncation") + + for key in sort!(collect(keys(data))) + key in _show_handled && continue + println(io, tab_space, rpad(string(key), 22), " = ", data[key]) end return nothing end + function Base.show(io::IO, info::AlgorithmInfo) - return print( - io, "AlgorithmInfo(converged = ", info.converged, ", normres = ", info.normres, - ", ϵ_max = ", info.ϵ_max, ", ϵ_total = ", info.ϵ_total, - ", numtrunc = ", info.numtrunc, ", numiter = ", info.numiter, ")" - ) + data = getfield(info, :data) + print(io, "AlgorithmInfo(") + join(io, (string(key, " = ", data[key]) for key in sort!(collect(keys(data)))), ", ") + return print(io, ")") end diff --git a/test/algorithms/approximate.jl b/test/algorithms/approximate.jl index b9b58e25f..5685f79fe 100644 --- a/test/algorithms/approximate.jl +++ b/test/algorithms/approximate.jl @@ -121,7 +121,7 @@ end got, info = approximate((O, ψ), Zipup(; trunc = trunc′, left_to_right)) @test norm(ref - got) / norm(ref) < 1.0e-10 @test info.ϵ_max < 1.0e-10 - @test isnothing(info.converged) && isnothing(info.normres) + @test !haskey(info, :converged) && isnothing(convergence_measure(info)) end @test norm(ψ - ψ_copy) < 1.0e-12 diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl index 3269a981c..1ca222bd5 100644 --- a/test/algorithms/groundstate.jl +++ b/test/algorithms/groundstate.jl @@ -37,7 +37,7 @@ verbosity_conv = 1 v = variance(ψ, H) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 @test info.numtrunc == 0 @@ -65,7 +65,7 @@ verbosity_conv = 1 v = variance(ψ, H) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 @@ -92,7 +92,7 @@ verbosity_conv = 1 v = variance(ψ, H) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 # the bond should have grown to the truncation target @@ -120,7 +120,7 @@ verbosity_conv = 1 v = variance(ψ, H) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 # the bond should have grown to the truncation target @@ -146,7 +146,7 @@ verbosity_conv = 1 v = variance(ψ, H) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 # the bond should have grown to the truncation target @@ -194,7 +194,7 @@ verbosity_conv = 1 v = variance(ψ, H) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 @test v < v₀ && v < 1.0e-2 end end @@ -229,7 +229,7 @@ end v = variance(ψ′, H, envs) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 end @@ -247,8 +247,8 @@ end v = variance(ψ, H, envs) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 - @test info.numtrunc == 0 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 + @test !haskey(info, :numtrunc) # single-site IDMRG never truncates @test v < v₀ @test v < 1.0e-2 end @@ -270,7 +270,7 @@ end v = variance(ψ, H, envs) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 @@ -299,7 +299,7 @@ end v = variance(ψ′, H, envs) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 end @@ -315,7 +315,7 @@ end v = variance(ψ, H, envs) # test using low variance - @test info.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 @test v < v₀ @test v < 1.0e-2 end diff --git a/test/misc/multifusion.jl b/test/misc/multifusion.jl index b2df5e81a..1631cca3a 100644 --- a/test/misc/multifusion.jl +++ b/test/misc/multifusion.jl @@ -85,8 +85,8 @@ module TestMultifusion v2 = variance(ψ2, H) E2 = expectation_value(ψ2, H, envs2) - @test info.normres ≈ 0 atol = 1.0e-3 - @test info2.normres ≈ 0 atol = 1.0e-3 + @test convergence_measure(info) ≈ 0 atol = 1.0e-3 + @test convergence_measure(info2) ≈ 0 atol = 1.0e-3 @test v < v₀ && v2 < v₀ @test isapprox(E, E2; atol = 1.0e-6) @@ -123,7 +123,7 @@ module TestMultifusion @test isapprox(E, E2; atol = 1.0e-6) @test isapprox(E, E3; atol = 1.0e-6) - for delta in [info.normres, info2.normres, info3.normres] + for delta in [convergence_measure(info), convergence_measure(info2), convergence_measure(info3)] @test delta ≈ 0 atol = 1.0e-3 end for var in [v, v2, v3] From 056df3e25a353978f18ae9d5599e848ce02972ca Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 28 Aug 2026 14:44:56 +0200 Subject: [PATCH 28/29] cut text by relying on docstring, less theory more implementation, format consistency [skip ci] --- docs/src/assets/mpskit.bib | 27 ----- docs/src/man/algorithms.md | 236 ++++++++----------------------------- 2 files changed, 49 insertions(+), 214 deletions(-) diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index 28d9eb56e..fd1352050 100644 --- a/docs/src/assets/mpskit.bib +++ b/docs/src/assets/mpskit.bib @@ -1000,33 +1000,6 @@ @article{hubig2015 url = {https://link.aps.org/doi/10.1103/PhysRevB.91.155115} } -@article{hubig2018, - title = {Error estimates for extrapolations with matrix-product states}, - author = {Hubig, C. and Haegeman, J. and Schollw\"ock, U.}, - journal = {Phys. Rev. B}, - volume = {97}, - issue = {4}, - pages = {045125}, - year = {2018}, - month = {Jan}, - publisher = {American Physical Society}, - doi = {10.1103/PhysRevB.97.045125}, - url = {https://link.aps.org/doi/10.1103/PhysRevB.97.045125} -} - -@article{lubich2015, - title = {Time Integration of Tensor Trains}, - author = {Lubich, Christian and Oseledets, Ivan V. and Vandereycken, Bart}, - journal = {SIAM Journal on Numerical Analysis}, - volume = {53}, - number = {2}, - pages = {917--941}, - year = {2015}, - publisher = {Society for Industrial and Applied Mathematics}, - doi = {10.1137/140976546}, - url = {https://doi.org/10.1137/140976546} -} - @article{li2024, title = {Time-{{Dependent Variational Principle}} with {{Controlled Bond Expansion}} for {{Matrix Product States}}}, author = {Li, Jheng-Wei and Gleis, Andreas and {von Delft}, Jan}, diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index d4f05e7a3..2d6ecee0d 100644 --- a/docs/src/man/algorithms.md +++ b/docs/src/man/algorithms.md @@ -280,245 +280,107 @@ isapprox(Es[1] - E₀, 2(g - 1); rtol=1e-2) # infinite analytical result The algorithms that solve for a state, particularly [`find_groundstate`](@ref), [`leading_boundary`](@ref), [`approximate`](@ref), [`timestep`](@ref) and [`time_evolve`](@ref), return an [`AlgorithmInfo`](@ref) as their last value, describing how they arrived at their result. [`excitations`](@ref) and [`changebonds`](@ref) report nothing. What limits their accuracy is covered below all the same. -A single bare number could not do this job, because the quantities involved are genuinely different. -A convergence measure and a truncation error answer different questions, are not comparable, and not every algorithm produces both. -Conversely, two algorithms that both report "a convergence measure" do not necessarily report the same measure. -[`AlgorithmInfo`](@ref) therefore carries a set of named entries, and each algorithm fills in only the ones it actually computes. -Entries are keyed by `Symbol`, and read either as a property (`info.max_truncation_error`) or by indexing (`info[:max_truncation_error]`). -The usual dictionary interface is supported, in particular `haskey`, `keys`, `values` and `get`. -Nothing is promised that was never measured: asking for an entry an algorithm did not report is an error naming what it did report, rather than a value that was never computed, and the name of the entry says which quantity it is rather than leaving that to prose. -Display the returned object (or call `keys(info)`) to see what a given algorithm reported. -The docstrings of every algorithm report these as well. ```@docs; canonical=false AlgorithmInfo ``` -The rest of this section explains what those quantities are, and - equally important - what they do not measure. +The rest of this section explains what quantities can be reported by the algorithms, and - equally important - what they do not measure. ### The error convention -**In theory:** A truncation error measures how much a factorisation changed the tensor it acted on: replacing ``A`` by its rank-restricted approximation ``\tilde{A}`` costs +Every factorisation in MPSKit reports ``\epsilon = \lVert A - \tilde{A} \rVert``, which is the 2-norm of the discarded singular values ([Schollwöck](@cite schollwoeck2011)). +Interpreting this truncation error as a "discarded weight" is accurate when the factorised object is normalised. -```math -\epsilon = \lVert A - \tilde{A} \rVert . -``` - -Because a truncated SVD keeps the largest singular values, the discarded part is orthogonal to the kept part, and ``\epsilon`` is exactly the 2-norm of the discarded singular values ([Schollwöck](@cite schollwoeck2011)). - -**In practice:** That is precisely what MPSKit computes: every factorisation reports the 2-norm of what it discarded, and nothing more. In particular, the value is absolute, i.e. it is not normalised by ``\lVert A \rVert``. - -**Where the two are often conflated:** In DMRG the same quantity is usually called the "discarded weight" and is read as a probability: when ``A`` is the bond tensor of a normalised state, the ``\sigma_\alpha^2`` are eigenvalues of the reduced density matrix across the cut which sum to 1 for a normalised state, so ``\epsilon^2`` is the probability of the discarded subspace. -Even though viewing it as a probability is intuitively useful, it is not what makes the value well-defined, and it does not always apply here. -In particular, under non-renormalising algorithms the state norm drifts away from 1 precisely as truncation accumulates. -The probability intuition holds only when the factorised object is normalised, which is not always the case. -The definition ``\epsilon = \lVert A - \tilde{A} \rVert`` holds in every case, though. - -What differs between algorithms is how these per-factorisation values are summed up (*aggregated*) into the numbers they report, which is the subject of the next subsection. +What differs between algorithms is how these per-factorisation values are summed up (*aggregated*) into the numbers they report. !!! warning A convergence measure and a truncation error are unrelated quantities. - [`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms fill in a convergence entry (`galerkin`, `gradientnorm`, `bondresidual` or `localchange` depending on the algorithm), which has no truncation interpretation at all, while the truncating algorithms fill in `max_truncation_error`/`total_truncation_error`, which say nothing about convergence. An algorithm that does both fills both, and they should not be compared with each other. + Convergence measures are covered below per algorithm. #### Aggregating truncation errors -This applies to every truncating algorithm. -A sweep of [`DMRG2`](@ref), [`IDMRG2`](@ref), [`TDVP2`](@ref), [`BUG`](@ref) or [`Zipup`](@ref) performs many local factorisations, each with its own ``\epsilon_k``, and there is no single number that answers every question one might ask of them. -Rather than pick one and hope the caller wants that one, [`AlgorithmInfo`](@ref) carries both aggregations under names that say what they are. +The per-factorisation errors are aggregated two ways, as a worst case (`max_truncation_error`) and in quadrature (`total_truncation_error`). +See the [`AlgorithmInfo`](@ref) docstring for what each is. -`max_truncation_error`, the largest single ``\epsilon_k``, is a worst case. -The point of it is that it is still a per-factorisation quantity: it is one of the ``\epsilon_k``, not a combination of them, so it does not grow just because the chain is longer or more sweeps or steps were taken. -That is what makes it the one to watch over the course of a run, or to compare between runs at different sizes. -An increase in `max_truncation_error` means individual bonds are being cut harder, whereas an increase in `total_truncation_error` may only mean there were more bonds to cut. -This is the same reason the ground state algorithms report a maximum over local gradient norms rather than a total. -It is also why `max_truncation_error` is the one [`time_evolve`](@ref) logs; the ground state algorithms log their convergence measure instead, since for them convergence rather than truncation is what the sweep is driving. +`max_truncation_error` is the entry a `trunc` setting most directly controls, though how directly depends on the strategy: -It is also the field that a `trunc` setting most directly controls, though how directly depends on the strategy: +- [`truncerror`](@extref MatrixAlgebraKit.truncerror) bounds the discarded weight of each factorisation, which is exactly ``\epsilon_k``, so `max_truncation_error` should come out at or below the tolerance you set. +- [`trunctol`](@extref MatrixAlgebraKit.trunctol) bounds each individual singular value instead. Discarding ``k`` of them leaves ``\epsilon_k \le \sqrt{k}\,\texttt{atol}``, so `max_truncation_error` lands near the tolerance but is not bounded by it. +- [`truncrank`](@extref MatrixAlgebraKit.truncrank) fixes the rank and says nothing about magnitudes at all. Here, `max_truncation_error` is not something you set but something you read off. It is thus the consequence of that choice of bond dimension. -* [`truncerror`](@extref MatrixAlgebraKit.truncerror) bounds the discarded weight of each factorisation, which is exactly ``\epsilon_k``, so `max_truncation_error` should come out at or below the tolerance you set. -* [`trunctol`](@extref MatrixAlgebraKit.trunctol) bounds each individual singular value instead. Discarding ``k`` of them leaves ``\epsilon_k \le \sqrt{k}\,\texttt{atol}``, so `max_truncation_error` lands near the tolerance but is not bounded by it. -* [`truncrank`](@extref MatrixAlgebraKit.truncrank) fixes the rank and says nothing about magnitudes at all. Here, `max_truncation_error` is not something you set but something you read off. It is thus the consequence of that choice of bond dimension. - -`total_truncation_error`, all of them combined in quadrature, +`total_truncation_error` sums the squares, ```math \epsilon_{\text{total}} = \sqrt{\textstyle\sum_k \epsilon_k^2} , ``` -is the one that adds up to something. -Note that it is the *squares* that are summed, and the root taken at the end, because ``\epsilon^2`` is the quantity that accumulates exactly: each local truncation is an orthogonal projection, so it removes exactly ``\epsilon_k^2`` from the squared norm of the tensor it acts on. -Summing the squares therefore tracks a running "cost" rather than merely a worst case. -The price is that it is extensive: it grows like ``\sqrt{N}`` in the number of truncations, so it is not comparable between different system sizes or sweep counts. - -Whether that running cost is also the error of the *final state* depends on what the algorithm does between truncations, so it is not a property of the aggregation itself. -It does hold for real-time evolution, which is worked out in [The norm as a record of truncation](@ref). -A variational sweep, by contrast, renormalises as it goes, so there `total_truncation_error` is a diagnostic of how hard the truncation is working rather than a norm deficit. - -Neither field is the distance to the untruncated solution. -Bounding that gives the linear sum ``\lVert \psi_{\text{untruncated}} - \psi \rVert \le \sum_k \epsilon_k``, which is a third quantity again, and always the largest of the three. - -##### Which factorisations are recorded - -Both aggregations and the `numtrunc` count are built from the set of per-factorisation errors the algorithm records. -This set is not the same for every algorithm. - -[`IDMRG`](@ref), [`IDMRG2`](@ref), [`TDVP2`](@ref), [`BUG`](@ref) and [`Zipup`](@ref) record every factorisation as it happens. -Thus, `numtrunc` is a count of factorisations, and a sweep that visits a bond twice contributes two entries. +which tracks a running cost rather than a worst case. +Whether that cost is also the error of the *final state* depends on what the algorithm does between truncations: in real time evolution, where truncations do not normalise by default, it is exactly the norm deficit of the state. -[`DMRG`](@ref) and [`DMRG2`](@ref) keep one slot per update position instead, overwritten as the sweep passes over it, and record the slots once at the end. -In these cases, `numtrunc` is the number of positions whose most recent cut discarded something, and it can never exceed `length(psi)` (single-site) or `length(psi) - 1` (two-site), no matter how many sweeps ran. -A reported `numtrunc` less than the number of bonds means some bonds last discarded a non-zero weight. -These are typically the bonds nearest the two ends. - -The difference here is made deliberately. -What a returned variational state still throws away at a given bond is the last cut made there, not the sum of every cut ever made there. -This means that recording the latest value per position is the more meaningful quantity for a sweeping ground-state algorithm. -Thus, `numtrunc` counts different things in the two families, and `total_truncation_error` correspondingly sums a different number of terms, so neither is directly comparable between, say, [`DMRG2`](@ref) and [`IDMRG2`](@ref). -Read `numtrunc` more as "how many recorded errors were non-zero", and less as a tally of SVD calls. - -The two error entries also read under the short aliases `ϵ_max` and `ϵ_total`. -They are only ever stored under the descriptive names, so `keys` and displaying/showing them returns one name per quantity. +Which factorisations an algorithm records into these differs per family, which is why `numtrunc` and `total_truncation_error` are not comparable across algorithms. ### Ground state accuracy [`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms report the quantity their `tol` is compared against, together with a `converged` flag. Because these are not the same quantity from one algorithm to the next, each is stored under a key that names it (`galerkin`, `gradientnorm`, `bondresidual` or `localchange`). [`convergence_measure`](@ref) returns whichever of them is present, for code that only wants the number. -**In theory:** Convergence is measured by the (norm of the) variational gradient: the component of ``H \lvert \psi \rangle`` that points away from the current state but still lies in the tangent space of the variational manifold. -It vanishes exactly at a variational fixed point, and its norm is what most of these algorithms report. -For the sweeping algorithms that norm is known as the Galerkin error. - -**In practice:** The sweeping algorithms ([`DMRG`](@ref), [`DMRG2`](@ref), [`VUMPS`](@ref), [`VOMPS`](@ref)) report the Galerkin error, computed per site as the norm of the local update projected onto the orthogonal complement of the current tensor, and then aggregated as the maximum over sites. -[`GradientGrassmann`](@ref) instead reports the gradient norm supplied by its optimiser, taken over the whole state at once in the preconditioned Grassmann metric. - -[`IDMRG`](@ref) and [`IDMRG2`](@ref) report neither, and this is easy to miss because their `tol` sits alongside the others. -Their `bondresidual` is a *fixed-point residual*: the change in the center bond tensor from one sweep to the next, ``\lVert C - C_{\text{prev}} \rVert``, with both tensors projected onto their common space when the bond space changed. -It measures how much a sweep still moves the state, which is weaker than measuring how far the state is from stationarity, since an algorithm crawling through a slow region reports a small value for the same reason a converged one does. -For [`IDMRG2`](@ref) the projection onto the common space also means that a change in bond dimension between sweeps is projected out of the measure rather than counted in it. - -**Why they differ:** The Galerkin and Grassmann measures are the same underlying gradient, not different physical quantities. -They differ in the metric it is measured in (the Grassmann gradient is preconditioned) and in how the per-site contributions are combined (a maximum versus a single global norm). -The maximum is a deliberate practical choice, as it keeps the reported number independent of system size. -The IDMRG residual is not that gradient at all: it certifies that the sweeps have stopped moving, not that there is nowhere left to move to. -The consequence of these mismatches is that a `tol` tuned for one algorithm is not a `tol` tuned for another, and this is worth keeping in mind when swapping algorithms at a fixed `tol`. - -In other words, convergence is only defined relative to the manifold you are optimising over. -A single-site algorithm at a fixed bond dimension can drive its convergence measure to machine precision and still be far from the true ground state, because the error that remains is the bond dimension itself, which no amount of further sweeping can address. -A small convergence measure certifies a fixed point, not an accurate state. +Importantly, they represent different things, and a `tol` tuned for one algorithm is not a `tol` tuned for another. + +A single-site algorithm at a fixed bond dimension can drive its convergence measure to machine precision and still be far from the true ground state. Growing the bond dimension is the job of the two-site algorithms ([`DMRG2`](@ref), [`IDMRG2`](@ref)) or of a bond expansion ([`DMRG`](@ref) with an `alg_expand`, or an expanding `alg_gauge` such as [`DMRG3S`](@ref)); see also [`changebonds`](@ref). Once an algorithm does truncate, the two error notions interact. -The Galerkin error cannot fall below the level set by the weight being discarded each sweep, so a truncating scheme converges once `galerkin` reaches the truncation error rather than the (unreachable) bare `tol`. -[`DMRG`](@ref)/[`DMRG2`](@ref) account for this: their stopping test is `ϵ ≤ max(tol, maximum(ϵ_trunc))`, which reduces to the plain `ϵ ≤ tol` when nothing is truncated. +In the case of the Galerkin error, it cannot fall below the level set by the weight being discarded each sweep, so a truncating scheme converges once `galerkin` reaches the truncation error rather than the (unreachable) bare `tol`. Neither measure is an error bar on an observable, and no cheap substitute for one exists. The energy variance ``\langle H^2 \rangle - \langle H \rangle^2`` is an independent and more demanding measure. Note what it actually quantifies, namely how far the state is from being an *exact eigenstate*, which is not the same thing as the error on some other observable. -The usual practical route is to compute a series of states at increasing bond dimension and extrapolate observables against the variance towards zero ([Hubig et al.](@cite hubig2018)). -This is an empirical extrapolation rather than a bound. ### Time evolution accuracy Unlike a ground state search, a time evolution has no convergence criterion to run to. There is no fixed point, and the error is made at every step. -Three sources behave differently and only one of them is reported. - -* **Truncation error.** - Whenever a bond is cut back down, the discarded singular values are lost from the state. - [`timestep`](@ref) and [`time_evolve`](@ref) report this through the `max_truncation_error`/`total_truncation_error` entries of the [`AlgorithmInfo`](@ref) they return. - It is the error you control through the algorithm's `trunc`, and the only one that is free to compute, since the truncating SVD produces it anyway. - It is non-zero for [`TDVP2`](@ref), for [`BUG`](@ref) with a `trunc`, and for [`TDVP`](@ref) with a bond expansion. - Plain single-site [`TDVP`](@ref) runs at fixed bond dimension and reports exactly `0`. - -* **Projection error.** - Single-site [`TDVP`](@ref) confines the evolution to the tangent space of a fixed-bond-dimension manifold, ``\lVert (1 - P_{T_\psi}) H \psi \rVert``. - The component of the exact evolution pointing off that manifold is simply dropped, and this happens even with no truncation and exact local solves. - It is not reported, since measuring it costs an extra effective-Hamiltonian application per site. - This is what a bond expansion (CBE) exists to reduce ([Li et al.](@cite li2024)). - -* **Splitting error.** - The evolution is not applied in one piece. - Rather, it is split into local terms that are integrated in sequence, which do not commute. - This is a Trotter-type error. - Note, however, that what is split differs per method: [`TDVP`](@ref) splits the tangent-space projector into site and bond terms, while the MPO methods ([`WI`](@ref), [`WII`](@ref), [`TaylorCluster`](@ref)) split the Hamiltonian in the more familiar sense. - For TDVP's symmetric back-and-forth sweep the result is globally ``O(dt^2)`` ([Lubich et al.](@cite lubich2015), [Paeckel et al.](@cite paeckel2019)), so it is controlled by `dt` alone. - This can only be estimated by comparing one step of `dt` against two of `dt / 2`. - -A trustworthy run needs all three under control, not just a small truncation error. -In practice: pick `dt` from a convergence check, pick `trunc` from the reported truncation error, and use a bond-adaptive scheme ([`TDVP2`](@ref), [`BUG`](@ref), or [`TDVP`](@ref) with `alg_expand`) whenever entanglement grows during the evolution, since a fixed bond dimension silently converts entanglement growth into projection error. - -Note that these do not all shrink together, so there is a sweet spot in `dt` rather than a "smaller is always better". -Decreasing `dt` reduces the splitting error, but it also means more steps to reach the same final time, and every step truncates again. -The accumulated truncation error grows with the number of steps, as does the compute time. -The projection error does not improve at all, since it represents a rate at which the exact solution leaves the manifold. -Shrinking `dt` in this case only samples this rate more finely. -The practical consequence is that below some `dt` the total error stops improving and eventually gets worse, and the remedy at that point is a larger bond dimension rather than a smaller step. - -#### The norm as a record of truncation - -How the per-factorisation errors are aggregated into `max_truncation_error`/`total_truncation_error` is described under [Aggregating truncation errors](@ref); what follows is specific to time evolution. - -By default none of the time evolution algorithms renormalize (`normalize = false`), which is deliberate. -In real time the local exponentials are unitary, so truncation is the only thing that changes the norm and it becomes a running record of what truncation has cost, - -```math -\lVert \psi \rVert^2 = \lVert \psi_0 \rVert^2 - \epsilon_{\text{total}}^2 . -``` +Time evolution has three distinct error sources, namely the truncation error, the projection error, and the splitting error. +Only the truncation error is reported in [`timestep`](@ref) and [`time_evolve`](@ref)'s [`AlgorithmInfo`](@ref). +This is non-zero for [`TDVP2`](@ref), for [`BUG`](@ref) with a `trunc`, and for [`TDVP`](@ref) with a bond expansion. -The reported `total_truncation_error` is the norm deficit, and this composes across steps. -This follows from the following two facts put together, one per half of a local update. -1) An SVD truncation is an orthogonal projection onto the kept Schmidt vectors and is 2-norm optimal at that rank ([Schollwöck](@cite schollwoeck2011)), so the kept and discarded parts are orthogonal. -By Pythagoras the squared norm drops by exactly the discarded weight, the usual way of quantifying truncation during a time evolution ([Paeckel et al.](@cite paeckel2019)). -2) The local exponentials of the projector-splitting sweep are unitary, so TDVP conserves the norm and the energy exactly when the local equations are solved exactly ([Paeckel et al.](@cite paeckel2019)), contributing nothing to the norm change. +The projection error is not reported, since measuring it costs an extra effective-Hamiltonian application per site. +This is what a bond expansion (CBE) exists to reduce ([Li et al.](@cite li2024)). +The splitting error is a Trotter-type error, and can only be estimated by comparing one step of `dt` against two of `dt / 2`. -These two hold for the local updates of a step, so composing over all steps gives the identity. +All three need to be under control, not just the reported one. +In practice: pick `dt` from a convergence check, pick `trunc` from the reported truncation error, and use a bond-adaptive scheme ([`TDVP2`](@ref), [`BUG`](@ref), or [`TDVP`](@ref) with `alg_expand`) whenever entanglement grows, since a fixed bond dimension silently converts entanglement growth into projection error. -!!! note - "Exactly" in the second fact is up to the tolerance of the local exponentials, and is thus in practice only approximate due to integrator tolerance. - -It is also specific to real time with `normalize = false`: - -* **Imaginary time** evolves with the non-unitary ``\exp(-H dt)``, which rescales the state on its own. - The norm then moves for two independent reasons, namely the physical decay of the weight and the truncation loss. - One cannot separate them from each other. - The reported errors still count only the truncation. -* **`normalize = true`** renormalizes at every local update, destroying the identity by construction. - This is usually what you want for imaginary-time evolution used as a ground state or thermal-state search. - The reported errors are unaffected. -* **No truncation at all** (plain single-site [`TDVP`](@ref), or [`BUG`](@ref) with a QR gauge) gives ``\epsilon = 0``, and in real time the norm is then conserved exactly. -* An **`InfiniteMPS`** is regauged to norm 1 per site structurally, so its norm carries no such information and `normalize` has no effect. +They do not shrink together, so there is a sweet spot in `dt` rather than "smaller is better". +This is because a smaller `dt` lowers the splitting error but takes more steps to reach the same time, and every step truncates again. ### Excitation accuracy [`excitations`](@ref) returns only `(energies, states)`: there is no error term, and none of the sources below is reported back to you. They are worth knowing about, because the dominant one is usually not the one the algorithm is working on. -* **Inherited ground state error.** - Every method here builds on a ground state you supply and treats it as exact. - Its error propagates straight into the gap, and since a gap is a difference of two large energies, it is typically the limiting factor. - A well-converged ground state (in the sense of its convergence measure and bond dimension) is necessary for a meaningful excitation calculation. - -* **Ansatz limitation.** - [`QuasiparticleAnsatz`](@ref) varies over the single-quasiparticle tangent space on top of a fixed ground state. - It is variational within that space and well suited to isolated quasiparticle branches, but multi-particle continua are not representable in it, so results there are not to be trusted. - How well it does on an isolated branch is controlled by the gaps around the targeted level. - The ansatz approximates an eigenvalue that is separated from the rest of the spectrum in its momentum sector with an error bounded exponentially in the size of the local operator's support, at a rate set by the gap below *and* above that eigenvalue ([Haegeman et al.](@cite haegeman2013)). - A branch that is nearly degenerate with the ground state, or that sits just below a continuum, therefore converges much more slowly than an isolated one. - -* **Eigensolver convergence.** - The eigenvalue problem is solved with KrylovKit, and a run that fails to converge `num` states emits a warning carrying the residual when the verbosity is set high enough. - That residual is neither returned nor thrown, so it is worth not running with warnings suppressed. - Its convergence is governed by the same spectral structure as the ansatz above: levels that are well separated converge quickly, while nearly degenerate ones converge slowly and are the ones most likely to come back unconverged. - -* **Penalty-based orthogonality** ([`FiniteExcited`](@ref)). - Higher states are found by minimising ``H + \lambda \sum_i |\psi_i\rangle\langle\psi_i|`` against the previously converged states, with ``\lambda`` the `weight` field. - A finite `weight` enforces orthogonality only approximately, so a residual overlap with a lower state biases the reported energy downwards. - Since the reported value is the expectation value of the bare `H`, this bias is invisible in the output. - Raising `weight` suppresses it at the cost of stretching the spectrum and slowing down the eigensolver's per-gap eigensolves. - -* **Truncation** ([`ChepigaAnsatz2`](@ref)). - The two-site excited state is split back to single-site tensors with a truncated SVD governed by `trunc`, and the resulting discarded weight is not reported. +- **Inherited ground state error.** + Every method builds on the ground state you supply and treats it as exact. + Since a gap is a difference of two large energies, that error propagates straight into it and is typically the limiting factor. + +- **Ansatz limitation.** + [`QuasiparticleAnsatz`](@ref) varies over the single-quasiparticle tangent space on top of a fixed ground state, so it is variational within that space and suited to isolated quasiparticle branches. + Its error is bounded exponentially in the support of the local operator, at a rate set by the gaps below *and* above the targeted eigenvalue ([Haegeman et al.](@cite haegeman2013)). + +- **Eigensolver convergence.** + The eigenvalue problem is solved with KrylovKit, and a run that fails to converge `num` states emits a warning carrying the residual when the verbosity is high enough. + That residual is neither returned nor thrown, so it is worth not suppressing warnings. + Nearly degenerate levels are the ones most likely to come back unconverged. + +- **Penalty-based orthogonality.** + [`FiniteExcited`](@ref) minimises ``H + \lambda \sum_i |\psi_i\rangle\langle\psi_i|`` against the previously converged states, with ``\lambda`` the `weight` field. + A finite `weight` enforces orthogonality only approximately, so a residual overlap with a lower state biases the energy downwards — invisibly, since the reported value is the expectation value of the bare `H`. + Raising `weight` suppresses the bias at the cost of stretching the spectrum and slowing the eigensolver. + +- **Truncation** ([`ChepigaAnsatz2`](@ref)). + The two-site excited state is split back to single-site tensors with a truncated SVD governed by `trunc`, and the discarded weight is not reported. ## `changebonds` From e5c8010170076c264889d9b25b51b7915a31c007 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 18 Sep 2026 18:37:10 +0200 Subject: [PATCH 29/29] one more if-else for spacelists --- test/groundstate/approximate.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/groundstate/approximate.jl b/test/groundstate/approximate.jl index 3162b553a..11c5590af 100644 --- a/test/groundstate/approximate.jl +++ b/test/groundstate/approximate.jl @@ -14,11 +14,11 @@ using Random verbosity_conv = 1 # fixtures for the `Zipup` testsets -zipup_spacelist = [ - (ℙ^4, ℙ^3, 4), - (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1), 8), -] -fast_tests && (zipup_spacelist = zipup_spacelist[1:1]) +zipup_spacelist = if fast_tests + [(ℙ^4, ℙ^3, 4)] +else + [(ℙ^4, ℙ^3, 4), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1), 8)] +end function _random_mpo_mps(pspace, Dspace, L; elt = ComplexF64) Random.seed!(1357)