Algorithms

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.

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:

MPSKit.find_groundstate Function

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 (see the automatic-selection notes below).

Arguments

• ψ₀::AbstractMPS: initial guess

• H::AbstractMPO: operator for which to find the ground state

• [environments]: MPS environment manager

• algorithm: optimization algorithm

Keyword Arguments

• tol::Float64 = 1.0e-10: tolerance for the convergence criterion

• maxiter::Int = 200: maximum number of iterations

• verbosity::Int = 3: display progress information

• trunc = nothing: if supplied, a truncation strategy that enables bond-dimension growth through a two-site algorithm (see below)

Automatic algorithm selection

When no algorithm is passed, the choice depends on the type of ψ₀:

• InfiniteMPS: VUMPS (with its tolerance floored at 1e-4), refined by GradientGrassmann when tol < 1e-4. If trunc is given, an IDMRG2 stage is prepended to grow the bond dimension.

• AbstractFiniteMPS: DMRG. If trunc is given, a DMRG2 stage is prepended to grow the bond dimension.

Because single-site DMRG preserves the bond dimension of ψ₀, passing a trunc (or an explicit two-site algorithm) is the usual way to converge from a low-bond-dimension initial guess such as a product state.

Returns

• ψ::AbstractMPS: converged ground state

• environments: environments corresponding to the converged state

• 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, DMRG2 and VUMPS, gradientnorm for GradientGrassmann, and bondresidual for IDMRG and IDMRG2. convergence_measure returns whichever of these is present, for code that only wants the number. A truncating algorithm additionally fills info.max_truncation_error/info.total_truncation_error with what its final sweep discarded. See AlgorithmInfo for the full vocabulary, and The error convention in the manual.

Examples

Ground state of a 4-site transverse-field Ising chain, H = -∑ XₖXₖ₊₁ - ∑ Zₖ, starting from a product state and letting DMRG2 grow the bond dimension:

julia> X = TensorMap(Float64[0 1; 1 0], ℂ^2, ℂ^2);

julia> Z = TensorMap(Float64[1 0; 0 -1], ℂ^2, ℂ^2);

julia> L = 4; lattice = fill(ℂ^2, L);

julia> H = FiniteMPOHamiltonian(lattice, ((i, i + 1) => -(X ⊗ X) for i in 1:(L - 1))) +
           FiniteMPOHamiltonian(lattice, ((i,) => -Z for i in 1:L));

julia> ψ₀ = FiniteMPS(ones(Float64, (ℂ^2)^L));

julia> ψ, envs, info = find_groundstate(ψ₀, H; verbosity = 0, trunc = truncrank(16));

julia> round(real(expectation_value(ψ, H)); digits = 4)
-4.7588

source

The returned error measures convergence to a variational fixed point, which is not the same as accuracy; see Ground state accuracy.

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. Here, we enumerate some of their properties in hopes of pointing you in the right direction. For convenience, the full list of algorithms is:

• DMRG

• DMRG2

• VUMPS

• Gradient descent

• TDVP

• Time evolution MPO

• Quasiparticle Ansatz

• Finite excitations

• "Chepiga Ansatz"

DMRG

Probably the most widely used algorithm for optimizing ground states with MPS is DMRG 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. This can become particularly relevant in the cases where symmetries are involved, since then finding a good distribution of charges is also required. To circumvent this, it is also possible to optimize over two sites at the same time with DMRG2, followed by a truncation back to the single site states. This can dynamically change the bond dimension but comes at an increase in cost.

MPSKit.DMRG Type

struct DMRG{A, F, E, G, B} <: MPSKit.Algorithm

Density Matrix Renormalization Group algorithm for finding the dominant eigenvector.

Each site update is, in order: (1) an optional bond expansion (alg_expand), (2) a single-site eigensolve, and (3) a gauge step (alg_gauge). With the defaults (alg_expand = nothing and alg_gauge = nothing, a non-truncating QR gauge derived from trunc = notrunc()) this is textbook single-site DMRG, which cannot change the bond dimension. Setting alg_expand to a bond-expansion algorithm (e.g. OptimalExpand, RandExpand, SketchedExpand) expands the bond with directions orthogonal to the current state ahead of each eigensolve, recovering Controlled Bond Expansion (CBE) DMRG. Setting alg_gauge to a bond-expanding gauge algorithm (e.g. DMRG3S) instead expands the bond as part of the gauge step, after the eigensolve. Either way, a truncating gauge (see below) is then desirable to cut the enlarged bond back down.

Choosing the gauge

By default, alg_gauge is built for you from trunc/alg_svd/alg_orth: trunc = notrunc() (the default) gives a QR decomposition (alg_orth, Householder by default), any other trunc gives a truncated SVD (alg_svd with that trunc).

DMRG()                            # QR gauge, no truncation
DMRG(; trunc = truncdim(50))   # truncated SVD gauge

To use a bond-expanding gauge such as DMRG3S, pass it directly as alg_gauge; trunc etc. are still routed through to build the inner gauge it wraps, exactly as above:

DMRG(; alg_gauge = DMRG3S(0.1, ExponentialDecay(0.7)), trunc = truncdim(50))

If alg_gauge is instead given with its inner gauge already set (e.g. DMRG3S(0.1, sched, some_gauge)), trunc/alg_svd/alg_orth must be left at their defaults — passing both is an error, since it leaves two conflicting sources for the same setting.

Fields

• tol::Float64: convergence tolerance on the Galerkin error (the tangent-space gradient norm), reported as the galerkin entry of the returned AlgorithmInfo. This acts as a floor: the stopping test is ϵ ≤ max(tol, maximum(ϵ_trunc)), which reduces to ϵ ≤ tol when nothing is truncated

• maxiter::Int64: maximal amount of iterations

• verbosity::Int64: setting for how much information is displayed

• alg_eigsolve::Any: algorithm used for the eigenvalue solvers

• finalize::Any: callback function applied after each iteration, of signature finalize(iter, ψ, H, envs) -> ψ, envs

• alg_expand::Any: algorithm used to expand the bond ahead of each local update, or nothing for none

• alg_gauge::Any: gauge algorithm applied after each local update: NoExpand for a plain gauge step (a QR algorithm with no truncation, or a truncated SVD), or an algorithm that additionally expands the bond beforehand (e.g. DMRG3S)

• backend::Any: backend for tensor contractions and index manipulations

See also

Used as the algorithm argument of find_groundstate and approximate.

source

MPSKit.DMRG2 Type

struct DMRG2{A, G, F, B} <: MPSKit.Algorithm

Two-site DMRG algorithm for finding the dominant eigenvector.

Fields

• tol::Float64: convergence tolerance on the Galerkin error (the tangent-space gradient norm), reported as the galerkin entry of the returned AlgorithmInfo. This acts as a floor: the stopping test is ϵ ≤ max(tol, maximum(ϵ_trunc)), which reduces to ϵ ≤ tol when nothing is truncated

• maxiter::Int64: maximal amount of iterations

• verbosity::Int64: setting for how much information is displayed

• alg_eigsolve::Any: algorithm used for the eigenvalue solvers

• alg_gauge::Any: factorization used for the post-update gauge: a truncated SVD (alg_svd with trunc)

• finalize::Any: callback function applied after each iteration, of signature finalize(iter, ψ, H, envs) -> ψ, envs

• backend::Any: backend for tensor contractions and index manipulations

See also

Used as the algorithm argument of find_groundstate and approximate.

source

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 bond dimension at a higher cost.

MPSKit.IDMRG Type

struct IDMRG{A, B} <: MPSKit.Algorithm

Single site infinite DMRG algorithm for finding the dominant eigenvector.

Fields

• tol::Float64: convergence tolerance, compared against the change in the center bond tensor over a sweep, reported as the bondresidual entry of the returned AlgorithmInfo. This is a fixed-point residual measuring how much a sweep still moves the state, and is not equivalent to the Galerkin error that DMRG and VUMPS report

• maxiter::Int64: maximal amount of iterations

