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 6d8f5ff89..b9437f025 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,16 @@ 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)) +- `find_groundstate(!)`, `leading_boundary`, `approximate(!)`, `timestep(!)` and `time_evolve(!)` + now have the return signature `(ψ, envs, info)` with `info` an `AlgorithmInfo` that contains + diagnostics about the run and convergence. See the updated docs or `AlgorithmInfo`'s + docstring for more information. ([#512](https://github.com/QuantumKitHub/MPSKit.jl/pull/512)) +- The manual has a new section on [Errors and accuracy](@ref) covering definitions and interpretations + of the various returned measures. ([#512](https://github.com/QuantumKitHub/MPSKit.jl/pull/512)) +- Convergence tolerances are now compared with `<=` everywhere, so an error exactly equal to `tol` + counts as converged for `DMRG`/`DMRG2` in `approximate`, the `IDMRG`/`IDMRG2` variants of + `approximate` and `leading_boundary`, and the gauge fixing of `InfiniteMPS`, as it already did for + the other algorithms. ([#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, @@ -132,11 +142,18 @@ When releasing a new version, move the "Unreleased" changes to a new version sec `TensorKit` 0.17.2 and `BlockTensorKit` <= 0.3.18, a permuted contraction into a sparse block tensor (e.g. an environment sweep near a `FiniteMPS` chain boundary when the operator and state have different `scalartype`s) could silently drop data instead of erroring. +- `approximate` with `IDMRG2` on a `MultilineMPS` with more than one row wrote the right-to-left + edge update into the wrong row and refreshed the unit cell edge of the first row only, so the + rows ended up mixed. ([#512](https://github.com/QuantumKitHub/MPSKit.jl/pull/512)) +- `leading_boundary` with `GradientGrassmann` ignored the algorithm's `hasconverged` and `shouldstop` + criteria, always using the optimizer's defaults instead. ([#512](https://github.com/QuantumKitHub/MPSKit.jl/pull/512)) ### Performance - 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), avoiding unnecessary re-orthogonalisation and complexification where applicable). ([#517](https://github.com/QuantumKitHub/MPSKit.jl/pull/517)) ## [0.13.11](https://github.com/QuantumKitHub/MPSKit.jl/compare/v0.13.10...v0.13.11) - 2026-05-04 diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md index 5ff913006..ddd518a48 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 ``` +See also [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,13 @@ WII TaylorCluster ``` +See also [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. +See also [Excitation accuracy](@ref). ```@docs; canonical=false excitations @@ -133,19 +143,16 @@ using TensorKit, MPSKit, MPSKitModels ### Quasiparticle Ansatz -The Quasiparticle Ansatz offers an approach to compute low-energy eigenstates in quantum -systems, playing a key role in both finite and infinite systems. It leverages localized -perturbations for approximations, as detailed in [haegeman2013](@cite). +The Quasiparticle Ansatz offers an approach to compute low-energy eigenstates in quantum systems, playing a key role in both finite and infinite systems. +It leverages localized perturbations for approximations, as detailed in [haegeman2013](@cite). #### Finite Systems: -In finite systems, we approximate low-energy states by altering a single tensor in the -Matrix Product State (MPS) for each site, and summing these across all sites. This method -introduces additional gauge freedoms, utilized to ensure orthogonality to the ground state. -Optimizing within this framework translates to solving an eigenvalue problem. For example, -in the transverse field Ising model, we calculate the first excited state as shown in the -provided code snippet, and check the accuracy against theoretical values. Some deviations -are expected, both due to finite-bond-dimension and finite-size effects. +In finite systems, we approximate low-energy states by altering a single tensor in the Matrix Product State (MPS) for each site, and summing these across all sites. +This method introduces additional gauge freedoms, utilized to ensure orthogonality to the ground state. +Optimizing within this framework translates to solving an eigenvalue problem. +For example, in the transverse field Ising model, we calculate the first excited state as shown in the provided code snippet, and check the accuracy against theoretical values. +Some deviations are expected, both due to finite-bond-dimension and finite-size effects. ```@example excitations # Model parameters @@ -164,10 +171,8 @@ isapprox(Es[1], 2(g - 1); rtol=1e-2) #### Infinite Systems: -The ansatz in infinite systems maintains translational invariance by perturbing every site -in the unit cell in a plane-wave superposition, requiring momentum specification. The -[Haldane gap](https://iopscience.iop.org/article/10.1088/0953-8984/1/19/001) computation in -the Heisenberg model illustrates this approach. +The ansatz in infinite systems maintains translational invariance by perturbing every site in the unit cell in a plane-wave superposition, requiring momentum specification. +The [Haldane gap](https://iopscience.iop.org/article/10.1088/0953-8984/1/19/001) computation in the Heisenberg model illustrates this approach. ```@example excitations # Setting up the model and momentum @@ -185,10 +190,9 @@ isapprox(Es[1], 0.41047925; atol=1e-4) #### Charged excitations: -When dealing with symmetric systems, the default optimization is for eigenvectors with -trivial total charge. However, quasiparticles with different charges can be obtained using -the sector keyword. For instance, in the transverse field Ising model, we consider an -excitation built up of flipping a single spin, aligning with `Z2Irrep(1)`. +When dealing with symmetric systems, the default optimization is for eigenvectors with trivial total charge. +However, quasiparticles with different charges can be obtained using the sector keyword. +For instance, in the transverse field Ising model, we consider an excitation built up of flipping a single spin, aligning with `Z2Irrep(1)`. ```@example excitations g = 10.0 @@ -206,8 +210,8 @@ QuasiparticleAnsatz ### Finite excitations -For finite systems we can also do something else - find the ground state of the Hamiltonian + -``\\text{weight} \sum_i | \\psi_i ⟩ ⟨ \\psi_i ``. This is also supported by calling +For finite systems we can instead find the ground state of the Hamiltonian where we project out previously obtained eigenstates, i.e. by projecting out ``\\text{weight} \sum_i | \\psi_i ⟩ ⟨ \\psi_i ``. +This is also supported by calling ```@example excitations # Model parameters @@ -229,17 +233,13 @@ FiniteExcited ### "Chepiga Ansatz" -Computing excitations in critical systems poses a significant challenge due to the diverging -correlation length, which requires very large bond dimensions. However, we can leverage this -long-range correlation to effectively identify excitations. In this context, the left/right -gauged MPS, serving as isometries, are effectively projecting the Hamiltonian into the -low-energy sector. This projection method is particularly effective in long-range systems, -where excitations are distributed throughout the entire system. Consequently, the low-lying -energy spectrum can be extracted by diagonalizing the effective Hamiltonian (without any -additional DMRG costs!). The states of these excitations are then represented by the ground -state MPS, with one site substituted by the corresponding eigenvector. This approach is -often referred to as the 'Chepiga ansatz', named after one of the authors of this paper -[chepiga2017](@cite). +Computing excitations in critical systems poses a significant challenge due to the diverging correlation length, which requires very large bond dimensions. +However, we can leverage this long-range correlation to effectively identify excitations. +In this context, the left/right gauged MPS, serving as isometries, are effectively projecting the Hamiltonian into the low-energy sector. +This projection method is particularly effective in long-range systems, where excitations are distributed throughout the entire system. +Consequently, the low-lying energy spectrum can be extracted by diagonalizing the effective Hamiltonian (without any additional DMRG costs!). +The states of these excitations are then represented by the ground state MPS, with one site substituted by the corresponding eigenvector. +This approach is often referred to as the 'Chepiga ansatz', named after one of the authors of this paper [chepiga2017](@cite). This is supported via the following syntax: @@ -254,14 +254,91 @@ Es, ϕs = excitations(H, ChepigaAnsatz(), ψ, envs; num=1) isapprox(Es[1] - E₀, 2(g - 1); rtol=1e-2) # infinite analytical result ``` -In order to improve the accuracy, a two-site version also exists, which varies two -neighbouring sites: +In order to improve the accuracy, a two-site version also exists, which varies two neighbouring sites: ```@example excitations 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. + +```@docs; canonical=false +AlgorithmInfo +``` + +The rest of this section explains what quantities can be reported by the algorithms, and how these can be interpreted. + +### Ground state accuracy + +[`find_groundstate`](@ref), [`leading_boundary`](@ref) and the iterative [`approximate`](@ref) algorithms stop once a convergence measure drops below their `tol`, and report that measure together with a `converged` flag. + +**Different measures for different algorithms.** +Which quantity is monitored depends on the algorithm, and these quantities cannot be compared to one another. +Each is therefore stored under its own name: `galerkin` for [`DMRG`](@ref), [`DMRG2`](@ref), [`VUMPS`](@ref) and [`VOMPS`](@ref), `gradientnorm` for [`GradientGrassmann`](@ref), `bondresidual` for [`IDMRG`](@ref) and [`IDMRG2`](@ref), and `localchange` for [`DMRG`](@ref) and [`DMRG2`](@ref) inside [`approximate`](@ref). +[`convergence_measure`](@ref) returns whichever of them is present, for code that only wants the number. +For the same reason, a `tol` that works well for one algorithm need not work well for another. + +**Interaction with a changing bond dimension.** +An algorithm that truncates changes the state by the discarded weight in every sweep, so its convergence measure cannot drop below the level set by what is discarded. +Rather than running into `maxiter` for a `tol` that cannot be reached, [`DMRG`](@ref) and [`DMRG2`](@ref) consider themselves converged once `galerkin ≤ max(tol, maximum(truncation_errors))`, which reduces to `galerkin ≤ tol` when nothing is truncated. +The reported `truncation_errors` are those of the final sweep, i.e. what the returned state still discards. + +**Convergence is not accuracy.** +For algorithms that remain on a fixed bond dimension MPS manifold, one can get the ground state *on that manifold* to converge to machine precision. +This however might still be far from the true ground state, whenever the bond dimension is not high enough. +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). +Likewise, none of these measures is an error bar on an observable. +An independent check is the energy [`variance`](@ref), ``\langle H^2 \rangle - \langle H \rangle^2``, which vanishes for an exact eigenstate, as is comparing results across bond dimensions. + +### Time evolution accuracy + +A time evolution has three distinct sources of error ([Paeckel et al.](@cite paeckel2019)): + +- The **truncation error**, from discarding singular values to keep the bond dimension bounded. + This is the only one that is reported, as the `truncation_errors` of [`timestep`](@ref) and [`time_evolve`](@ref). + It is non-zero for [`TDVP2`](@ref), for [`BUG`](@ref) with a `trunc`, and for [`TDVP`](@ref) with a bond expansion. +- The **projection error**, from restricting the evolution to MPS of the current bond dimension. + Measuring it costs an extra two-site effective-Hamiltonian application per bond, so it is not reported. + A fixed bond dimension turns entanglement growth into projection error, which is what a bond expansion (CBE) exists to reduce ([Li et al.](@cite li2024)). +- The **time-step error**, from splitting each step into a sequence of local updates, similar to a Trotter error. + It decreases with `dt`, and can be estimated by comparing one step of `dt` against two steps of `dt / 2`. + +Each of these has its own control: `dt` for the time-step error, `trunc` for the truncation error, and a bond-adaptive scheme ([`TDVP2`](@ref), [`BUG`](@ref), or [`TDVP`](@ref) with `alg_expand`) for the projection error whenever entanglement grows. + +These controls are not independent. +With a truncation strategy [`trunc::TruncationStrategy`](@extref MatrixAlgebraKit.TruncationStrategy), every step can discard weight, however small `dt` is, and reaching a fixed final time `T` takes `T / dt` steps. +Halving `dt` therefore lowers the time-step error, but doubles the number of truncations, so the accumulated truncation error can grow. +For a given `trunc` there is an optimal `dt`, below which taking smaller steps makes the result worse rather than better. + +### 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 +351,18 @@ state. changebonds ``` +All of these are controlled by a `trunc`, and the weight they discard is measured the same way as the `truncation_errors` of an [`AlgorithmInfo`](@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..92cd83616 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,46 @@ 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. + - `IDMRG2` fills `truncation_errors` with what its final sweep discarded at each bond, i.e. what + the returned state is still throwing away 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 + `truncation_errors` instead. + + See [`AlgorithmInfo`](@ref) for the full list and [Ground state accuracy](@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!` | @@ -93,13 +120,13 @@ function approximate( envs = environments(ψ, toapprox...) ) envs′ = Multiline([envs]) - multi, envs, δ = approximate( + multi, envs, info = approximate( convert(MultilineMPS, ψ), (convert(MultilineMPO, toapprox[1]), convert(MultilineMPS, toapprox[2])), algorithm, envs′ ) ψ = convert(InfiniteMPS, multi) - return ψ, envs, δ + return ψ, envs, _singleline_info(info) end # dispatch to in-place method @@ -115,11 +142,11 @@ function approximate( algorithm::Union{IDMRG, IDMRG2}, envs = environments(ψ, toapprox...) ) envs′ = Multiline([envs]) - multi, envs, δ = approximate( + multi, envs, info = approximate( convert(MultilineMPS, ψ), (convert(MultilineMPO, toapprox[1]), convert(MultilineMPS, toapprox[2])), algorithm, envs′ ) ψ = convert(InfiniteMPS, multi) - return ψ, envs, δ + return ψ, envs, _singleline_info(info) end diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl index 47dcfdfef..8de9325f5 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 + iter = 0 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) @@ -21,7 +22,7 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environment # finalize ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ), typeof(envs)} - if ϵ < alg.tol + if ϵ <= alg.tol @infov 2 logfinish!(log, iter, ϵ) break end @@ -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 + iter = 0 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) @@ -56,7 +58,7 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments # finalize ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ), typeof(envs)} - if ϵ < alg.tol + if ϵ <= alg.tol @infov 2 logfinish!(log, iter, ϵ) break end @@ -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..d9a8326f0 100644 --- a/src/algorithms/approximate/idmrg.jl +++ b/src/algorithms/approximate/idmrg.jl @@ -5,7 +5,7 @@ function approximate!( allocator = default_allocator(ψ, SerialScheduler()) log = IterLog("IDMRG") ϵ::Float64 = 2 * alg.tol - local iter + iter = 0 LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ) @@ -42,7 +42,7 @@ function approximate!( ϵ = norm(C_current - ψ.C[:, 0]) - if ϵ < alg.tol + if ϵ <= alg.tol @infov 2 logfinish!(log, iter, ϵ) break end @@ -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,7 +72,8 @@ function approximate!( ϵ::Float64 = 2 * alg.tol log = IterLog("IDMRG2") O, ϕ = toapprox - local iter + iter = 0 + ϵ_truncs = PeriodicMatrix(zeros(real(scalartype(ψ)), size(ψ))) LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ) @@ -86,7 +87,7 @@ 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, ϵ_truncs[row + 1, site] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) normalize!(c) ψ.AL[row + 1, site] = al @@ -100,14 +101,14 @@ function approximate!( end # update the edge - ψ.AL[1, end] = ψ.AC[1, end] / ψ.C[1, end] - ψ.AC[1, 1] = _mul_tail(ψ.AL[1, 1], ψ.C[1, 1]) + ψ.AL[:, end] .= ψ.AC[:, end] ./ ψ.C[:, end] + ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1]) for row in 1:size(ψ, 1) AC2′ = AC2_projection( 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, ϵ_truncs[row + 1, end] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) normalize!(c) ψ.AL[row + 1, end] = al @@ -132,7 +133,7 @@ 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, ϵ_truncs[row + 1, site] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) normalize!(c) ψ.AL[row + 1, site] = al @@ -145,22 +146,22 @@ function approximate!( end # update the edge - ψ.AC[1, end] = _mul_front(ψ.C[1, end - 1], ψ.AR[1, end]) - ψ.AR[1, 1] = _transpose_front(ψ.C[1, end] \ _transpose_tail(ψ.AC[1, 1])) + ψ.AC[:, end] .= _mul_front.(ψ.C[:, end - 1], ψ.AR[:, end]) + ψ.AR[:, 1] .= _transpose_front.(ψ.C[:, end] .\ _transpose_tail.(ψ.AC[:, 1])) for row in 1:size(ψ, 1) AC2′ = AC2_projection( 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, ϵ_truncs[row + 1, end] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd) normalize!(c) - ψ.AL[row, end] = al - ψ.C[row, end] = complex(c) - ψ.AR[row, 1] = _transpose_front(ar) + ψ.AL[row + 1, end] = al + ψ.C[row + 1, end] = complex(c) + ψ.AR[row + 1, 1] = _transpose_front(ar) - ψ.AR[row, end] = _transpose_front(ψ.C[row, end - 1] \ _transpose_tail(al * c)) - ψ.AC[row, 1] = _transpose_front(c * ar) + ψ.AR[row + 1, end] = _transpose_front(ψ.C[row + 1, end - 1] \ _transpose_tail(al * c)) + ψ.AC[row + 1, 1] = _transpose_front(c * ar) end transfer_leftenv!(envs, ψ, toapprox, 1) transfer_rightenv!(envs, ψ, toapprox, 0) @@ -175,7 +176,7 @@ function approximate!( return norm(e2' * c2 * e2 - e1' * c1 * e1) end - if ϵ < alg.tol + if ϵ <= alg.tol @infov 2 logfinish!(log, iter, ϵ) break end @@ -193,5 +194,6 @@ function approximate!( copy!(ψ, ψ′) # ensure output destination is unchanged recalculate!(envs, ψ, toapprox) - return ψ, envs, ϵ + info = AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, truncation_errors = parent(ϵ_truncs), 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..5d6609f46 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 @@ -73,11 +73,12 @@ function approximate!(ψ::FiniteMPS, (O, ϕ)::Tuple{Any, <:FiniteMPS}, alg::Zipu throw(SpaceMismatch("MPO input physical space does not match MPS physical space at site $i")) end - return if alg.left_to_right + ψ, ϵ_truncs = if alg.left_to_right zip_left_right!(ψ, O, ϕ, alg.alg_zipup, alg.alg_zipdown) else zip_right_left!(ψ, O, ϕ, alg.alg_zipup, alg.alg_zipdown) end + return ψ, AlgorithmInfo(; truncation_errors = ϵ_truncs) end function approximate(Oϕ::Tuple{Any, <:FiniteMPS}, alg::Zipup) @@ -96,8 +97,8 @@ 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 the truncation error `ϵ` of every bond, the 2-norm of the singular values discarded +by the last cut made there. """ zip_left_right! @doc (@doc zip_left_right!) zip_right_left! @@ -115,14 +116,13 @@ 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(ψ))) + ϵ_truncs = zeros(real(scalartype(ψ)), N - 1) # 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 + AL, Fₗ, ϵ_truncs[i] = left_gauge(Aᶻ, alg_zipup) # right factor doubles as the next left fuser ψ.ALs[i] = AL - ϵ = max(ϵ, ϵᵢ) 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ᵣ) @@ -130,12 +130,11 @@ function zip_left_right!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown # zip down from right to left, truncating in a locally gauged basis if !isnothing(alg_zipdown) for i in N:-1:2 - ψ, ϵᵢ = right_gauge!(ψ, i, ψ.AC[i], alg_zipdown) - ϵ = max(ϵ, ϵᵢ) + ψ, ϵ_truncs[i - 1] = right_gauge!(ψ, i, ψ.AC[i], alg_zipdown) end end - return ψ, ϵ + return ψ, ϵ_truncs end function zip_right_left!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown = nothing) @@ -149,14 +148,13 @@ 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(ψ))) + ϵ_truncs = zeros(real(scalartype(ψ)), N - 1) # 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 + Fᵣ, AR, ϵ_truncs[i - 1] = _right_gauge_zip(Aᶻ, alg_zipup) # left factor doubles as the next right fuser ψ.ARs[i] = AR - ϵ = max(ϵ, ϵᵢ) 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)) @@ -165,12 +163,11 @@ function zip_right_left!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown # zip down from left to right, truncating in a locally gauged basis if !isnothing(alg_zipdown) for i in 1:(N - 1) - ψ, ϵᵢ = left_gauge!(ψ, i, ψ.AC[i], alg_zipdown) - ϵ = max(ϵ, ϵᵢ) + ψ, ϵ_truncs[i] = left_gauge!(ψ, i, ψ.AC[i], alg_zipdown) end end - return ψ, ϵ + return ψ, ϵ_truncs end # `right_gauge` for the tensors of a right-to-left zip-up sweep: these are already partitioned across diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index 0afe339df..4db4ac557 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -81,9 +81,9 @@ function changebonds!(mpo::FiniteMPO, alg::SvdCut) O = transpose(mpo[N], ((1,), (3, 4, 2))) for i in (N - 1):-1:1 U, S, Vᴴ = svd_trunc!(O; trunc = alg.trunc, alg = alg.alg_svd) - _warn_empty_bond(i, space(S, 1)) n = sqrt(norm(S)) - iszero(n) && (n = one(n)) + _warn_empty_bond(i, space(S, 1), n) + n <= sqrt(eps(real(one(n)))) && (n = one(n)) @inbounds mpo[i + 1] = transpose(scale!(Vᴴ, n), ((1, 4), (2, 3))) if i > 1 @plansor O[-1; -3 -4 -2] := mpo[i][-1 -2; -3 2] * U[2; 1] * S[1; -4] / n @@ -103,8 +103,12 @@ function _extract_norm!(t) return log(n) end -function _warn_empty_bond(bond::Int, V) - dim(V) == 0 && @warn "`SvdCut` truncated the bond between sites $bond and $(bond + 1) down to zero dimensions; the resulting operator is identically zero. Loosen `trunc` or check the scale of the input operator." +function _warn_empty_bond(bond::Int, V, n) + if dim(V) == 0 + @warn "`SvdCut` truncated the bond between sites $bond and $(bond + 1) down to zero dimensions; the resulting operator is identically zero. Loosen `trunc` or check the scale of the input operator." + elseif n <= sqrt(eps(real(one(n)))) + @warn "`SvdCut` truncated the bond between sites $bond and $(bond + 1) down to a near-zero norm ($n); the rescaled tensor may lose numerical precision. Loosen `trunc` or check the scale of the input operator." + end return nothing end diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl index 1fb1654f6..46146377b 100644 --- a/src/algorithms/derivatives/mpo_derivatives.jl +++ b/src/algorithms/derivatives/mpo_derivatives.jl @@ -221,16 +221,8 @@ function prepare_operator!!(H::MPO_C_Hamiltonian{<:MPSTensor, <:MPSTensor}) return prepared_operator_type(typeof(H))(leftenv, rightenv, H.backend, H.allocator) end -# Scratch space of `prepare_operator!!` -# ------------------------------------ -# The environments of a prepared derivative are the dense, fused form of an environment-operator -# contraction. Neither that contraction nor its dense copy is kept - only the repartitioned -# result is - so both are taken from the allocator and handed straight back. -# -# `GL_O` and `O_GR` are allocated here rather than by `:=`, which always allocates the result of -# a `@plansor` block on the heap. - -# A `TensorMap` is already dense, so `repartition` is the only tensor that has to be kept. +# Only the repartitioned result of `prepare_operator!!`'s environment-operator contraction is +# kept, so the contraction itself and its dense copy are taken from the allocator and freed here. @inline function _fuse_env(t::TensorMap, N₁::Int, N₂::Int, backend, allocator) return repartition(fuse_legs(t, N₁, N₂), 2, 2; copy = true, backend, allocator) end diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index b6aa8ef70..7b5bb3491 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -60,7 +60,7 @@ $(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, see [Ground state accuracy](@ref)" tol::Float64 "maximal amount of iterations" @@ -171,7 +171,7 @@ $(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, see [Ground state accuracy](@ref)" tol::Float64 "maximal amount of iterations" @@ -232,7 +232,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) @@ -253,10 +253,15 @@ _num_updates(::DMRG2, ψ) = length(ψ) - 1 _sweep_ranges(::DMRG, ψ) = (1:(length(ψ) - 1), length(ψ):-1:2) _sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1) +# per-bond truncation errors of the last cut from the per-update slots: single-site slot `pos` +# holds the backward cut of bond `pos - 1`, which is the last one made there +_bond_truncation_errors(::DMRG, ϵ_truncs) = ϵ_truncs[2:end] +_bond_truncation_errors(::DMRG2, ϵ_truncs) = copy(ϵ_truncs) + 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 +278,10 @@ Currently supported for the finite-system algorithms [`DMRG`](@ref) and [`DMRG2` - `ψ::AbstractFiniteMPS`: converged ground state - `environments`: environments corresponding to the converged state -- `ϵ::Float64`: final convergence error upon terminating the algorithm +- `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. `info.truncation_errors` holds, per bond, what the last cut there discarded + (all zero for a gauge that does not truncate). See [`AlgorithmInfo`](@ref). """ function find_groundstate!( ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs = environments(ψ, H, ψ) @@ -300,10 +308,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, ψ) + iter = 0 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 +361,12 @@ function find_groundstate_sweep!( end end end - return ψ, envs, ϵ_global + + info = AlgorithmInfo(; + converged = ϵ_global <= max(alg.tol, maximum(ϵ_truncs)), galerkin = ϵ_global, + truncation_errors = _bond_truncation_errors(alg, ϵ_truncs), 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..dc52266a9 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.truncation_errors` with what its final sweep discarded at each bond. + See [`AlgorithmInfo`](@ref) for the full vocabulary, and [Ground state accuracy](@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..365c20bee 100644 --- a/src/algorithms/groundstate/gradient_grassmann.jl +++ b/src/algorithms/groundstate/gradient_grassmann.jl @@ -14,7 +14,7 @@ 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 on the norm of the gradient, see [Ground state accuracy](@ref) - `maxiter = Defaults.maxiter`: maximum amount of iterations - `verbosity = Defaults.verbosity - 1`: level of information display - `hasconverged = OptimKit.DefaultHasConverged(tol)`: convergence criterium @@ -79,7 +79,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!(ψ) @@ -109,7 +109,7 @@ function _find_groundstate(ψ, H, alg::GradientGrassmann, envs, scheduler, timer () -> GrassmannMPS.precondition(state, g), timeroutput, "precondition", ) - x, _, _, _, normgradhistory = optimize( + x, f, g, _, normgradhistory = optimize( fg, ψ, alg.method; retract, transport!, precondition, GrassmannMPS.inner, @@ -125,5 +125,17 @@ function _find_groundstate(ψ, H, alg::GradientGrassmann, envs, scheduler, timer @infov 4 TimerReport(timeroutput) end - return x, envs, normgradhistory[end] + info = _optimkit_info(alg, x, f, g, normgradhistory) + return x, envs, info +end + +# `optimize` returns its full history as `[fhistory normgradhistory]`, with the initial point +# before the first iteration as its first row. Convergence is judged by the same `hasconverged` +# the optimizer stopped on. +function _optimkit_info(alg::GradientGrassmann, x, f, g, normgradhistory) + normres = normgradhistory[end, 2] + return AlgorithmInfo(; + converged = alg.hasconverged(x, f, g, normres), gradientnorm = normres, + numiter = size(normgradhistory, 1) - 1 + ) end diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 519933630..b6da37a17 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -12,7 +12,7 @@ $(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 on the change of the center bond tensor over a sweep, see [Ground state accuracy](@ref)" tol::Float64 = Defaults.tol "maximal amount of iterations" @@ -45,7 +45,7 @@ $(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 on the change of the center bond tensor over a sweep, see [Ground state accuracy](@ref)" tol::Float64 = Defaults.tol "maximal amount of iterations" @@ -71,22 +71,24 @@ Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_bounda end # Internal state of the IDMRG algorithm -struct IDMRGState{S, O, E, T, TO, A} +struct IDMRGState{S, O, E, V, T, TO, A} mps::S operator::O envs::E iter::Int ϵ::Float64 # TODO: Could be any <:Real + truncation_errors::V # per bond, 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_errors::V, 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 + ) where {S, O, E, V, T, TO, A} + return IDMRGState{S, O, E, V, T, TO, A}( + mps, operator, envs, iter, ϵ, truncation_errors, T(energy), timeroutput, allocator ) end @@ -113,7 +115,8 @@ function _find_groundstate_idmrg(mps, operator, alg::alg_type, envs) where {alg_ end end - state = IDMRGState(mps, operator, envs, iter, ϵ, E, timeroutput, allocator) + ϵ_truncs = zeros(real(scalartype(mps)), length(mps)) + state = IDMRGState(mps, operator, envs, iter, ϵ, ϵ_truncs, E, timeroutput, allocator) it = IterativeSolver(alg, state) return LoggingExtras.withlevel(; alg.verbosity) do @@ -134,15 +137,21 @@ 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, ψ′; alg.backend) - return ψ′, envs, it.state.ϵ + info = AlgorithmInfo(; + converged = it.state.ϵ <= alg.tol, bondresidual = it.state.ϵ, + truncation_errors = alg isa IDMRG2 ? it.state.truncation_errors : nothing, + numiter = it.state.iter + ) + return ψ′, envs, info end end function Base.iterate( - it::IterativeSolver{alg_type}, state::IDMRGState{<:Any, <:Any, <:Any, T} = it.state + it::IterativeSolver{alg_type}, state::IDMRGState{<:Any, <: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) + ϵ_truncs = zero(state.truncation_errors) # fresh each sweep, filled by the sweep itself + mps, envs, C_old, E_new = @timeit timeroutput "localupdate" localupdate_step!(it, state, ϵ_truncs) # error criterion C = mps.C[0] @@ -163,14 +172,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, ϵ, ϵ_truncs, E_new, timeroutput, state.allocator, ) return (mps, envs, ϵ, ΔE), it.state end function localupdate_step!( - it::IterativeSolver{<:IDMRG}, state + it::IterativeSolver{<:IDMRG}, state, ϵ_truncs ) alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = state.iter, g_global = state.ϵ) return _localupdate_sweep_idmrg!( @@ -180,12 +189,12 @@ function localupdate_step!( end function localupdate_step!( - it::IterativeSolver{<:IDMRG2}, state + it::IterativeSolver{<:IDMRG2}, state, ϵ_truncs ) 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, ϵ_truncs; it.backend, state.allocator, ) end @@ -229,7 +238,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, ϵ_truncs::AbstractVector; backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator() ) # @timeit wraps its body in try-finally, which is a new lexical scope: declare locals @@ -243,7 +252,7 @@ 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, ϵ_truncs[pos] = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) normalize!(c) ψ.AL[pos] = al @@ -266,7 +275,7 @@ 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, ϵ_truncs[end] = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) normalize!(c) ψ.AL[end] = al @@ -294,7 +303,7 @@ 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, ϵ_truncs[pos] = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) normalize!(c) ψ.AL[pos] = al @@ -318,7 +327,7 @@ 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, ϵ_truncs[end] = svd_trunc!(ac2′; trunc = alg_trunc, alg = alg_svd) normalize!(c) ψ.AL[end] = al diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index aa51f44b5..f766edc0d 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -17,7 +17,7 @@ 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 on the Galerkin error, see [Ground state accuracy](@ref)" tol::Float64 = Defaults.tol "maximal amount of iterations" @@ -85,12 +85,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..347411f57 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -28,7 +28,7 @@ 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 on the largest change of a center tensor over a sweep" 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..6ad788fe4 100644 --- a/src/algorithms/statmech/gradient_grassmann.jl +++ b/src/algorithms/statmech/gradient_grassmann.jl @@ -7,7 +7,7 @@ function leading_boundary( # read the scheduler here rather than in `fg`, so that the allocator it selects is inferable scheduler = Defaults.scheduler[] fg(x) = GrassmannMPS.fg(x, operator, envs; alg.backend, scheduler) - x, _, _, _, normgradhistory = optimize( + x, f, g, _, normgradhistory = optimize( fg, state, alg.method; GrassmannMPS.transport!, @@ -17,7 +17,11 @@ function leading_boundary( GrassmannMPS.add!, GrassmannMPS.precondition, alg.finalize!, + alg.hasconverged, + alg.shouldstop, isometrictransport = true ) - return x, envs, normgradhistory[end] + + info = _optimkit_info(alg, x, f, g, normgradhistory) + return x, envs, info end diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 84edab763..29e7b3b8a 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -4,7 +4,7 @@ function leading_boundary( allocator = default_allocator(ψ, SerialScheduler()) log = IterLog("IDMRG") ϵ::Float64 = 2 * alg.tol - local iter + iter = 0 LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ, expectation_value(ψ, operator, envs)) @@ -43,7 +43,7 @@ function leading_boundary( ϵ = norm(C_current - ψ.C[:, 0]) - if ϵ < alg.tol + if ϵ <= alg.tol @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, operator, envs)) break end @@ -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,7 +69,8 @@ function leading_boundary( size(ψ, 2) < 2 && throw(ArgumentError("unit cell should be >= 2")) ϵ::Float64 = 2 * alg.tol log = IterLog("IDMRG2") - local iter + iter = 0 + ϵ_truncs = PeriodicMatrix(zeros(real(scalartype(ψ)), size(ψ))) LoggingExtras.withlevel(; alg.verbosity) do @infov 2 loginit!(log, ϵ) @@ -84,7 +85,7 @@ 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, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) normalize!(c) ψ.AL[row + 1, site] = al @@ -108,7 +109,7 @@ 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, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) normalize!(c) ψ.AL[row + 1, site] = al @@ -133,7 +134,7 @@ 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, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) normalize!(c) ψ.AL[row + 1, site] = al @@ -156,7 +157,7 @@ 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, ϵ_truncs[row + 1, end] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd) normalize!(c) ψ.AL[row + 1, end] = al @@ -180,7 +181,7 @@ function leading_boundary( return norm(e2' * c2 * e2 - e1' * c1 * e1) end - if ϵ < alg.tol + if ϵ <= alg.tol @infov 2 logfinish!(log, iter, ϵ) break end @@ -196,5 +197,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_errors = parent(ϵ_truncs), numiter = iter) end diff --git a/src/algorithms/statmech/leading_boundary.jl b/src/algorithms/statmech/leading_boundary.jl index 218a3f8f7..c9796af8d 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,21 @@ 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. [`IDMRG2`](@ref) additionally reports `truncation_errors`, what its final sweep + discarded at each bond. + + See [`AlgorithmInfo`](@ref), [`find_groundstate`](@ref), and the manual on + [Ground state accuracy](@ref). """ leading_boundary # TODO: alg selector @@ -31,20 +46,20 @@ optimization algorithm will be attempted based on the supplied keywords. function leading_boundary(state::InfiniteMPS, operator::InfiniteMPO, alg) state_multi = convert(MultilineMPS, state) operator_multi = convert(MultilineMPO, operator) - state_multi′, envs_multi′, err = leading_boundary( + state_multi′, envs_multi′, info = leading_boundary( state_multi, operator_multi, alg, ) state′ = convert(InfiniteMPS, state_multi′) - return state′, only(envs_multi′), err + return state′, only(envs_multi′), _singleline_info(info) end function leading_boundary(state::InfiniteMPS, operator::InfiniteMPO, alg, envs) state_multi = convert(MultilineMPS, state) operator_multi = convert(MultilineMPO, operator) envs_multi = Multiline([envs]) - state_multi′, envs_multi′, err = leading_boundary( + state_multi′, envs_multi′, info = leading_boundary( state_multi, operator_multi, alg, envs_multi ) state′ = convert(InfiniteMPS, state_multi′) - return state′, only(envs_multi′), err + return state′, only(envs_multi′), _singleline_info(info) end diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl index 9427f5488..87c5cb8cb 100644 --- a/src/algorithms/statmech/vomps.jl +++ b/src/algorithms/statmech/vomps.jl @@ -18,7 +18,7 @@ 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 on the Galerkin error, see [Ground state accuracy](@ref)" tol::Float64 = Defaults.tol "maximal amount of iterations" @@ -73,11 +73,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..24df2c962 100644 --- a/src/algorithms/timestep/bug.jl +++ b/src/algorithms/timestep/bug.jl @@ -23,10 +23,15 @@ 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 all bond cuts of the step, + whereas in imaginary time the norm also carries the physical decay of the weight. Pass `normalize = true` to `timestep`/`time_evolve` to renormalize after every half-sweep instead. +!!! warning + The `truncation_errors` reported by [`timestep`](@ref) are those of the cuts in step (i). Step (iii) + augments the basis without truncating, so within a sweep the bond dimension temporarily exceeds `trunc`. + !!! 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 +144,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 +159,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 +177,27 @@ function timestep!( # the sweep is serial, so a single allocator serves all local updates allocator = default_allocator(ψ, SerialScheduler()) + ϵ_truncs = zeros(real(scalartype(ψ)), L - 1) + # 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 ) + site < L && (ϵ_truncs[site] = ϵ) 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 ) + site > 1 && (ϵ_truncs[site - 1] = ϵ) end - return ψ, envs + return ψ, envs, AlgorithmInfo(; truncation_errors = ϵ_truncs) end # copying version diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index a3da18c52..af4707f2d 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -11,11 +11,14 @@ 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 the squared norm drops by precisely the + weight discarded by all cuts of the step, whereas in imaginary time the norm also carries the + physical decay of the weight. 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 +132,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_errors` (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 +154,8 @@ function _timestep_finite!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs, allocator; imaginary_evolution::Bool, normalize::Bool ) + ϵ_truncs = zeros(real(scalartype(ψ)), length(ψ) - 1) + # sweep left to right for i in 1:(length(ψ) - 1) # 1. optionally expand the bond ahead of the local update (CBE) @@ -161,7 +169,7 @@ 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) + _, ϵ_truncs[i] = left_gauge!(ψ, i, AC, alg.alg_gauge; normalize) # 4. evolve the bond tensor backward Hc = C_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator) @@ -190,7 +198,7 @@ function _timestep_finite!( # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 (norm preserved by # default; `normalize` renormalizes) - right_gauge!(ψ, i, AC, alg.alg_gauge; normalize) + _, ϵ_truncs[i - 1] = right_gauge!(ψ, i, AC, alg.alg_gauge; normalize) # 4. evolve the bond tensor backward Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator) @@ -207,13 +215,14 @@ function _timestep_finite!( imaginary_evolution ) - return ψ, envs + return ψ, envs, AlgorithmInfo(; truncation_errors = ϵ_truncs) end """ $(TYPEDEF) Two-site MPS time-evolution algorithm based on the Time-Dependent Variational Principle. +See [`TDVP`](@ref) for more information. # Fields @@ -266,16 +275,19 @@ 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) + + ϵ_truncs = zeros(real(scalartype(ψ)), length(ψ) - 1) + # 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 + _, ϵ_truncs[i] = gauge2!(ψ, i, Val(:right), ac2′, alg_gauge; normalize) if i != (length(ψ) - 1) Hac = AC_hamiltonian(i + 1, ψ, H, ψ, envs; alg.backend, allocator) @@ -292,10 +304,7 @@ 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)) + _, ϵ_truncs[i - 1] = gauge2!(ψ, i - 1, Val(:left), ac2′, alg_gauge; normalize) if i != 2 Hac = AC_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator) @@ -306,7 +315,7 @@ function _timestep2_finite!( end end - return ψ, envs + return ψ, envs, AlgorithmInfo(; truncation_errors = ϵ_truncs) end # copying version diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 68014ec7f..0214b8e88 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,13 @@ through each of the time points obtained by iterating t_span. - `ψ`: the time-stepped state - `envs`: the updated environment manager +- `info::AlgorithmInfo`: `numiter`, the number of steps taken, and for an algorithm that + truncates `truncation_errors`, holding the per-bond truncation errors of every step as reported + by [`timestep`](@ref), one vector per step. See [`AlgorithmInfo`](@ref) and + [Time evolution accuracy](@ref) in the manual for what this does not measure. + +The largest truncation error of each step is logged at `verbosity ≥ 3`, and that of the whole +evolution at `verbosity ≥ 2`. """ function time_evolve end, function time_evolve! end @@ -35,27 +42,42 @@ 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)))) + truncation_errors = [] + ϵ_max = 0.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)} - @infov 3 logiter!(log, iter, 0.0, t) + # the log shows the largest per-bond error, or zero for a non-truncating algorithm + ϵ_step = 0.0 + if haskey(info_step, :truncation_errors) + push!(truncation_errors, info_step.truncation_errors) + ϵ_step = Float64(maximum(info_step.truncation_errors; init = 0.0)) + end + ϵ_max = max(ϵ_max, ϵ_step) + @infov 3 logiter!(log, iter, ϵ_step, t) end - @infov 2 logfinish!(log, length(t_span), 0.0, t_span[end]) + @infov 2 logfinish!(log, length(t_span), ϵ_max, t_span[end]) end - return ψ, envs + info = AlgorithmInfo(; + numiter = length(t_span) - 1, + truncation_errors = isempty(truncation_errors) ? nothing : copy(truncation_errors) + ) + 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 +103,20 @@ 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 finite-system step sweeps over every bond twice, and `info.truncation_errors[i]` is the +truncation error of the last cut made at bond `i`, i.e. of the second half-sweep. +The entries are non-zero only for algorithms that truncate ([`TDVP2`](@ref), [`BUG`](@ref) with a +`trunc`, and [`TDVP`](@ref) with a bond expansion), and exactly `0` for a step that happened to +discard nothing. Infinite one-site [`TDVP`](@ref) never truncates and reports no +`truncation_errors` 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/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 44a9fbb79..317750329 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -22,7 +22,7 @@ function entropy(spectrum::TensorKit.SectorVector{T}) where {T} for (c, b) in pairs(spectrum) s = sum(b; init = zero(S)) do x x² = x^2 - return x < tol ? zero(x) : x² * log(x²) + return x <= tol ? zero(x) : x² * log(x²) end S += oftype(S, dim(c) * s) end diff --git a/src/states/ortho.jl b/src/states/ortho.jl index 515f899f2..5a456f120 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 @@ -238,7 +238,7 @@ function uniform_leftorth!( # iteratively solve for (AL, C) in it iter, ϵ = it.iter, it.ϵ - if ϵ < it.tol + if ϵ <= it.tol @infov 2 logfinish!(log, iter, ϵ) return AL, C elseif iter > it.maxiter @@ -306,7 +306,7 @@ function uniform_rightorth!( # iteratively solve for (AR, C) in it iter, ϵ = it.iter, it.ϵ - if ϵ < it.tol + if ϵ <= it.tol @infov 2 logfinish!(log, iter, ϵ) return AR, C elseif iter > it.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..1da88d34e --- /dev/null +++ b/src/utility/algorithminfo.jl @@ -0,0 +1,172 @@ +""" +$(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 + + - `truncation_errors`: the truncation error of every bond, i.e. the 2-norm of the singular values + discarded by the most recent factorisation of that bond. Entry `i` is the bond between sites + `i` and `i + 1`, so a finite MPS of length `L` has `L - 1` entries and an infinite unit cell of + length `L` has `L`, the last one being the bond across the unit cell boundary. A `Multiline` + algorithm reports a matrix instead, indexed as `[row, bond]`. + +A sweep that cuts a bond more than once keeps only the last cut, which is what the returned state +still throws away there. An algorithm that truncates reports this entry even when it discarded +nothing (all entries zero), whereas the entry being absent 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(; 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. + +`numiter` defaults to `1` for the single-shot algorithms. +""" +function AlgorithmInfo(; kwargs...) + data = Dict{Symbol, Any}() + for (key, value) in kwargs + isnothing(value) || (data[key] = value) + end + get!(data, :numiter, 1) + 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 + +# the single-line wrappers of the `Multiline` algorithms report a vector instead of a 1-row matrix +function _singleline_info(info::AlgorithmInfo) + data = copy(getfield(info, :data)) + haskey(data, :truncation_errors) && (data[:truncation_errors] = vec(data[:truncation_errors])) + return AlgorithmInfo(data) +end + +# dictionary interface +Base.getindex(info::AlgorithmInfo, key::Symbol) = getfield(info, :data)[key] +Base.haskey(info::AlgorithmInfo, key::Symbol) = haskey(getfield(info, :data), key) +Base.get(info::AlgorithmInfo, key::Symbol, default) = get(getfield(info, :data), key, default) +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) + haskey(data, key) && return data[key] + return _no_entry_error(info, key) +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 + +# custom show +# convergence measures are displayed first, followed by everything else in alphabetical order +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", 18), " = ", data[:converged], + " after ", numiter, " iterations" + ) + elseif !isnothing(numiter) + println(io, tab_space, numiter, " iteration", numiter == 1 ? "" : "s") + end + + handled = (:converged, :numiter, convergence_keys...) + others = sort!(filter(∉(handled), collect(keys(data)))) + compact = IOContext(io, :compact => true, :limit => true) + for key in (filter(in(keys(data)), convergence_keys)..., others...) + print(io, tab_space, rpad(string(key), 18), " = ") + show(compact, data[key]) + println(io) + 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/algorithms/anyonic_braiding.jl b/test/algorithms/anyonic_braiding.jl index 22f59773b..8ea2d1121 100644 --- a/test/algorithms/anyonic_braiding.jl +++ b/test/algorithms/anyonic_braiding.jl @@ -75,8 +75,8 @@ end # even though the eigenproblem is degenerate over the ancilla leg A = randn(T, V^L, V^L) Oh = FiniteMPO(A + A') - ψ3, _, ϵ3 = find_groundstate(ρ_mps, Oh, DMRG(; tol = 1.0e-10, maxiter = 50)) - @test ϵ3 ≤ 1.0e-10 + ψ3, _, info3 = find_groundstate(ρ_mps, Oh, DMRG(; tol = 1.0e-10, maxiter = 50)) + @test info3.galerkin ≤ 1.0e-10 # prepared AC2, reached by two-site TDVP Ohd = convert(TensorMap, Oh) diff --git a/test/groundstate/approximate.jl b/test/groundstate/approximate.jl index 31f5d1362..3b1ae7f8f 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) @@ -53,10 +53,14 @@ end end ψ3, _ = approximate(ψ0, (W1, ψ), IDMRG(; verbosity)) - ψ4, _ = approximate(ψ0, (sW2, ψ), IDMRG2(; trunc = truncrank(12), verbosity)) - ψ5, _ = timestep(ψ, H, 0.0, dt, TDVP()) + ψ4, _, info4 = approximate(ψ0, (sW2, ψ), IDMRG2(; trunc = truncrank(12), verbosity)) + ψ5, _, info5 = timestep(ψ, H, 0.0, dt, TDVP()) ψ6 = changebonds(W1 * ψ, SvdCut(; trunc = truncrank(12))) + @test info4.truncation_errors isa AbstractVector + @test length(info4.truncation_errors) == length(ψ0) + @test !haskey(info5, :truncation_errors) # infinite one-site TDVP never truncates + @test abs(dot(ψ1, ψ5)) ≈ 1.0 atol = dt @test abs(dot(ψ3, ψ5)) ≈ 1.0 atol = dt @test abs(dot(ψ6, ψ5)) ≈ 1.0 atol = dt @@ -66,6 +70,28 @@ end @test dim(space(nW1[1], 1)) == 1 end + # every row of a multiline problem has to land in its own row of the result, including at the + # unit cell edge, which used to write the last row's update into the wrong row + @testset "multiline IDMRG2 rows" begin + Random.seed!(42) + dt = 1.0e-2 + Wa = make_time_mpo(force_planar(repeat(transverse_field_ising(; g = 4), 2)), dt, WII()) + Wb = make_time_mpo(force_planar(repeat(transverse_field_ising(; g = 0.5), 2)), dt, WII()) + ϕa, ϕb = (InfiniteMPS([ℙ^2, ℙ^2], [ℙ^8, ℙ^8]) for _ in 1:2) + ψ0a, ψ0b = (InfiniteMPS([ℙ^2, ℙ^2], [ℙ^10, ℙ^10]) for _ in 1:2) + alg = IDMRG2(; trunc = truncrank(10), verbosity, maxiter = 50) + + multi, = approximate( + MultilineMPS([ψ0a, ψ0b]), (MultilineMPO([Wa, Wb]), MultilineMPS([ϕa, ϕb])), alg + ) + ref_a, = approximate(ψ0a, (Wa, ϕa), alg) + ref_b, = approximate(ψ0b, (Wb, ϕb), alg) + + # row `r + 1` of the result approximates `O[r] * ϕ[r]` + @test abs(dot(multi[2], ref_a)) ≈ 1 atol = 1.0e-8 + @test abs(dot(multi[1], ref_b)) ≈ 1 atol = 1.0e-8 + end + finite_algs = [DMRG(; verbosity), DMRG2(; verbosity, trunc = truncrank(10))] @testset "finitemps1 ≈ finitemps2" for alg in finite_algs a = FiniteMPS(10, ℂ^2, ℂ^10) @@ -119,9 +145,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, ϵ = 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 length(info.truncation_errors) == length(ψ) - 1 + @test maximum(info.truncation_errors) < 1.0e-10 + @test !haskey(info, :converged) && isnothing(convergence_measure(info)) end @test norm(ψ - ψ_copy) < 1.0e-12 @@ -140,11 +168,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 maximum(info.truncation_errors) < 1.0e-10 end end @@ -156,11 +184,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 maximum(info_one_sweep.truncation_errors) > 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 +201,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.truncation_errors ≈ info_ref.truncation_errors # 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.truncation_errors ≈ info_ref.truncation_errors # 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.truncation_errors ≈ info_ref.truncation_errors 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 527947e05..58e349692 100644 --- a/test/groundstate/groundstate.jl +++ b/test/groundstate/groundstate.jl @@ -31,19 +31,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.truncation_errors == zeros(L - 1) end @testset "DMRG2" begin @@ -51,19 +52,22 @@ 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 length(info.truncation_errors) == L - 1 + @test maximum(info.truncation_errors) > 0 end @testset "CBEDMRG" begin @@ -74,17 +78,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 @@ -102,17 +106,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 @@ -128,17 +132,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 @@ -150,17 +154,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(; tol, 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 @@ -172,13 +176,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), @@ -231,7 +235,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) @@ -241,7 +245,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 @@ -251,11 +255,11 @@ end ψ₀ = InfiniteMPS(fill(ℙ^2, 3), fill(ℙ^D, 3)) v₀ = variance(ψ₀, H) - ψ′, envs, δ = find_groundstate(ψ₀, H, VUMPS(; tol, verbosity = verbosity_conv, maxiter = 20)) + ψ′, envs, info = find_groundstate(ψ₀, H, VUMPS(; tol, verbosity = verbosity_conv, maxiter = 20)) 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 @@ -265,15 +269,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, :truncation_errors) # single-site IDMRG never truncates @test v < v₀ @test v < 1.0e-2 end @@ -283,11 +288,11 @@ end ψ₀ = InfiniteMPS(fill(ℙ^2, 3), fill(ℙ^D, 3)) v₀ = variance(ψ₀, H) - ψ, envs, δ = find_groundstate(ψ₀, H, IDMRG(; tol, verbosity = verbosity_conv, maxiter = 20)) + ψ, envs, info = find_groundstate(ψ₀, H, IDMRG(; tol, verbosity = verbosity_conv, maxiter = 20)) 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 @@ -299,19 +304,22 @@ 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 length(info.truncation_errors) == length(ψ) + @test maximum(info.truncation_errors) > 0 end # Regression: IDMRG2 used to error on a non-abelian (SU2) unit cell due to a space mismatch. @@ -336,7 +344,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) @@ -348,7 +356,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 @@ -359,12 +367,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 @@ -394,13 +402,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 @@ -408,28 +416,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 @@ -456,26 +464,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 @@ -487,26 +495,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 d050381d7..0e2d9b09f 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..3baebcf91 100644 --- a/test/timeevolution/timestep.jl +++ b/test/timeevolution/timestep.jl @@ -263,6 +263,49 @@ 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 errors are 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.truncation_errors == zeros(L - 1) + @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 length(tight.truncation_errors) == L - 1 + @test maximum(tight.truncation_errors) > 1.0e-4 # throw away real weight + # throw away less weight with a more forgiving truncation + @test maximum(loose.truncation_errors) < maximum(tight.truncation_errors) + # in real time with `normalize = false` the norm loss covers every cut of the step, + # of which the reported last cut per bond is a subset + @test norm(ψ₀)^2 - norm(ψ)^2 >= sum(abs2, tight.truncation_errors) - 1.0e-12 + end + + @testset "history 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 length(total.truncation_errors) == nsteps + @test all(==(L - 1) ∘ length, total.truncation_errors) + @test first(total.truncation_errors) ≈ step.truncation_errors + 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))]