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.

Groundstates

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, ϵ)
find_groundstate(ψ₀, H, algorithm, [environments]) -> (ψ, environments, ϵ)

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

• ϵ::Float64: final convergence error upon terminating the algorithm, i.e. the quantity compared against the algorithm's tol. It measures distance from a variational fixed point and is not a truncation error. See the manual under The error convention for more information. Which measure it is depends on the algorithm: the sweeping algorithms (DMRG, DMRG2, VUMPS, IDMRG, IDMRG2) report the Galerkin error, whereas GradientGrassmann reports the norm of the Riemannian gradient from its optimizer.

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, ϵ = 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 groundstates 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, compared against the Galerkin error (the tangent-space gradient norm)

• 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, compared against the Galerkin error (the tangent-space gradient norm)

• 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 Galerkin error (the tangent-space gradient norm)

• 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 Galerkin error (the tangent-space gradient norm)

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

• 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

• Zauner-Stauber et al. Phys. Rev. B 97 (2018)

• Vanderstraeten et al. SciPost Phys. Lect. Notes 7 (2019)

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. This is also the ϵ returned by find_groundstate for this algorithm.

• 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

• Hauru et al. SciPost Phys. 10 (2021)

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, ϵ)
timestep!(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ₀, envs, ϵ)

Time-step the state ψ₀ with Hamiltonian H over a given time step dt at time t, solving the Schroedinger equation:  .

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

• ϵ: the truncation error of the step (see below)

Truncation error

ϵ is the truncation error of the step, i.e. the 2-norm of the sum of squared singular values discarded by the local gauge factorizations. In other words, ϵ² is the truncated ("discarded") weight.

It is nonzero only for algorithms that truncate such as TDVP2, BUG with a trunc, and TDVP with a bond expansion, while it is exactly 0 for one-site TDVP, which runs at fixed bond dimension. A zero ϵ does not mean the step was exact, but that this particular error channel is absent. In particular, the projection and time-discretisation errors are never included in ϵ.

In real time with normalize = false, ϵ is exactly the norm lost to truncation,   .

See Time evolution accuracy in the manual for what the other error sources are, why the per-bond errors combine in a squared manner, and the precise statement and caveats of the norm identity.

Examples

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, ϵ)
time_evolve!(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ₀, envs, ϵ)

Time-evolve the initial state ψ₀ with Hamiltonian H over a given time span by stepping through each of the time points obtained by iterating t_span.

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

• ϵ: the truncation error accumulated over the whole evolution, i.e. the per-step errors of timestep,  . In real time with normalize = false this is exactly the norm lost to truncation over the whole evolution,   . This is zero for algorithms that never truncate. See timestep for what it does and does not measure.

The per-step error is logged at verbosity ≥ 3 and the accumulated total at verbosity ≥ 2.

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 ϵ the truncation error 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

• Haegeman et al. Phys. Rev. Lett. 107 (2011)

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

• Haegeman et al. Phys. Rev. Lett. 107 (2011)

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 ϵ the truncation error returned by timestep. In imaginary time the norm also carries the physical decay of the weight and no longer isolates the truncation. ϵ is exactly zero when not truncating. Pass normalize = true to timestep/time_evolve to renormalize after every half-sweep instead.

Warning

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

Tip

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

References

• Ceruti et al. BIT Numer. Math. 62 (2022)

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

• Zaletel et al. Phys. Rev. B 91 (2015)

• Paeckel et al. Ann. of Phys. 411 (2019)

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

• Van Damme et al. SciPost Phys. 17 (2024)

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 (Haegeman et al., 2013).

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

• Haegeman et al. Phys. Rev. Let. 111 (2013)

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 (Chepiga and Mila, 2017).

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

Most algorithms in MPSKit report an error alongside their result, and the manual pages above refer to it as ϵ throughout. That single name covers genuinely different quantities, and the differences matter. This section clarifies the differences, and explains what they do and, just as important, what they don't measure.

The error convention

Where ϵ is a truncation error, it is the 2-norm of the discarded singular values of a single factorisation, so that ϵ² is the discarded weight and the squared norm of the factorised tensor drops by exactly ϵ². The name "discarded weight" refers to it truly representing a probability: the are the eigenvalues of the reduced density matrix across the cut, i.e. the statistical weights of the Schmidt states, summing to 1 for a normalised state. So ϵ² is the probability weight thrown away and ϵ is the corresponding amplitude, which is why squares appear wherever these errors are combined. Note that ϵ is absolute rather than relative: it is not divided by the norm of the state, which under a non-renormalising algorithm drifts away from 1 precisely as truncation accumulates.

What differs between algorithms is how the per-bond values are aggregated. For this reason, values between algorithms are not directly comparable.

Warning

Not every ϵ is a truncation error. The ϵ returned by find_groundstate, leading_boundary and the iterative approximate algorithms is a convergence measure, with no truncation interpretation at all. The two are unrelated quantities that happen to share a name; see the two sections below.

Ground-state accuracy