• verbosity::Int64: setting for how much information is displayed

• alg_gauge::Any: algorithm used for gauging the MPS

• alg_eigsolve::Any: algorithm used for the eigenvalue solvers

• backend::Any: backend for tensor contractions and index manipulations

See also

Used as the algorithm argument of find_groundstate, leading_boundary, and approximate.

source

MPSKit.IDMRG2 Type

struct IDMRG2{A, S, B} <: MPSKit.Algorithm

Two-site infinite DMRG algorithm for finding the dominant eigenvector.

Fields

• tol::Float64: convergence tolerance, compared against the change in the center bond tensor over a sweep, reported as the bondresidual entry of the returned AlgorithmInfo. This is a fixed-point residual measuring how much a sweep still moves the state, and is not equivalent to the Galerkin error that DMRG2 reports

• maxiter::Int64: maximal amount of iterations

• verbosity::Int64: setting for how much information is displayed

• alg_gauge::Any: algorithm used for gauging the MPS

• alg_eigsolve::Any: algorithm used for the eigenvalue solvers

• alg_svd::Any: algorithm used for the singular value decomposition

• trunc::MatrixAlgebraKit.TruncationStrategy: algorithm used for truncation of the two-site update

• backend::Any: backend for tensor contractions and index manipulations

See also

Used as the algorithm argument of find_groundstate, leading_boundary, and approximate.

source

VUMPS

VUMPS is an (I)DMRG inspired algorithm that can be used to variationally find the ground state as a Uniform (infinite) Matrix Product State. In particular, a local update is followed by a re-gauging procedure that effectively replaces the entire network with the newly updated tensor. Compared to IDMRG, this often achieves a higher rate of convergence, since updates are felt throughout the system immediately. Nevertheless, this algorithm only works whenever the state is injective, i.e. there is a unique ground state. Since VUMPS is a single-site algorithm, it cannot alter the bond dimension.

MPSKit.VUMPS Type

struct VUMPS{F, B} <: MPSKit.Algorithm

Variational optimization algorithm for uniform matrix product states, based on the combination of DMRG with matrix product state tangent space concepts.

Fields

• tol::Float64: convergence tolerance, compared against the Galerkin error (the tangent-space gradient norm), reported as the galerkin entry of the returned AlgorithmInfo

• maxiter::Int64: maximal amount of iterations

• verbosity::Int64: setting for how much information is displayed

• alg_gauge::Any: algorithm used for gauging the InfiniteMPS

• alg_eigsolve::Any: algorithm used for the eigenvalue solvers

• alg_environments::Any: algorithm used for the MPS environments

• finalize::Any: callback function applied after each iteration, of signature finalize(iter, ψ, H, envs) -> ψ, envs

• backend::Any: backend for tensor contractions and index manipulations

See also

Used as the algorithm argument of find_groundstate and leading_boundary.

References

DocumenterCitations.CitationSiteNode("zauner-stauber2018-cite-1")

• DocumenterCitations.CitationSiteNode("vanderstraeten2019-cite-1")

source

Gradient descent

Both finite and infinite matrix product states can be parametrized by a set of isometric tensors, which we can optimize over. Making use of the geometry of the manifold (a Grassmann manifold), we can greatly outperform naive optimization strategies. Compared to the other algorithms, quite often the convergence rate in the tail of the optimization procedure is higher, such that often the fastest method combines a different algorithm far from convergence with this algorithm close to convergence. Since this is again a single-site algorithm, there is no way to alter the bond dimension.

MPSKit.GradientGrassmann Type

struct GradientGrassmann{O<:OptimKit.OptimizationAlgorithm, F, B} <: MPSKit.Algorithm

Variational gradient-based optimization algorithm that keeps the MPS in left-canonical form, as points on a Grassmann manifold. The optimization is then a Riemannian gradient descent with a preconditioner to induce the metric from the Hilbert space inner product.

Constructors

GradientGrassmann(; kwargs...)

Keyword Arguments

• method = ConjugateGradient: instance of optimization algorithm, or type of optimization algorithm to construct

• finalize!: finalizer algorithm

• tol = Defaults.tol: convergence tolerance, compared against the norm of the Riemannian (Grassmann) gradient reported by the optimizer, which find_groundstate returns as the gradientnorm entry of its AlgorithmInfo.

• maxiter = Defaults.maxiter: maximum amount of iterations

• verbosity = Defaults.verbosity - 1: level of information display

• hasconverged = OptimKit.DefaultHasConverged(tol): convergence criterium

• shouldstop = OptimKit.DefaultShouldStop(maxiter): stopping criterium

• backend = Defaults.backend(): backend for tensor contractions and index manipulations

Fields

• method::OptimKit.OptimizationAlgorithm: optimization algorithm

• finalize!::Any: callback function applied after each iteration, of signature finalize!(x, f, g, numiter) -> x, f, g

• backend::Any: backend for tensor contractions and index manipulations

• hasconverged::Any: function indicating whether the optimization has converged, of signature hasconverged(x, f, g, normgrad) -> Bool

• shouldstop::Any: function indicating whether the optimization should terminate, of signature shouldstop(x, f, g, numfg, numiter, t) -> Bool

See also

Used as the algorithm argument of find_groundstate and leading_boundary.

References

DocumenterCitations.CitationSiteNode("hauru2021-cite-1")

source

Time evolution

Given a particular state, it can also often be useful to examine the evolution of certain properties over time. To that end, there are two main approaches to solving the Schrödinger equation in MPSKit.

MPSKit.timestep Function

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:  .

Arguments

• ψ₀::AbstractMPS: initial state

• H::AbstractMPO: operator that generates the time evolution (can be time-dependent).

• t::Number: starting time of time-step

• dt::Number: time-step magnitude

• alg: algorithm to use for the time evolution, e.g. TDVP or TDVP2.

• envs: MPS environment manager

Keyword Arguments

• imaginary_evolution::Bool = false: if true, the time evolution is done with an imaginary time step instead, (i.e.   instead of  ). This can be useful to compute the ground state of a Hamiltonian, or to compute finite-temperature properties of a system.

• normalize::Bool = false: if true, the state is renormalized after every step, which can be useful to retain numerical stability when the norm loss is not information that is needed.

Returns

• ψ: the time-stepped state

• envs: the updated environment manager

• info::AlgorithmInfo: what the step truncated (see below)

Truncation error

A step performs many local factorisations, each discarding some weight. Rather than collapse those into one number, info reports both aggregations under names that say what they are: info.max_truncation_error is the largest single one (size-independent, comparable against trunc and across runs) and info.total_truncation_error sums them in squares.

Both are non-zero only for algorithms that truncate (TDVP2, BUG with a trunc, and TDVP with a bond expansion). A finite-system step that truncates but happened to discard nothing reports them as exactly 0, whereas infinite one-site TDVP never truncates and reports no truncation entries at all. Neither case means the step was exact, but rather that this particular source of error is either absent or idle.

See AlgorithmInfo for the entries, and Time evolution accuracy in the manual for the other error sources.

Examples

Real-time evolution of the |+···+⟩ product state under a transverse field H = ∑ Zₖ. Each spin precesses independently, so ⟨Xₖ(t)⟩ = cos(2t); after a step dt = 0.1 this is cos(0.2) ≈ 0.980067. The initial state must be complex, since real-time evolution multiplies by -i:

julia> X = TensorMap(ComplexF64[0 1; 1 0], ℂ^2, ℂ^2);

julia> Z = TensorMap(ComplexF64[1 0; 0 -1], ℂ^2, ℂ^2);

julia> ψ₀ = FiniteMPS(ones(ComplexF64, (ℂ^2)^4));

julia> H = FiniteMPOHamiltonian(fill(ℂ^2, 4), ((i,) => Z for i in 1:4));

julia> ψ, envs = timestep(ψ₀, H, 0.0, 0.1, TDVP());

julia> round(real(expectation_value(ψ, 2 => X)); digits = 6)
0.980067

source

MPSKit.time_evolve Function

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.

