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/assets/mpskit.bib b/docs/src/assets/mpskit.bib
index 44e762933..fd1352050 100644
--- a/docs/src/assets/mpskit.bib
+++ b/docs/src/assets/mpskit.bib
@@ -999,3 +999,31 @@ @article{hubig2015
   doi = {10.1103/PhysRevB.91.155115},
   url = {https://link.aps.org/doi/10.1103/PhysRevB.91.155115}
 }
+
+@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/changelog.md b/docs/src/changelog.md
index 893a670e7..ac050c7f7 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
@@ -45,6 +45,20 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   virtual channels between terms that start out with the same operators, up to a scalar factor,
   and add up terms that are linearly dependent. The resulting Hamiltonian is unchanged, but its
   bond dimension is generally smaller ([#518](https://github.com/QuantumKitHub/MPSKit.jl/pull/518))
+- 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
+  `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,
+  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,
@@ -120,6 +134,17 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - Reorganised the test suite to reduce CI wall time, as well as added the `--fast` test flag
   to test fewer sector and scalar types. ([#517](https://github.com/QuantumKitHub/MPSKit.jl/pull/517))
 
+- `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. ([#512](https://github.com/QuantumKitHub/MPSKit.jl/pull/512))
+  - 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.
+    ([#512](https://github.com/QuantumKitHub/MPSKit.jl/pull/512))
+
 ## [0.13.11](https://github.com/QuantumKitHub/MPSKit.jl/compare/v0.13.10...v0.13.11) - 2026-05-04
 
 ### Added
diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md
index 5ff913006..d68509466 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,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.
@@ -32,7 +34,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.
@@ -48,7 +50,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
@@ -100,6 +102,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
@@ -118,10 +125,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
@@ -262,6 +274,113 @@ Es, ϕs = excitations(H, ChepigaAnsatz2(), ψ, envs; num=1)
 isapprox(Es[1] - E₀, 2(g - 1); rtol=1e-2) # infinite analytical result
 ```
 
+## Errors and accuracy
+
+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.
+
+```@docs; canonical=false
+AlgorithmInfo
+```
+
+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
+
+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.
+
+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.
+    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
+
+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` is the entry 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.
+
+`total_truncation_error` sums the squares,
+
+```math
+\epsilon_{\text{total}} = \sqrt{\textstyle\sum_k \epsilon_k^2} ,
+```
+
+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.
+
+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.
+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.
+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.
+
+### 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.
+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 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`.
+
+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.
+
+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 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`
 
 Many of the previously mentioned algorithms do not possess a way to dynamically change to
@@ -274,15 +393,18 @@ state.
 changebonds
 ```
 
+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
 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/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)
 
diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index f900884ae..c10608acd 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, convergence_measure
 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/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl
index 7cfb5dc12..d788e6dff 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,
@@ -29,19 +29,47 @@ 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 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
   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)
+- `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` 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:
+    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.
+
 # Algorithms
 
 Each algorithm below only supports a subset of the general interface. Check this table before
 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/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl
index 47dcfdfef..5ca2a7045 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, localchange = ϵ, 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, localchange = ϵ, numiter = iter)
 end
diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl
index 7503ac76d..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, ϵ
+    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, numiter = iter)
 end
 
 function approximate!(
@@ -72,11 +72,12 @@ 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, ϵ)
         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, bondresidual = ϵ, truncation = acc, numiter = iter)
+    return ψ, envs, info
 end
diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl
index 77e6c8cde..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, ϵ
+                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, ϵ
+                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 06b6ec401..6e5f4bc57 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.
+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!
@@ -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
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index b6aa8ef70..2ca815147 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -60,7 +60,9 @@ $(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 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"
@@ -171,7 +173,9 @@ $(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 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"
@@ -232,7 +236,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)
@@ -256,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.
@@ -273,7 +277,14 @@ 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
+- `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, ψ)
@@ -300,10 +311,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
@@ -352,7 +364,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)), galerkin = ϵ_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/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl
index 29d37271f..67dc3e8ab 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,7 +38,14 @@ 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
+- `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
 
@@ -57,7 +64,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/groundstate/gradient_grassmann.jl b/src/algorithms/groundstate/gradient_grassmann.jl
index 550edaac4..2a889b294 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, 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
@@ -79,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} 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!(ψ)
@@ -125,5 +127,10 @@ function _find_groundstate(ψ, H, alg::GradientGrassmann, envs, scheduler, timer
         @infov 4 TimerReport(timeroutput)
     end
 
-    return x, envs, normgradhistory[end]
+    normres = normgradhistory[end, 2] # full history returned as [fhistory normgradhistory]
+    info = AlgorithmInfo(;
+        converged = normres <= alg.method.gradtol, gradientnorm = normres,
+        numiter = size(normgradhistory, 1) - 1 # history starts with initial point before first iteration
+    )
+    return x, envs, info
 end
diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl
index 1b9c3c494..31a6aa2d3 100644
--- a/src/algorithms/groundstate/idmrg.jl
+++ b/src/algorithms/groundstate/idmrg.jl
@@ -12,7 +12,10 @@ $(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 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"
@@ -45,7 +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
-    "tolerance for convergence criterium"
+    "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"
@@ -77,16 +83,18 @@ struct IDMRGState{S, O, E, T, TO, A}
     envs::E
     iter::Int
     ϵ::Float64 # TODO: Could be any <:Real
+    truncation::TruncationAccumulator{Float64} # of the most recent sweep only
     energy::T
     timeroutput::TO
     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::TO, allocator::A,
     ) where {S, O, E, T, TO, A}
     return IDMRGState{S, O, E, T, TO, A}(
-        mps, operator, envs, iter, ϵ, T(energy), timeroutput, allocator
+        mps, operator, envs, iter, ϵ, truncation, T(energy), timeroutput, allocator
     )
 end
 
@@ -113,7 +121,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
@@ -134,7 +143,12 @@ 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, bondresidual = it.state.ϵ,
+            truncation = alg isa IDMRG2 ? it.state.truncation : nothing,
+            numiter = it.state.iter
+        )
+        return ψ′, envs, info
     end
 end
 
@@ -142,7 +156,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]
@@ -163,14 +178,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!(
@@ -180,12 +195,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
@@ -229,7 +244,7 @@ function _localupdate_sweep_idmrg!(
 end
 
 function _localupdate_sweep_idmrg2!(
-        ψ, H, envs, alg_eigsolve, alg_trunc, alg_svd, timeroutput;
+        ψ, H, envs, alg_eigsolve, alg_trunc, alg_svd, timeroutput, acc::TruncationAccumulator;
         backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
     )
     # @timeit wraps its body in try-finally, which is a new lexical scope: declare locals
@@ -243,7 +258,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
@@ -266,7 +282,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
@@ -294,7 +311,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
@@ -318,7 +336,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/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl
index 7be8c096a..b4ea716eb 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), reported as the `galerkin` entry of the returned [`AlgorithmInfo`](@ref)"
     tol::Float64 = Defaults.tol
 
     "maximal amount of iterations"
@@ -82,12 +83,12 @@ function dominant_eigsolve(
             if ϵ ≤ alg.tol
                 @infov 4 TimerReport(timeroutput)
                 @infov 2 logfinish!(log, it.iter, ϵ, expectation_value(mps, operator, envs))
-                return mps, envs, ϵ
+                return mps, envs, AlgorithmInfo(; converged = true, galerkin = ϵ, numiter = it.iter)
             end
             if it.iter ≥ alg.maxiter
                 @infov 4 TimerReport(timeroutput)
                 @warnv 1 logcancel!(log, it.iter, ϵ, expectation_value(mps, operator, envs))
-                return mps, envs, ϵ
+                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/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/gradient_grassmann.jl b/src/algorithms/statmech/gradient_grassmann.jl
index 8084cbda2..67e1422d7 100644
--- a/src/algorithms/statmech/gradient_grassmann.jl
+++ b/src/algorithms/statmech/gradient_grassmann.jl
@@ -19,5 +19,11 @@ function leading_boundary(
         alg.finalize!,
         isometrictransport = true
     )
-    return x, envs, normgradhistory[end]
+
+    normres = normgradhistory[end, 2] # full history returned as [fhistory normgradhistory]
+    info = AlgorithmInfo(;
+        converged = normres <= alg.method.gradtol, gradientnorm = normres,
+        numiter = size(normgradhistory, 1) - 1 # history starts with initial point before first iteration
+    )
+    return x, envs, info
 end
diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl
index 84edab763..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, ϵ
+    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, numiter = iter)
 end
 
 function leading_boundary(
@@ -69,11 +69,12 @@ 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, ϵ)
         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
@@ -156,7 +160,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
@@ -196,5 +201,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, bondresidual = ϵ, truncation = acc, numiter = iter)
 end
diff --git a/src/algorithms/statmech/leading_boundary.jl b/src/algorithms/statmech/leading_boundary.jl
index 218a3f8f7..6960bfa3f 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.
@@ -14,7 +14,8 @@ 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 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
 
@@ -22,7 +23,20 @@ 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
+- `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 `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).
 """ leading_boundary
 
 # TODO: alg selector
diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl
index 9427f5488..626d17b05 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), reported as the `galerkin` entry of the returned [`AlgorithmInfo`](@ref)"
     tol::Float64 = Defaults.tol
 
     "maximal amount of iterations"
@@ -73,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, galerkin = ϵ, 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, 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 00916ab90..990d85edf 100644
--- a/src/algorithms/timestep/bug.jl
+++ b/src/algorithms/timestep/bug.jl
@@ -23,10 +23,19 @@ 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_{\\text{total}}^2``, 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 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
+    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.
+
 !!! 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 +148,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 +163,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 +181,27 @@ function timestep!(
     # the sweep is serial, so a single allocator serves all local updates
     allocator = default_allocator(ψ, SerialScheduler())
 
+    acc = TruncationAccumulator(ψ)
+
     # 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
         )
+        push_error!(acc, ϵ)
     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
         )
+        push_error!(acc, ϵ)
     end
 
-    return ψ, envs
+    return ψ, envs, AlgorithmInfo(; truncation = acc)
 end
 
 # copying version
diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl
index a3da18c52..83ccbdb1f 100644
--- a/src/algorithms/timestep/tdvp.jl
+++ b/src/algorithms/timestep/tdvp.jl
@@ -11,11 +11,17 @@ 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_{\\text{total}}^2``,
+    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.
+
+    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
 
@@ -129,7 +135,10 @@ function _timestep_infinite(
     end
 
     recalculate!(envs, ψ′, H)
-    return ψ′, envs
+    # 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
 
 function timestep!(
@@ -148,6 +157,8 @@ function _timestep_finite!(
         ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs, allocator;
         imaginary_evolution::Bool, normalize::Bool
     )
+    acc = TruncationAccumulator(ψ)
+
     # sweep left to right
     for i in 1:(length(ψ) - 1)
         # 1. optionally expand the bond ahead of the local update (CBE)
@@ -161,7 +172,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)
+        push_error!(acc, ϵ)
 
         # 4. evolve the bond tensor backward
         Hc = C_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator)
@@ -190,7 +202,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)
+        push_error!(acc, ϵ)
 
         # 4. evolve the bond tensor backward
         Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator)
@@ -207,13 +220,14 @@ function _timestep_finite!(
         imaginary_evolution
     )
 
-    return ψ, envs
+    return ψ, envs, AlgorithmInfo(; truncation = acc)
 end
 
 """
 $(TYPEDEF)
 
 Two-site MPS time-evolution algorithm based on the Time-Dependent Variational Principle.
+See [`TDVP`](@ref) for more information.
 
 # Fields
 
@@ -266,16 +280,20 @@ 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)
+
+    acc = TruncationAccumulator(ψ)
+
     # 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 norm of the discarded singular values is the truncation error
+        _, ϵ = gauge2!(ψ, i, Val(:right), ac2′, alg_gauge; normalize)
+        push_error!(acc, ϵ)
 
         if i != (length(ψ) - 1)
             Hac = AC_hamiltonian(i + 1, ψ, H, ψ, envs; alg.backend, allocator)
@@ -292,10 +310,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)
+        push_error!(acc, ϵ)
 
         if i != 2
             Hac = AC_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator)
@@ -306,7 +322,7 @@ function _timestep2_finite!(
         end
     end
 
-    return ψ, envs
+    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 68014ec7f..6483952dd 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,6 +26,15 @@ 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, 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_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
 
@@ -35,27 +44,33 @@ 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))))
+        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)
                 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)}
+                info = _combine(info, info_step)
 
-                @infov 3 logiter!(log, iter, 0.0, t)
+                # log the size-independent error measure
+                # 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), 0.0, 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
+        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 ψ``.
@@ -81,6 +96,23 @@ solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``.
 
 - `ψ`: the time-stepped state
 - `envs`: the updated environment manager
+- `info::AlgorithmInfo`: what the step truncated (see below)
+
+# Truncation error
+
+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_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). 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 entries, and [Time evolution accuracy](@ref) in the manual
+for the other error sources.
 
 # Examples
 
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 68b062d35..b69665a27 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
diff --git a/src/utility/algorithminfo.jl b/src/utility/algorithminfo.jl
new file mode 100644
index 000000000..7021768cf
--- /dev/null
+++ b/src/utility/algorithminfo.jl
@@ -0,0 +1,281 @@
+"""
+$(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 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.
+
+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.
+
+## The vocabulary
+
+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.
+
+### Convergence
+
+  - `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.
+
+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
+    data::Dict{Symbol, Any}
+end
+
+"""
+    AlgorithmInfo(; truncation = nothing, kwargs...)
+
+Build an [`AlgorithmInfo`](@ref) from the entries an algorithm actually produced. Every keyword
+becomes an entry.
+
+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(; 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
+
+# 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}
+
+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`, 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)
+    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:")
+    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
+
+    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)
+    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/groundstate/approximate.jl b/test/groundstate/approximate.jl
index 31f5d1362..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)
@@ -119,9 +119,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 !haskey(info, :converged) && isnothing(convergence_measure(info))
         end
 
         @test norm(ψ - ψ_copy) < 1.0e-12
@@ -140,11 +141,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
 
@@ -156,11 +157,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)
@@ -173,26 +174,26 @@ end
         (pspace, Dspace, Dcut) in zipup_spacelist, left_to_right in (true, false)
         O, ψ = _random_mpo_mps(pspace, Dspace, 6)
         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/groundstate/groundstate.jl b/test/groundstate/groundstate.jl
index dd42b951b..c407199ba 100644
--- a/test/groundstate/groundstate.jl
+++ b/test/groundstate/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 convergence_measure(info) ≈ 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 convergence_measure(info) ≈ 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 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
@@ -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 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
@@ -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 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
@@ -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 convergence_measure(info) ≈ 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 convergence_measure(info) ≈ 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 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
@@ -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 convergence_measure(info) ≈ 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
 
     # Regression: IDMRG2 used to error on a non-abelian (SU2) unit cell due to a space mismatch.
@@ -291,7 +301,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)
@@ -303,7 +313,7 @@ end
         v = variance(ψ′, H, envs)
 
         # test using low variance
-        @test sum(δ) ≈ 0 atol = 1.0e-3
+        @test convergence_measure(info) ≈ 0 atol = 1.0e-3
         @test v < v₀
         @test v < 1.0e-2
     end
@@ -314,12 +324,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 convergence_measure(info) ≈ 0 atol = 1.0e-3
         @test v < v₀
         @test v < 1.0e-2
     end
@@ -349,13 +359,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
@@ -363,28 +373,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
@@ -411,26 +421,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
@@ -442,26 +452,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/symmetries/multifusion.jl b/test/symmetries/multifusion.jl
index 3096062a2..1631cca3a 100644
--- a/test/symmetries/multifusion.jl
+++ b/test/symmetries/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 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)
@@ -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 [convergence_measure(info), convergence_measure(info2), convergence_measure(info3)]
             @test delta ≈ 0 atol = 1.0e-3
         end
         for var in [v, v2, v3]
diff --git a/test/timeevolution/timestep.jl b/test/timeevolution/timestep.jl
index 124eb6f27..23b2a587c 100644
--- a/test/timeevolution/timestep.jl
+++ b/test/timeevolution/timestep.jl
@@ -263,6 +263,54 @@ 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
+        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.numiter == 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.ϵ_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 "aggregation over an evolution" begin
+        alg = TDVP2(; trunc = truncrank(2))
+        nsteps = 4
+        _, _, step = timestep(ψ₀, H, 0.0, dt, alg)
+        ψ, _, total = time_evolve(ψ₀, H, 0:dt:(nsteps * dt), alg)
+
+        @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
+        @test total.ϵ_max >= step.ϵ_max
+        @test total.ϵ_max <= total.ϵ_total
+    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))]