find_groundstate, leading_boundary and the iterative approximate algorithms return the quantity their tol is compared against. For the sweeping algorithms (DMRG, DMRG2, VUMPS, IDMRG, IDMRG2) this is the Galerkin error: the norm of the local gradient projected orthogonally to the current state. It vanishes exactly at a variational fixed point. GradientGrassmann instead reports the norm of the Riemannian gradient from its optimizer. Both vanish at a fixed point and both are gradient norms. However, since they are taken in different metrics, their magnitudes are not directly comparable. In particular, a tol tuned for one is not a tol tuned for the other.

In other words, convergence is only defined relative to the manifold you are optimising over. A single-site algorithm at a fixed bond dimension can drive its ϵ to machine precision and still be far from the true ground state, because the error that remains is the bond dimension itself, which no amount of further sweeping can address. A small ϵ certifies a fixed point, not an accurate state. Growing the bond dimension is the job of the two-site algorithms (DMRG2, 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 ϵ 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. For that, the standard route is the energy variance   and extrapolation of observables towards zero variance (Hubig et al.).

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 return this as their third value ϵ, the norm of the discarded component, so that ϵ² is the discarded weight. It is the error you control through the algorithm's trunc, and the only one that is free to compute, since the truncating SVD produces it anyway. It is non-zero for TDVP2, for BUG with a trunc, and for TDVP with a bond expansion. Plain single-site TDVP runs at fixed bond dimension and returns 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 as measuring it costs an extra effective-Hamiltonian application per site. This is what a bond expansion (CBE) exists to reduce (Li et al.).

• Time-discretization error. The projector splitting is globally for the symmetric back-and-forth sweep (Lubich et al., Paeckel et al.), so it is controlled by dt alone. This can only be estimated by comparing one step of dt against two of dt / 2.

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

Summing local errors in squares and relation to the norm

is the quantity that adds exactly. Each local truncation is an orthogonal projection, so it removes exactly from the squared norm, and the substeps between truncations preserve the norm. Summing the squares therefore tracks a conserved "cost". The alternative is to sum :   , a distance to the untruncated solution, which is a different and always larger quantity.

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

The reported ϵ is the norm deficit, and this composes across steps. This follows from the following two facts put together, one per half of a local update. 2. An SVD truncation is an orthogonal projection onto the kept Schmidt vectors and is 2-norm optimal at that rank (Schollwöck), so the kept and discarded parts are orthogonal.

By Pythagoras the squared norm drops by exactly the discarded weight, the usual way of quantifying truncation during a time evolution (Paeckel et al.). 2. The local exponentials of the projector-splitting sweep are unitary, so TDVP conserves the norm and the energy exactly when the local equations are solved exactly (Paeckel et al.), 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. ϵ still counts only the truncation.

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

• No truncation at all (plain single-site TDVP, 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 ϵ 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. For infinite systems the momentum superposition itself is exact, so momentum is a good quantum number and no error enters through it.

• Eigensolver convergence. The local eigenvalue problem is solved with KrylovKit, and a run that fails to converge num states emits a warning on the residual when the verbosity is set high enough. This residual is neither returned nor thrown, so it is worth not running with warnings suppressed.

• Penalty-based orthogonality (FiniteExcited). 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 explained in the ϵ convention 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, ϵ)
leading_boundary(ψ₀, O, algorithm, environments) -> (ψ, environments, ϵ)

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 Galerkin error (see Returns below)

• maxiter::Int: maximum amount of iterations

• verbosity::Int: display progress information

Returns

• ψ::AbstractMPS: converged leading boundary MPS

• environments: environments corresponding to the converged boundary

• ϵ::Float64: final convergence error upon terminating the algorithm, i.e. the Galerkin error. It is not a truncation error; see find_groundstate 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, ϵ)
approximate(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, ϵ)
approximate!(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, ϵ)
approximate(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ)
approximate!(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ)
approximate((O, ψ), algorithm) -> (ψ′, ϵ)
approximate!(ψ₀, (O, ψ), algorithm) -> (ψ, ϵ)

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 Galerkin error (see Returns below)

• maxiter::Int: maximum amount of iterations

• verbosity::Int: display progress information

• trunc: if supplied, a truncated two-site sweep (DMRG2/IDMRG2) is prepended to refine the bond dimension before the single-site algorithm polishes the result.

Returns

• ψ: the approximated state

• environments: environments corresponding to the result (not returned by Zipup, which uses none)

• ϵ::Float64: an error measure whose meaning depends on the algorithm:

  • for the iterative algorithms (DMRG, DMRG2, IDMRG, IDMRG2, VOMPS) it is the final convergence error, i.e. the Galerkin error compared against tol. This is the same quantity find_groundstate returns, measuring distance from the variational fixed point.

  • for Zipup it is instead a truncation error, i.e. the largest 2-norm of the discarded singular values over all bonds and sweeps. Zipup is a single non-iterative sweep, so there is no convergence measure to report and no tol to compare against.

  The two are not comparable, and a small ϵ means different things in each case. See the manual on the ϵ convention under The error convention.

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

• Jeckelmann. Phys. Rev. B 66 (2002)

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

• Jeckelmann. Phys. Rev. B 66 (2002)

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