Arguments

• ψ₀::AbstractMPS: initial state

• H::AbstractMPO: operator that generates the time evolution (can be time-dependent).

• t_span::AbstractVector{<:Number}: time points over which the time evolution is stepped

• alg: algorithm to use for the time evolution, e.g. TDVP or TDVP2.

• envs: MPS environment manager

Keyword Arguments

• verbosity::Int = 0: verbosity level for logging

• imaginary_evolution::Bool = false: if true, the time evolution is done with an imaginary time step instead, (i.e.   instead of  ). This can be useful to compute the ground state of a Hamiltonian, or to compute finite-temperature properties of a system.

• normalize::Bool = false: if true, the state is renormalized after every step, which can be useful to retain numerical stability when the norm loss is not information that is needed.

Returns

• ψ: the time-stepped state

• envs: the updated environment manager

• info::AlgorithmInfo: the truncation performed over the whole evolution, accumulated from the individual steps, together with numiter, the number of steps taken. An algorithm that never truncates reports no truncation entries at all. See AlgorithmInfo and Time evolution accuracy in the manual for the difference and when to use which reported error measure, and timestep for what neither measures.

max_truncation_error is logged per step at verbosity ≥ 3 and for the whole evolution at verbosity ≥ 2. The size-independent measure is used here for the same reason as the ground state algorithms.

source

MPSKit.make_time_mpo Function

make_time_mpo(H::MPOHamiltonian, dt::Number, alg; kwargs...) -> O::MPO

Construct an MPO that approximates  .

Keyword Arguments

• imaginary_evolution::Bool = false: if true, the time evolution is done with an imaginary time step instead, (i.e.   instead of  ). This can be useful to compute the ground state of a Hamiltonian, or to compute finite-temperature properties of a system.

source

TDVP

The first is focused around approximately solving the equation for a small timestep, and repeating this until the desired evolution is achieved. 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 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 evolves two sites at a time and splits the result back apart with a truncated SVD.

• TDVP 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 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.

MPSKit.TDVP Type

struct TDVP{A, E, G, F, B} <: MPSKit.Algorithm

Single site MPS time-evolution algorithm based on the Time-Dependent Variational Principle.

For finite MPS, setting alg_expand to a bond-expansion algorithm (e.g. OptimalExpand, SketchedExpand) expands the bond with directions orthogonal to the current state ahead of each local integration, recovering Controlled Bond Expansion (CBE) TDVP and lifting the fixed-bond limitation of plain single-site TDVP. A truncating trunc is then required to cut the enlarged bond back down (selecting the truncated-SVD gauge). The expansion is 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. In real time this is exact, namely the squared norm drops by precisely the truncated ("discarded") weight,   , with total_truncation_error from the AlgorithmInfo returned by timestep. In imaginary time the norm also carries the physical decay of the weight, so it no longer isolates the truncation. Without trunc nothing is discarded at all and the norm is conserved exactly in real time.

Pass normalize = true to timestep/time_evolve to renormalize at every step instead, like a ground state search. This is independent of imaginary_evolution. CBE is only available for finite MPS.

Fields

• integrator::Any: algorithm used in the exponential solvers

• tolgauge::Float64: tolerance for gauging algorithm

• gaugemaxiter::Int64: maximal amount of iterations for gauging algorithm

• alg_expand::Any: algorithm used to expand the bond ahead of each local update, or nothing for none (finite CBE-TDVP)

• alg_gauge::Any: factorization used for the post-update gauge: a QR algorithm (no truncation) or a truncated SVD

• finalize::Any: callback function applied after each iteration, of signature finalize(t, ψ, H, envs) -> ψ, envs

• backend::Any: backend for tensor contractions and index manipulations

See also

Used as the algorithm argument of timestep, timestep! and time_evolve.

References

DocumenterCitations.CitationSiteNode("haegeman2011-cite-1")

source

MPSKit.TDVP2 Type

struct TDVP2{A, S, F, B} <: MPSKit.Algorithm

Two-site MPS time-evolution algorithm based on the Time-Dependent Variational Principle. See TDVP for more information.

Fields

• integrator::Any: algorithm used in the exponential solvers

• tolgauge::Float64: tolerance for gauging algorithm

• gaugemaxiter::Int64: maximal amount of iterations for gauging algorithm

• alg_svd::Any: algorithm used for the singular value decomposition

• trunc::MatrixAlgebraKit.TruncationStrategy: algorithm used for truncation of the two-site update

• finalize::Any: callback function applied after each iteration, of signature finalize(t, ψ, H, envs) -> ψ, envs

• backend::Any: backend for tensor contractions and index manipulations

See also

Used as the algorithm argument of timestep, timestep! and time_evolve.

References

DocumenterCitations.CitationSiteNode("haegeman2011-cite-2")

source

MPSKit.BUG Type

struct BUG{A, O, G, F, B} <: MPSKit.Algorithm

Single-site time-evolution algorithm for finite MPS, based on the Basis-Update & Galerkin (BUG) integrator, an unconventional robust integrator for dynamical low-rank approximation with an observed second-order convergence.

Unlike TDVP, BUG advances both the basis (K-step) and the core (Galerkin C-step) tensors forward in time, with no backward-in-time substep. This makes it a more natural choice for imaginary-time (dissipative) evolution, where the backward core step of the conventional projector-splitting integrator TDVP can become unstable for large timesteps.

Fields

• integrator::Any: algorithm used in the exponential solvers

• alg_orth::Any: algorithm used to orthonormalize the augmented basis [U₀ │ K₁] after each local update

• alg_gauge::Any: factorization used to gauge and truncate the bond ahead of each local update

• finalize::Any: callback function applied after each iteration, of signature finalize(t, ψ, H, envs) -> ψ, envs

• backend::Any: backend for tensor contractions and index manipulations

Algorithm

Each half-sweep visits every site in turn and, at each site, (i) splits off the bond ahead of it (in the sweep direction) with alg_gauge, truncating it back to trunc, (ii) evolves the connecting tensor over dt/2, and (iii) augments the basis with the new directions discovered by the evolved tensor (old basis first, [U₀ │ K₁], orthonormalized with alg_orth)

Notably, this last step does not include any truncation, and is meant to truncate the previous half-sweep's augmentation. As a result, a truncation scheme truncrank(D) will result in a final MPS of dimension 2D. To restore a maximal dimension of D, apply changebonds with an SvdCut algorithm.

Note

By default the state is not renormalized, as the (loss of) norm accumulates useful information. In real time the squared norm lost is exactly the weight discarded by the bond cuts,   , with total_truncation_error from the AlgorithmInfo returned by timestep. In imaginary time the norm also carries the physical decay of the weight and no longer isolates the truncation. Both reported errors are exactly zero when not truncating. Pass normalize = true to timestep/time_evolve to renormalize after every half-sweep instead.

Warning

The reported errors count only the cuts in step (i). Step (iii) augments the basis without truncating, so within a sweep the bond dimension temporarily exceeds trunc and the norm identity above applies to the state returned at the end of the step.

Tip

Pass a TimerOutputs.TimerOutput as timeroutput to timestep/timestep! to obtain a breakdown of the time spent in the three steps above (cut_bond, AC_integrate, augment) over all local updates.

References

DocumenterCitations.CitationSiteNode("ceruti2022-cite-1")

source

Time evolution MPO

The other approach instead tries to first approximately represent the evolution operator, and only then attempts to apply this operator to the initial state. Typically the first step happens through make_time_mpo, while the second can be achieved through approximate. Here, there are several algorithms available

MPSKit.WI Constant

const WI = TaylorCluster(; N = 1, extension = false, compression = false)

First order Taylor expansion for a time-evolution MPO.

source

MPSKit.WII Type

struct WII <: MPSKit.Algorithm

Generalization of the Euler approximation of the operator exponential for MPOs.

Fields

• tol::Float64: tolerance of the Arnoldi exponentiation used to build each local block

• maxiter::Int64: maximal number of iterations of that exponentiation

See also

Used as the algorithm argument of make_time_mpo.

References

DocumenterCitations.CitationSiteNode("zaletel2015-cite-1")

• DocumenterCitations.CitationSiteNode("paeckel2019-cite-1")

source

MPSKit.TaylorCluster Type

struct TaylorCluster <: MPSKit.Algorithm

Algorithm for constructing the Nth order time evolution MPO using the Taylor cluster expansion.

Fields

• N::Int64: order of the Taylor expansion

• extension::Bool: include higher-order corrections

• compression::Bool: approximate compression of corrections, accurate up to order N

See also

Used as the algorithm argument of make_time_mpo.

References

DocumenterCitations.CitationSiteNode("vandamme2024-cite-1")

source

Time evolution has three distinct error sources, only one of which is reported back to the user. See Time evolution accuracy.

Excitations

It might also be desirable to obtain information beyond the lowest energy state of a given system, and study the dispersion relation. While it is typically not feasible to resolve states in the middle of the energy spectrum, there are several ways to target a few of the lowest-lying energy states. None of these report an error. For what limits their accuracy, see Excitation accuracy.

MPSKit.excitations Function

excitations(
        H, algorithm::QuasiparticleAnsatz, ψ::FiniteQP, [left_environments],
        [right_environments]; num = 1
    ) -> (energies, states)
excitations(
        H, algorithm::QuasiparticleAnsatz, ψ::InfiniteQP, [left_environments],
        [right_environments]; num = 1
    ) -> (energies, states)
excitations(
        H, algorithm::FiniteExcited, ψs::NTuple{<:Any, <:FiniteMPS};
        num = 1, init
    ) -> (energies, states)
excitations(
        H, algorithm::ChepigaAnsatz, ψ::FiniteMPS, [envs];
        num = 1, pos = length(ψ) ÷ 2
    ) -> (energies, states)
excitations(
        H, algorithm::ChepigaAnsatz2, ψ::FiniteMPS, [envs];
        num = 1, pos = length(ψ) ÷ 2
    ) -> (energies, states)

Compute the first excited states and their energy gap above a ground state.

Arguments

• H::AbstractMPO: operator for which to find the excitations

• algorithm: optimization algorithm

• ψ::QP: initial quasiparticle guess

• ψs::NTuple{N, <:FiniteMPS}: N first excited states

• [left_environments]: left ground state environment

• [right_environments]: right ground state environment

Keyword Arguments

• num::Int: number of excited states to compute

• solver: algorithm for the linear solver of the quasiparticle environments

• init: initial excited state guess; defaults to a copy of the first state in ψs

• pos: position of perturbation

source

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 ( DocumenterCitations.CitationSiteNode("haegeman2013-cite-1")

).

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.

# Model parameters
g = 10.0
L = 16
H = transverse_field_ising(FiniteChain(L); g)

# Finding the ground state
ψ₀ = FiniteMPS(L, ℂ^2, ℂ^32)
ψ, = find_groundstate(ψ₀, H; verbosity=0)

# Computing excitations using the Quasiparticle Ansatz
Es, ϕs = excitations(H, QuasiparticleAnsatz(), ψ; num=1)
isapprox(Es[1], 2(g - 1); rtol=1e-2)

true

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 computation in the Heisenberg model illustrates this approach.

# Setting up the model and momentum
momentum = π
H = heisenberg_XXX()

# Ground state computation
ψ₀ = InfiniteMPS(ℂ^3, ℂ^48)
ψ, = find_groundstate(ψ₀, H; verbosity=0)

# Excitation calculations
Es, ϕs = excitations(H, QuasiparticleAnsatz(), momentum, ψ)
isapprox(Es[1], 0.41047925; atol=1e-4)

true

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).

g = 10.0
L = 16
H = transverse_field_ising(Z2Irrep, FiniteChain(L); g)
ψ₀ = FiniteMPS(L, Z2Space(0 => 1, 1 => 1), Z2Space(0 => 16, 1 => 16))
ψ, = find_groundstate(ψ₀, H; verbosity=0)
Es, ϕs = excitations(H, QuasiparticleAnsatz(), ψ; num=1, sector=Z2Irrep(1))
isapprox(Es[1], 2(g - 1); rtol=1e-2) # infinite analytical result

true

MPSKit.QuasiparticleAnsatz Type

struct QuasiparticleAnsatz{A, E} <: MPSKit.Algorithm

Optimization algorithm for quasi-particle excitations on top of MPS groundstates.

Constructors

QuasiparticleAnsatz()
QuasiparticleAnsatz(; kwargs...)
QuasiparticleAnsatz(alg)

Create a QuasiparticleAnsatz algorithm with the given eigensolver, or by passing the keyword arguments to Arnoldi.

Fields

• alg::Any: algorithm used for the eigenvalue solvers

• alg_environments::Any: algorithm used for the quasiparticle environments

See also

Used as the algorithm argument of excitations.

References

DocumenterCitations.CitationSiteNode("haegeman2013-cite-2")

source

Finite excitations

For finite systems we can also do something else - find the ground state of the Hamiltonian +       . This is also supported by calling

# Model parameters
g = 10.0
L = 16
H = transverse_field_ising(FiniteChain(L); g)

# Finding the ground state
ψ₀ = FiniteMPS(L, ℂ^2, ℂ^32)
ψ, = find_groundstate(ψ₀, H; verbosity=0)

Es, ϕs = excitations(H, FiniteExcited(), ψ; num=1)
isapprox(Es[1], 2(g - 1); rtol=1e-2)

false

MPSKit.FiniteExcited Type

struct FiniteExcited{A} <: MPSKit.Algorithm

Variational optimization algorithm for excitations of finite MPS by minimizing the energy of

Fields

• gsalg::Any: optimization algorithm

• weight::Float64: energy penalty for enforcing orthogonality with previous states

See also

Used as the algorithm argument of excitations.

source

"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 ( DocumenterCitations.CitationSiteNode("chepiga2017-cite-1")

).

This is supported via the following syntax:

g = 10.0
L = 16
H = transverse_field_ising(FiniteChain(L); g)
ψ₀ = FiniteMPS(L, ComplexSpace(2), ComplexSpace(32))
ψ, envs, = find_groundstate(ψ₀, H; verbosity=0)
E₀ = real(sum(expectation_value(ψ, H, envs)))
Es, ϕs = excitations(H, ChepigaAnsatz(), ψ, envs; num=1)
isapprox(Es[1] - E₀, 2(g - 1); rtol=1e-2) # infinite analytical result

true

In order to improve the accuracy, a two-site version also exists, which varies two neighbouring sites:

Es, ϕs = excitations(H, ChepigaAnsatz2(), ψ, envs; num=1)
isapprox(Es[1] - E₀, 2(g - 1); rtol=1e-2) # infinite analytical result

true

Errors and accuracy

The algorithms that solve for a state, particularly find_groundstate, leading_boundary, approximate, timestep and time_evolve, return an AlgorithmInfo as their last value, describing how they arrived at their result. excitations and changebonds report nothing. What limits their accuracy is covered below all the same. A single bare number could not do this job, because the quantities involved are genuinely different. A convergence measure and a truncation error answer different questions, are not comparable, and not every algorithm produces both. Conversely, two algorithms that both report "a convergence measure" do not necessarily report the same measure. AlgorithmInfo therefore carries a set of named entries, and each algorithm fills in only the ones it actually computes. Entries are keyed by Symbol, and read either as a property (info.max_truncation_error) or by indexing (info[:max_truncation_error]). The usual dictionary interface is supported, in particular haskey, keys, values and get. Nothing is promised that was never measured: asking for an entry an algorithm did not report is an error naming what it did report, rather than a value that was never computed, and the name of the entry says which quantity it is rather than leaving that to prose. Display the returned object (or call keys(info)) to see what a given algorithm reported. The docstrings of every algorithm report these as well.

MPSKit.AlgorithmInfo Type

struct AlgorithmInfo

Information about how an algorithm arrived at its result, returned as the last value by find_groundstate, leading_boundary, approximate, timestep and time_evolve.

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, DMRG2, VUMPS and VOMPS when solving for a state.

• gradientnorm: the norm of the Riemannian (Grassmann) gradient, as supplied by the optimiser. Reported by GradientGrassmann.

• 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 and IDMRG2.

• localchange: the largest relative change of a local tensor over a sweep. Reported by DMRG and DMRG2 inside approximate.

convergence_measure returns whichever of these is present, for code that only wants "the number that was compared against tol" without caring which one it is.

Truncation

Both truncation entries are built from the same per-factorisation quantity, namely the 2-norm of the singular values a single local factorisation discarded, but aggregate it differently, because no single aggregation answers every question:

• max_truncation_error: the largest of them. It is still a per-factorisation quantity rather than a combination of them, so it does not grow with system size or iteration count, which is what makes it comparable between runs. It is also the entry a trunc setting most directly controls, though how directly depends on the strategy.

• total_truncation_error: all of them combined in quadrature, . This grows with system size and iteration count, so unlike max_truncation_error it is not comparable between runs.

• numtrunc: how many of the recorded errors were non-zero, i.e. how many actually discarded anything.

The two error entries also read under the short aliases ϵ_max and ϵ_total (info.ϵ_max, info[:ϵ_total], haskey(info, :ϵ_max)). They are only ever stored under the descriptive names, so keys and displaying/showing them returns one name per quantity.

Which factorisations get recorded is not the same for every algorithm, and this is worth knowing before comparing numtrunc (or total_truncation_error) between them:

• IDMRG, IDMRG2, TDVP2, BUG and Zipup record every factorisation as it happens, so numtrunc is a count of factorisations. A sweep that visits a bond twice contributes twice.

• DMRG and DMRG2 instead keep one slot per update position, overwritten as the sweep passes, and record those slots once at the end. numtrunc is therefore the number of positions whose most recent cut discarded something . This is never more than length(ψ) for DMRG or length(ψ) - 1 for DMRG2, however many sweeps ran and however many SVDs each performed.

The sweeping choice is deliberate: what the returned state still throws away at a bond is the last cut made there, not the sum of every cut ever made there. It does mean numtrunc counts different things in the two families, so read it as "how many recorded errors were non-zero" rather than as a tally of SVD calls.

See Aggregating truncation errors for how the two relate to a trunc setting, and for the per-strategy caveats.

An algorithm that truncates reports all three even on a run where it happened to discard nothing, so max_truncation_error == 0 means "truncated, but cut nothing away", whereas the entries being absent altogether means the algorithm never truncates. Neither says the result is exact. See the manual on Errors and accuracy for what is not measured here.

source

The rest of this section explains what those quantities are, and - equally important - what they do not measure.

The error convention

In theory: A truncation error measures how much a factorisation changed the tensor it acted on: replacing by its rank-restricted approximation costs

Because a truncated SVD keeps the largest singular values, the discarded part is orthogonal to the kept part, and is exactly the 2-norm of the discarded singular values ( DocumenterCitations.CitationSiteNode("schollwoeck2011-cite-1")

).

In practice: That is precisely what MPSKit computes: every factorisation reports the 2-norm of what it discarded, and nothing more. In particular, the value is absolute, i.e. it is not normalised by .

Where the two are often conflated: In DMRG the same quantity is usually called the "discarded weight" and is read as a probability: when is the bond tensor of a normalised state, the are eigenvalues of the reduced density matrix across the cut which sum to 1 for a normalised state, so is the probability of the discarded subspace. Even though viewing it as a probability is intuitively useful, it is not what makes the value well-defined, and it does not always apply here. In particular, under non-renormalising algorithms the state norm drifts away from 1 precisely as truncation accumulates. The probability intuition holds only when the factorised object is normalised, which is not always the case. The definition    holds in every case, though.

What differs between algorithms is how these per-factorisation values are summed up (aggregated) into the numbers they report, which is the subject of the next subsection.

Warning

A convergence measure and a truncation error are unrelated quantities. find_groundstate, leading_boundary and the iterative approximate algorithms fill in a convergence entry (galerkin, gradientnorm, bondresidual or localchange depending on the algorithm), which has no truncation interpretation at all, while the truncating algorithms fill in max_truncation_error/total_truncation_error, which say nothing about convergence. An algorithm that does both fills both, and they should not be compared with each other.

Aggregating truncation errors

This applies to every truncating algorithm. A sweep of DMRG2, IDMRG2, TDVP2, BUG or Zipup performs many local factorisations, each with its own , and there is no single number that answers every question one might ask of them. Rather than pick one and hope the caller wants that one, AlgorithmInfo carries both aggregations under names that say what they are.

max_truncation_error, the largest single , is a worst case. The point of it is that it is still a per-factorisation quantity: it is one of the , not a combination of them, so it does not grow just because the chain is longer or more sweeps or steps were taken. That is what makes it the one to watch over the course of a run, or to compare between runs at different sizes. An increase in max_truncation_error means individual bonds are being cut harder, whereas an increase in total_truncation_error may only mean there were more bonds to cut. This is the same reason the ground state algorithms report a maximum over local gradient norms rather than a total. It is also why max_truncation_error is the one time_evolve logs; the ground state algorithms log their convergence measure instead, since for them convergence rather than truncation is what the sweep is driving.

It is also the field that a trunc setting most directly controls, though how directly depends on the strategy:

• truncerror bounds the discarded weight of each factorisation, which is exactly , so max_truncation_error should come out at or below the tolerance you set.

• trunctol bounds each individual singular value instead. Discarding of them leaves   , so max_truncation_error lands near the tolerance but is not bounded by it.

• truncrank fixes the rank and says nothing about magnitudes at all. Here, max_truncation_error is not something you set but something you read off. It is thus the consequence of that choice of bond dimension.

total_truncation_error, all of them combined in quadrature,

is the one that adds up to something. Note that it is the squares that are summed, and the root taken at the end, because is the quantity that accumulates exactly: each local truncation is an orthogonal projection, so it removes exactly from the squared norm of the tensor it acts on. Summing the squares therefore tracks a running "cost" rather than merely a worst case. The price is that it is extensive: it grows like in the number of truncations, so it is not comparable between different system sizes or sweep counts.

Whether that running cost is also the error of the final state depends on what the algorithm does between truncations, so it is not a property of the aggregation itself. It does hold for real-time evolution, which is worked out in The norm as a record of truncation. A variational sweep, by contrast, renormalises as it goes, so there total_truncation_error is a diagnostic of how hard the truncation is working rather than a norm deficit.

Neither field is the distance to the untruncated solution. Bounding that gives the linear sum   , which is a third quantity again, and always the largest of the three.

Which factorisations are recorded

Both aggregations and the numtrunc count are built from the set of per-factorisation errors the algorithm records. This set is not the same for every algorithm.

IDMRG, IDMRG2, TDVP2, BUG and Zipup record every factorisation as it happens. Thus, numtrunc is a count of factorisations, and a sweep that visits a bond twice contributes two entries.

DMRG and DMRG2 keep one slot per update position instead, overwritten as the sweep passes over it, and record the slots once at the end. In these cases, numtrunc is the number of positions whose most recent cut discarded something, and it can never exceed length(psi) (single-site) or length(psi) - 1 (two-site), no matter how many sweeps ran. A reported numtrunc less than the number of bonds means some bonds last discarded a non-zero weight. These are typically the bonds nearest the two ends.

The difference here is made deliberately. What a returned variational state still throws away at a given bond is the last cut made there, not the sum of every cut ever made there. This means that recording the latest value per position is the more meaningful quantity for a sweeping ground-state algorithm. Thus, numtrunc counts different things in the two families, and total_truncation_error correspondingly sums a different number of terms, so neither is directly comparable between, say, DMRG2 and IDMRG2. Read numtrunc more as "how many recorded errors were non-zero", and less as a tally of SVD calls.

The two error entries also read under the short aliases ϵ_max and ϵ_total. They are only ever stored under the descriptive names, so keys and displaying/showing them returns one name per quantity.

Ground state accuracy

find_groundstate, leading_boundary and the iterative approximate algorithms report the quantity their tol is compared against, together with a converged flag. Because these are not the same quantity from one algorithm to the next, each is stored under a key that names it (galerkin, gradientnorm, bondresidual or localchange). convergence_measure returns whichever of them is present, for code that only wants the number. In theory: Convergence is measured by the (norm of the) variational gradient: the component of that points away from the current state but still lies in the tangent space of the variational manifold. It vanishes exactly at a variational fixed point, and its norm is what most of these algorithms report. For the sweeping algorithms that norm is known as the Galerkin error.

In practice: The sweeping algorithms (DMRG, DMRG2, VUMPS, VOMPS) report the Galerkin error, computed per site as the norm of the local update projected onto the orthogonal complement of the current tensor, and then aggregated as the maximum over sites. GradientGrassmann instead reports the gradient norm supplied by its optimiser, taken over the whole state at once in the preconditioned Grassmann metric.

IDMRG and IDMRG2 report neither, and this is easy to miss because their tol sits alongside the others. Their bondresidual is a fixed-point residual: the change in the center bond tensor from one sweep to the next,  , with both tensors projected onto their common space when the bond space changed. It measures how much a sweep still moves the state, which is weaker than measuring how far the state is from stationarity, since an algorithm crawling through a slow region reports a small value for the same reason a converged one does. For IDMRG2 the projection onto the common space also means that a change in bond dimension between sweeps is projected out of the measure rather than counted in it.

Why they differ: The Galerkin and Grassmann measures are the same underlying gradient, not different physical quantities. They differ in the metric it is measured in (the Grassmann gradient is preconditioned) and in how the per-site contributions are combined (a maximum versus a single global norm). The maximum is a deliberate practical choice, as it keeps the reported number independent of system size. The IDMRG residual is not that gradient at all: it certifies that the sweeps have stopped moving, not that there is nowhere left to move to. The consequence of these mismatches is that a tol tuned for one algorithm is not a tol tuned for another, and this is worth keeping in mind when swapping algorithms at a fixed tol.

In other words, convergence is only defined relative to the manifold you are optimising over. A single-site algorithm at a fixed bond dimension can drive its convergence measure to machine precision and still be far from the true ground state, because the error that remains is the bond dimension itself, which no amount of further sweeping can address. A small convergence measure certifies a fixed point, not an accurate state. Growing the bond dimension is the job of the two-site algorithms (DMRG2, IDMRG2) or of a bond expansion (DMRG with an alg_expand, or an expanding alg_gauge such as DMRG3S); see also changebonds.

Once an algorithm does truncate, the two error notions interact. The Galerkin error cannot fall below the level set by the weight being discarded each sweep, so a truncating scheme converges once galerkin reaches the truncation error rather than the (unreachable) bare tol. DMRG/DMRG2 account for this: their stopping test is ϵ ≤ max(tol, maximum(ϵ_trunc)), which reduces to the plain ϵ ≤ tol when nothing is truncated.

Neither measure is an error bar on an observable, and no cheap substitute for one exists. The energy variance   is an independent and more demanding measure. Note what it actually quantifies, namely how far the state is from being an exact eigenstate, which is not the same thing as the error on some other observable. The usual practical route is to compute a series of states at increasing bond dimension and extrapolate observables against the variance towards zero ( DocumenterCitations.CitationSiteNode("hubig2018-cite-1")

). This is an empirical extrapolation rather than a bound.

Time evolution accuracy

Unlike a ground state search, a time evolution has no convergence criterion to run to. There is no fixed point, and the error is made at every step. Three sources behave differently and only one of them is reported.

• Truncation error. Whenever a bond is cut back down, the discarded singular values are lost from the state. timestep and time_evolve report this through the max_truncation_error/total_truncation_error entries of the AlgorithmInfo they return. It is the error you control through the algorithm's trunc, and the only one that is free to compute, since the truncating SVD produces it anyway. It is non-zero for TDVP2, for BUG with a trunc, and for TDVP with a bond expansion. Plain single-site TDVP runs at fixed bond dimension and reports exactly 0.

• Projection error. Single-site TDVP confines the evolution to the tangent space of a fixed-bond-dimension manifold,  . The component of the exact evolution pointing off that manifold is simply dropped, and this happens even with no truncation and exact local solves. It is not reported, since measuring it costs an extra effective-Hamiltonian application per site. This is what a bond expansion (CBE) exists to reduce ( DocumenterCitations.CitationSiteNode("li2024-cite-1")

  ).

• Splitting error. The evolution is not applied in one piece. Rather, it is split into local terms that are integrated in sequence, which do not commute. This is a Trotter-type error. Note, however, that what is split differs per method: TDVP splits the tangent-space projector into site and bond terms, while the MPO methods (WI, WII, TaylorCluster) split the Hamiltonian in the more familiar sense. For TDVP's symmetric back-and-forth sweep the result is globally ( DocumenterCitations.CitationSiteNode("lubich2015-cite-1")

  , DocumenterCitations.CitationSiteNode("paeckel2019-cite-2")

  ), so it is controlled by dt alone. This can only be estimated by comparing one step of dt against two of dt / 2.

A trustworthy run needs all three under control, not just a small truncation error. In practice: pick dt from a convergence check, pick trunc from the reported truncation error, and use a bond-adaptive scheme (TDVP2, BUG, or TDVP with alg_expand) whenever entanglement grows during the evolution, since a fixed bond dimension silently converts entanglement growth into projection error.

Note that these do not all shrink together, so there is a sweet spot in dt rather than a "smaller is always better". Decreasing dt reduces the splitting error, but it also means more steps to reach the same final time, and every step truncates again. The accumulated truncation error grows with the number of steps, as does the compute time. The projection error does not improve at all, since it represents a rate at which the exact solution leaves the manifold. Shrinking dt in this case only samples this rate more finely. The practical consequence is that below some dt the total error stops improving and eventually gets worse, and the remedy at that point is a larger bond dimension rather than a smaller step.

The norm as a record of truncation

How the per-factorisation errors are aggregated into max_truncation_error/total_truncation_error is described under Aggregating truncation errors; what follows is specific to time evolution.

By default none of the time evolution algorithms renormalize (normalize = false), which is deliberate. In real time the local exponentials are unitary, so truncation is the only thing that changes the norm and it becomes a running record of what truncation has cost,

The reported total_truncation_error is the norm deficit, and this composes across steps. This follows from the following two facts put together, one per half of a local update. 2. An SVD truncation is an orthogonal projection onto the kept Schmidt vectors and is 2-norm optimal at that rank ( DocumenterCitations.CitationSiteNode("schollwoeck2011-cite-2")

), so the kept and discarded parts are orthogonal.

By Pythagoras the squared norm drops by exactly the discarded weight, the usual way of quantifying truncation during a time evolution ( DocumenterCitations.CitationSiteNode("paeckel2019-cite-3")

). 2. The local exponentials of the projector-splitting sweep are unitary, so TDVP conserves the norm and the energy exactly when the local equations are solved exactly ( DocumenterCitations.CitationSiteNode("paeckel2019-cite-4")

), contributing nothing to the norm change.

These two hold for the local updates of a step, so composing over all steps gives the identity.

Note

"Exactly" in the second fact is up to the tolerance of the local exponentials, and is thus in practice only approximate due to integrator tolerance.

It is also specific to real time with normalize = false:

• Imaginary time evolves with the non-unitary  , which rescales the state on its own. The norm then moves for two independent reasons, namely the physical decay of the weight and the truncation loss. One cannot separate them from each other. The reported errors still count only the truncation.

• normalize = true renormalizes at every local update, destroying the identity by construction. This is usually what you want for imaginary-time evolution used as a ground state or thermal-state search. The reported errors are unaffected.

• No truncation at all (plain single-site TDVP, or BUG with a QR gauge) gives  , and in real time the norm is then conserved exactly.

• An InfiniteMPS is regauged to norm 1 per site structurally, so its norm carries no such information and normalize has no effect.

Excitation accuracy

excitations returns only (energies, states): there is no error term, and none of the sources below is reported back to you. They are worth knowing about, because the dominant one is usually not the one the algorithm is working on.

• Inherited ground state error. Every method here builds on a ground state you supply and treats it as exact. Its error propagates straight into the gap, and since a gap is a difference of two large energies, it is typically the limiting factor. A well-converged ground state (in the sense of its convergence measure and bond dimension) is necessary for a meaningful excitation calculation.

• Ansatz limitation. QuasiparticleAnsatz varies over the single-quasiparticle tangent space on top of a fixed ground state. It is variational within that space and well suited to isolated quasiparticle branches, but multi-particle continua are not representable in it, so results there are not to be trusted. How well it does on an isolated branch is controlled by the gaps around the targeted level. The ansatz approximates an eigenvalue that is separated from the rest of the spectrum in its momentum sector with an error bounded exponentially in the size of the local operator's support, at a rate set by the gap below and above that eigenvalue ( DocumenterCitations.CitationSiteNode("haegeman2013-cite-3")

  ). A branch that is nearly degenerate with the ground state, or that sits just below a continuum, therefore converges much more slowly than an isolated one.

• Eigensolver convergence. The eigenvalue problem is solved with KrylovKit, and a run that fails to converge num states emits a warning carrying the residual when the verbosity is set high enough. That residual is neither returned nor thrown, so it is worth not running with warnings suppressed. Its convergence is governed by the same spectral structure as the ansatz above: levels that are well separated converge quickly, while nearly degenerate ones converge slowly and are the ones most likely to come back unconverged.

• Penalty-based orthogonality (FiniteExcited). Higher states are found by minimising   against the previously converged states, with the weight field. A finite weight enforces orthogonality only approximately, so a residual overlap with a lower state biases the reported energy downwards. Since the reported value is the expectation value of the bare H, this bias is invisible in the output. Raising weight suppresses it at the cost of stretching the spectrum and slowing down the eigensolver's per-gap eigensolves.

• Truncation (ChepigaAnsatz2). The two-site excited state is split back to single-site tensors with a truncated SVD governed by trunc, and the resulting discarded weight is not reported.

changebonds

Many of the previously mentioned algorithms do not possess a way to dynamically change to bond dimension. This is often a problem, as the optimal bond dimension is often not a priori known, or needs to increase because of entanglement growth throughout the course of a simulation. changebonds exposes a way to change the bond dimension of a given state.

MPSKit.changebonds Function

changebonds(ψ::AbstractMPS, H, alg, envs) -> ψ′, envs′
changebonds(ψ::AbstractMPS, alg) -> ψ′

Change the bond dimension of ψ using the algorithm alg, and return the new ψ and the new envs. For AbstractInfiniteMPS, changebonds returns new environments without modifying the one provided. changedbonds! can modify both the provided state and environments, depending on the algorithm. For FiniteMPS, changebonds also modifies the environments.

See also: SvdCut, RandExpand, VUMPSSvdCut, OptimalExpand

Examples

Growing the bond dimension of a product state with OptimalExpand, which expands each bond with directions orthogonal to the current state (using the environments of H):

julia> Z = TensorMap(Float64[1 0; 0 -1], ℂ^2, ℂ^2);

julia> ψ = FiniteMPS(ones(Float64, (ℂ^2)^4));

julia> H = FiniteMPOHamiltonian(fill(ℂ^2, 4), ((i, i + 1) => Z ⊗ Z for i in 1:3));

julia> dim(left_virtualspace(ψ, 3))
1

julia> ψ′, envs = changebonds(ψ, H, OptimalExpand(; trunc = truncrank(4)));

julia> dim(left_virtualspace(ψ′, 3))
2

Note

A bond is only expanded if there is something to expand it with. If the projection of the two-site update onto the orthogonal complement of the current state vanishes — for instance when the state is already an exact eigenstate of the local terms, or when the operator does not couple into a symmetry sector yet — that bond is left untouched. Replacing Z ⊗ Z by X ⊗ X above illustrates this: ones(Float64, (ℂ^2)^4) is an eigenstate of every X ⊗ X term, so every bond stays at dimension 1.

source

All of these are controlled by a trunc, and the weight they discard is measured the same way as described under The error convention. 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: 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 algorithm in combination with approximate. Since the output of this method might have a 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.

• OptimalExpand: This algorithm is based on the idea of expanding the bond dimension by investigating the two-site derivative, and adding the most important blocks which are orthogonal to the current state. From the point of view of a local two-site update, this procedure is optimal, but it requires to evaluate a two-site derivative, which can be costly when the physical space is large. The state will remain unchanged, but a one-site scheme will now be able to push the optimization further. The subspace used for expansion can be truncated through trunc, which dictates how many singular values will be added.

• RandExpand: This algorithm similarly adds blocks orthogonal to the current state, but does not attempt to select the most important ones, and rather just selects them at random. The advantage here is that this is much cheaper than the optimal expand, and if the bond dimension is grown slow enough, this still obtains a very good expansion scheme. Again, The state will remain unchanged and a one-site scheme will now be able to push the optimization further. The subspace used for expansion can be truncated through trunc, which dictates how many orthogonal vectors will be added.

• VUMPSSvdCut: This algorithm is based on the VUMPS algorithm, and consists of performing a two-site update, and then truncating the state back down. Because of the two-site update, this can again become expensive, but the algorithm has the option of both expanding as well as truncating the bond dimension. Here, trunc controls the truncation of the full state after the two-site update.

Leading boundary

For statistical mechanics partition functions we want to find the approximate leading boundary MPS. Again this can be done with VUMPS:

th = nonsym_ising_mpo()
ts = InfiniteMPS([ℂ^2],[ℂ^20]);
(ts,envs,_) = leading_boundary(ts,th,VUMPS(maxiter=400,verbosity=false));

If the mpo satisfies certain properties (positive and hermitian), it may also be possible to use GradientGrassmann.

MPSKit.leading_boundary Function

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.

Arguments

• ψ₀::AbstractMPS: initial guess

• O::AbstractMPO: operator for which to find the leading_boundary

• [environments]: MPS environment manager

• algorithm: optimization algorithm

Keyword Arguments

• 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

Returns

• ψ::AbstractMPS: converged leading boundary MPS

• environments: environments corresponding to the converged boundary

• 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 and VOMPS report galerkin, the maximum over sites of the local update projected onto the orthogonal complement of the current tensor.

  • GradientGrassmann reports gradientnorm from its optimiser.

  • IDMRG and IDMRG2 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 returns whichever is present, for code that only wants the number. See AlgorithmInfo, find_groundstate, and the manual on the ϵ convention under The error convention and Ground state accuracy.

source

approximate

Often, it is useful to approximate a given MPS by another, typically by one of a different bond dimension. This is achieved by approximating an application of an MPO to the initial state, by a new state.

MPSKit.approximate Function

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, ψ₀ is approximated directly to ψ (i.e. O is taken to be the identity).

Not every algorithm supports every combination of arguments below — see the per-algorithm notes at the end of this docstring before picking one.

Arguments

• ψ₀::AbstractMPS: initial guess of the approximated state

• (O::AbstractMPO, ψ::AbstractMPS): operator O and state ψ to be approximated

• ψ::AbstractMPS: state to be approximated directly (without an operator)

• algorithm: approximation algorithm. See below for a list of available algorithms.

• [environments]: MPS environment manager

Keyword Arguments

The keyword-based call (no explicit algorithm) is a convenience method that picks an algorithm for you based on the type of ψ₀ (DMRG/DMRG2 for a finite MPS, VOMPS/IDMRG/ IDMRG2 for an infinite MPS) and only accepts the (O, ψ) tuple form of toapprox. Once you pass an explicit algorithm, keywords are no longer accepted here — configure the algorithm struct itself instead (e.g. DMRG(; tol, maxiter, verbosity)).

• tol::Float64: convergence tolerance, compared against the 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.

    convergence_measure returns whichever of these is present, for code that only wants the number.

  • the two-site ones (DMRG2, IDMRG2) which involve a truncated SVD fill the truncation fields with what their final sweep discarded, i.e. what the returned state is still throwing away per sweep rather than what the early, unconverged sweeps did.

  • Zipup is a single non-iterative sweep, so it has no convergence measure at all: it reports no converged entry and none of the convergence entries, and fills the truncation entries instead.

  See AlgorithmInfo for the full list and The error convention 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 (ψ, info); its ψ₀ is a write destination, not an initial guess, and it may be omitted.

Algorithm	Scheme	State ψ₀	bare ψ allowed?	approximate!
DMRG	single-site, fixes bond dim	AbstractFiniteMPS	✅	✅
DMRG2	two-site, truncates via trunc	AbstractFiniteMPS	✅	✅
Zipup	streaming MPO-MPS compression	FiniteMPS destination, optional	❌ (tuple only)	✅
IDMRG	single-site, thermodynamic limit	InfiniteMPS / MultilineMPS	❌ (tuple only)	✅
IDMRG2	two-site, thermodynamic limit, needs unit cell ≥ 2	InfiniteMPS / MultilineMPS	❌ (tuple only)	✅
VOMPS	tangent-space truncation	InfiniteMPS / MultilineMPS	❌ (tuple only)	❌ (out-of-place only)

InfiniteMPS/InfiniteMPO inputs are converted internally to MultilineMPS/MultilineMPO for IDMRG, IDMRG2, and VOMPS; you can also pass those types directly.

source

Varia

What follows is a medley of lesser known (or used) algorithms and don't entirely fit under one of the above categories.

Dynamical DMRG

Dynamical DMRG has been described in other papers and is a way to find the propagator. The basic idea is that to calculate    , one can variationally find    and then the propagator simply equals  .

MPSKit.propagator Function

propagator(ψ₀::AbstractFiniteMPS, z::Number, H::MPOHamiltonian, alg::DynamicalDMRG; init = copy(ψ₀)) -> (g, ψ)

Calculate the action of the propagator   using the dynamical DMRG algorithm.

Returns

• g: approximation of the propagator matrix element   

• ψ: MPS approximation of  

source

MPSKit.DynamicalDMRG Type

struct DynamicalDMRG{F<:MPSKit.DDMRG_Flavour, S, B} <: MPSKit.Algorithm

A dynamical DMRG method for calculating dynamical properties and excited states, based on a variational principle for dynamical correlation functions.

Fields

• flavour::MPSKit.DDMRG_Flavour: flavour of the algorithm to use, either of type NaiveInvert or Jeckelmann

• solver::Any: algorithm used for the linear solvers

• tol::Float64: convergence tolerance, compared against the largest change in a center tensor over a sweep, maxᵢ ‖ACᵢ′ - ACᵢ‖. This represents a measure of how much the sweep still moves the state, not a residual of the linear system (that is controlled by solver)

• maxiter::Int64: maximal amount of iterations

• verbosity::Int64: setting for how much information is displayed

• backend::Any: backend for tensor contractions and index manipulations

See also

Used as the algorithm argument of propagator.

References

DocumenterCitations.CitationSiteNode("jeckelmann2002-cite-1")

source

MPSKit.NaiveInvert Type

struct NaiveInvert <: MPSKit.DDMRG_Flavour

An alternative approach to the dynamical DMRG algorithm, without quadratic terms but with a less controlled approximation. This algorithm minimizes the following cost function

Returns the approximation of    and  .

See also

Jeckelmann for the original approach.

source

MPSKit.Jeckelmann Type

struct Jeckelmann <: MPSKit.DDMRG_Flavour

The original flavour of dynamical DMRG, which minimizes functional (14) from Jeckelmann2002. Writing   and  , this is

which attains its minimum at

Together with equation (11) from that same paper we can determine the full propagator  .

Returns the approximation of    and  .

See also

NaiveInvert for a less costly but less accurate alternative.

References

DocumenterCitations.CitationSiteNode("jeckelmann2002-cite-2")

source

fidelity susceptibility

The fidelity susceptibility measures how much the ground state changes when tuning a parameter in your Hamiltonian. Divergences occur at phase transitions, making it a valuable measure when no order parameter is known.

MPSKit.fidelity_susceptibility Function

fidelity_susceptibility(
    state::Union{FiniteMPS, InfiniteMPS}, H₀::T,
    Vs::AbstractVector{T}, [henvs = environments(state, H₀, state)];
    maxiter = Defaults.maxiter,
    tol = Defaults.tol
) where {T <: MPOHamiltonian}

Computes the fidelity susceptibility of a the ground state state of a base Hamiltonian H₀ with respect to a set of perturbing Hamiltonians Vs. Each of the perturbing Hamiltonians can be interpreted as corresponding to a tuning parameter in a 'total' Hamiltonian    .

Returns a matrix containing the overlaps of the elementary excitations on top of state corresponding to each of the perturbing Hamiltonians.

source

Boundary conditions

You can impose periodic or open boundary conditions on an infinite Hamiltonian, to generate a finite counterpart. In particular, for periodic boundary conditions we still return an MPO that does not form a closed loop, such that it can be used with regular matrix product states. This is straightforward to implement but, and while this effectively squares the bond dimension, it is still competitive with more advanced periodic MPS algorithms.

MPSKit.open_boundary_conditions Function

open_boundary_conditions(mpo::InfiniteMPO, L::Int) -> FiniteMPO

Convert an infinite MPO into a finite MPO of length L, by applying open boundary conditions.

source

open_boundary_conditions(mpo::InfiniteMPOHamiltonian, L::Int) -> FiniteMPOHamiltonian

Convert an infinite MPO into a finite MPO of length L, by applying open boundary conditions.

source

MPSKit.periodic_boundary_conditions Function

periodic_boundary_conditions(mpo::AbstractInfiniteMPO, L::Int)

Convert an infinite MPO into a finite MPO of length L, by mapping periodic boundary conditions onto an open system.

source

Exact diagonalization

As a side effect, our code supports exact diagonalization. The idea is to construct a finite matrix product state with maximal bond dimension, and then optimize the middle site. Because we never truncate the bond dimension, this single site effectively parametrizes the entire Hilbert space.

MPSKit.exact_diagonalization Function

exact_diagonalization(
        H::FiniteMPOHamiltonian;
        sector = rightunit(H), num::Int = 1, which::Symbol = :SR,
        alg = Defaults.alg_eigsolve(; dynamic_tols = false),
        backend = Defaults.backend()
    ) -> vals, state_vecs, convhist

Use KrylovKit.eigsolve to perform exact diagonalization on a FiniteMPOHamiltonian to find its eigenvectors as FiniteMPS of maximal rank, essentially equivalent to dense eigenvectors.

Arguments

• H::FiniteMPOHamiltonian: the Hamiltonian to diagonalize.

Keyword Arguments

• sector = rightunit(H): the total charge of the eigenvectors, which is chosen trivial by default.

• num::Int = 1: the number of eigenvectors to find.

• which::Symbol = :SR: the kind eigenvalues to find, see KrylovKit.eigsolve.

• alg = Defaults.alg_eigsolve(; dynamic_tols = false): the diagonalization algorithm to use, see KrylovKit.eigsolve.

• backend = Defaults.backend(): backend for tensor contractions and index manipulations.

Valid sector values

The total charge of the eigenvectors is imposed by adding a charged auxiliary space as the leftmost virtualspace of each eigenvector. Specifically, this is achieved by passing left = Vect[typeof(sector)](sector => 1) to the FiniteMPS constructor. As such, the only valid sector values (i.e. sector values for which the corresponding eigenstates have valid fusion channels) are those that occur in the dual of the fusion of all the physical spaces in the system.

source