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.truncation_errors with what its final sweep discarded at each bond. See AlgorithmInfo for the full vocabulary, and Ground state accuracy 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

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

• 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, see Ground state accuracy

• 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, see Ground state accuracy

• 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 on the change of the center bond tensor over a sweep, see Ground state accuracy

• 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 on the change of the center bond tensor over a sweep, see Ground state accuracy

• 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 on the Galerkin error, see Ground state accuracy

• 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 on the norm of the gradient, see Ground state accuracy

• 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, 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 finite-system step sweeps over every bond twice, and info.truncation_errors[i] is the truncation error of the last cut made at bond i, i.e. of the second half-sweep. The entries are non-zero only for algorithms that truncate (TDVP2, BUG with a trunc, and TDVP with a bond expansion), and exactly 0 for a step that happened to discard nothing. Infinite one-site TDVP never truncates and reports no truncation_errors at all. Neither case means the step was exact, but rather that this particular source of error is either absent or idle.

See AlgorithmInfo 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: numiter, the number of steps taken, and for an algorithm that truncates truncation_errors, holding the per-bond truncation errors of every step as reported by timestep, one vector per step. See AlgorithmInfo and Time evolution accuracy in the manual for what this does not measure.

The largest truncation error of each step is logged at verbosity ≥ 3, and that of the whole evolution at verbosity ≥ 2.

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 the squared norm drops by precisely the weight discarded by all cuts of the step, whereas in imaginary time the norm also carries the physical decay of the weight. Without trunc nothing is discarded at all and the norm is conserved exactly in real time.

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

Fields

• 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 all bond cuts of the step, whereas in imaginary time the norm also carries the physical decay of the weight. Pass normalize = true to timestep/time_evolve to renormalize after every half-sweep instead.

Warning

The truncation_errors reported by timestep are those of the cuts in step (i). Step (iii) augments the basis without truncating, so within a sweep the bond dimension temporarily exceeds trunc.

Tip

Pass a TimerOutputs.TimerOutput as timeroutput to timestep/timestep! to obtain a breakdown of the time spent in the three steps above (cut_bond, AC_integrate, augment) 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

MPSKit.ClusterExpansion Type

struct ClusterExpansion <: MPSKit.Algorithm

Algorithm for constructing a finite or infinite time evolution MPO using the nonperturbative cluster expansion. Exact exponentials are matched on clusters of up to N sites. For nearest-neighbor Hamiltonians the error is O(dt^N).

Fields

• N::Int64: maximum cluster size (positive integer, default: 2)

• tol::Float64: relative singular-value cutoff for shared SVD support (0 < tol < 1, default: 1.0e-12)

Algorithm

Exact evolution within clusters of up to N sites is included to all orders in dt, together with products of disjoint clusters. The remaining error comes from connected contributions extending beyond N sites: for nearest-neighbor Hamiltonians the one-step error is O(dt^N). Infinite Hamiltonians include clusters crossing the unit-cell boundary; finite chains become exact up to roundoff when N reaches the chain length.

Compared with TaylorCluster, which matches a specified Taylor order, ClusterExpansion controls the spatial extent of the exact evolution. Its N is a cluster size, not a Taylor order. Resumming smaller clusters can reduce the error prefactor and permit larger time steps. Compared with the compact, local construction of WII, increasing N systematically incorporates larger clusters, at the cost of exponentially growing cluster exponentials and MPO bond dimensions. For full-rank uniform physical spaces of dimension d, the latter scale as 1 + d^2 + ⋯ + d^(2 floor(N/2)); mixed spaces or rank-deficient residuals can require extra channels.

This is useful for evolution with fewer, larger steps and for growing the MPS bond dimension from weakly entangled states. Larger steps can also increase compression costs, so check convergence and runtime by varying dt and N. MPS compression introduces a separate error, and the approximate MPO need not be unitary. tol only identifies shared SVD support: smaller singular components are retained in complementary channels rather than discarded.

See also

Used as the algorithm argument of make_time_mpo.

References

• Vanhecke et al. Phys. Rev. A 103 (2021)

source

See also Time evolution accuracy.

Nonperturbative cluster expansion

ClusterExpansion constructs evolution MPOs for finite and infinite Hamiltonians by matching exact exponentials on clusters of up to N sites (Vanhecke et al., 2021). Unlike TaylorCluster, N specifies the cluster size rather than a Taylor order. For nearest-neighbor interactions, the one-step error is O(dt^N), a finite chain is matched up to roundoff when N reaches its length.

U = make_time_mpo(H, dt, ClusterExpansion(; N = 3))
ψ, envs, info = approximate(ψ, (U, ψ); trunc = truncrank(32))

Matching small clusters to all orders can permit larger time steps and help grow the bond dimension from weakly entangled states. The construction costs more than the compact WII approximation, and increasing N requires exponentially larger cluster exponentials and MPO bonds. Check convergence by varying dt and N, as well as the MPS compression settings.

Excitations

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

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

The following examples use the shared model constructors, built from TensorKitTensors operators and MPSKit Hamiltonians.

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(; 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; 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, B} <: 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

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

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 instead find the ground state of the Hamiltonian where we project out previously obtained eigenstates, i.e. by projecting out       . This is also supported by calling

# Model parameters
g = 10.0
L = 16
H = transverse_field_ising(; 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(; 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.

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, linear_solve, 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 residual projected onto the local variational space. For an eigenproblem this is the local update projected onto the orthogonal complement of the current tensor, as reported by DMRG, DMRG2, VUMPS and VOMPS when solving for a state. For a linear system it is the local residual ‖(a₀ + a₁·A)·x − b‖ relative to ‖b‖, as reported by DMRGSolve and DMRGSolve2 in linear_solve.

• 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

• truncation_errors: the truncation error of every bond, i.e. the 2-norm of the singular values discarded by the most recent factorisation of that bond. Entry i is the bond between sites i and i + 1, so a finite MPS of length L has L - 1 entries and an infinite unit cell of length L has L, the last one being the bond across the unit cell boundary. A Multiline algorithm reports a matrix instead, indexed as [row, bond].

A sweep that cuts a bond more than once keeps only the last cut, which is what the returned state still throws away there. An algorithm that truncates reports this entry even when it discarded nothing (all entries zero), whereas the entry being absent means the algorithm never truncates. Neither says the result is exact. See the manual on Errors and accuracy for what is not measured here.

source

The rest of this section explains what quantities can be reported by the algorithms, and how these can be interpreted.

Ground state accuracy

find_groundstate, leading_boundary and the iterative approximate algorithms stop once a convergence measure drops below their tol, and report that measure together with a converged flag.

Different measures for different algorithms. Which quantity is monitored depends on the algorithm, and these quantities cannot be compared to one another. Each is therefore stored under its own name: galerkin for DMRG, DMRG2, VUMPS and VOMPS, gradientnorm for GradientGrassmann, bondresidual for IDMRG and IDMRG2, and localchange for DMRG and DMRG2 inside approximate. convergence_measure returns whichever of them is present, for code that only wants the number. For the same reason, a tol that works well for one algorithm need not work well for another.

Interaction with a changing bond dimension. An algorithm that truncates changes the state by the discarded weight in every sweep, so its convergence measure cannot drop below the level set by what is discarded. Rather than running into maxiter for a tol that cannot be reached, DMRG and DMRG2 consider themselves converged once galerkin ≤ max(tol, maximum(truncation_errors)), which reduces to galerkin ≤ tol when nothing is truncated. The reported truncation_errors are those of the final sweep, i.e. what the returned state still discards.

Convergence is not accuracy. For algorithms that remain on a fixed bond dimension MPS manifold, one can get the ground state on that manifold to converge to machine precision. This however might still be far from the true ground state, whenever the bond dimension is not high enough. Growing the bond dimension is the job of the two-site algorithms (DMRG2, IDMRG2) or of a bond expansion (DMRG with an alg_expand, or an expanding alg_gauge such as DMRG3S); see also changebonds. Likewise, none of these measures is an error bar on an observable. An independent check is the energy variance,  , which vanishes for an exact eigenstate, as is comparing results across bond dimensions.

Time evolution accuracy

A time evolution has three distinct sources of error (Paeckel et al.):

• The truncation error, from discarding singular values to keep the bond dimension bounded. This is the only one that is reported, as the truncation_errors of timestep and time_evolve. It is non-zero for TDVP2, for BUG with a trunc, and for TDVP with a bond expansion.

• The projection error, from restricting the evolution to MPS of the current bond dimension. Measuring it costs an extra two-site effective-Hamiltonian application per bond, so it is not reported. A fixed bond dimension turns entanglement growth into projection error, which is what a bond expansion (CBE) exists to reduce (Li et al.).

• The time-step error, from splitting each step into a sequence of local updates, similar to a Trotter error. It decreases with dt, and can be estimated by comparing one step of dt against two steps of dt / 2.

Each of these has its own control: dt for the time-step error, trunc for the truncation error, and a bond-adaptive scheme (TDVP2, BUG, or TDVP with alg_expand) for the projection error whenever entanglement grows.

These controls are not independent. With a truncation strategy trunc::TruncationStrategy, every step can discard weight, however small dt is, and reaching a fixed final time T takes T / dt steps. Halving dt therefore lowers the time-step error, but doubles the number of truncations, so the accumulated truncation error can grow. For a given trunc there is an optimal dt, below which taking smaller steps makes the result worse rather than better.

Excitation accuracy

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

• Inherited ground state error. Every method builds on the ground state you supply and treats it as exact. Since a gap is a difference of two large energies, that error propagates straight into it and is typically the limiting factor.

• Ansatz limitation. QuasiparticleAnsatz varies over the single-quasiparticle tangent space on top of a fixed ground state, so it is variational within that space and suited to isolated quasiparticle branches. Its error is bounded exponentially in the support of the local operator, at a rate set by the gaps below and above the targeted eigenvalue (Haegeman et al.).

• Eigensolver convergence. The eigenvalue problem is solved with KrylovKit, and a run that fails to converge num states emits a warning carrying the residual when the verbosity is high enough. That residual is neither returned nor thrown, so it is worth not suppressing warnings. Nearly degenerate levels are the ones most likely to come back unconverged.

• Penalty-based orthogonality. FiniteExcited minimises   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 energy downwards — invisibly, since the reported value is the expectation value of the bare H. Raising weight suppresses the bias at the cost of stretching the spectrum and slowing the eigensolver.

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

changebonds

Many of the previously mentioned algorithms do not possess a way to dynamically change to 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 the truncation_errors of an AlgorithmInfo. 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. IDMRG2 additionally reports truncation_errors, what its final sweep discarded at each bond. See AlgorithmInfo, find_groundstate, and the manual on Ground state accuracy.

source

The associated eigenvalue, i.e. the partition function density, is obtained with leading_eigenvalue. For a single-line InfiniteMPS/InfiniteMPO pair this is equivalent to expectation_value.

MPSKit.leading_eigenvalue Function

leading_eigenvalue(ψ, O, [environments]) -> λ

Eigenvalue of the transfer operator O for the boundary MPS ψ, accumulated over one unit cell of the network. See the manual at One fixed point, many lines on how to interpret the role of the lines in a MultilineMPS and how the eigenvalue is accumulated over them.

Arguments

• ψ::Union{InfiniteMPS, MultilineMPS}: the boundary MPS

• O::Union{InfiniteMPO, MultilineMPO}: the transfer operator

• environments: the environments to use, calculated if not given

Returns

• λ::Number: the eigenvalue for one unit cell. Complex in general, since a transfer operator need not be Hermitian.

This is not an expectation value

For a MultilineMPO with more than one row this is not an overlap. The contraction pairs row i of the operator with line i of ψ as the ket and line i + 1 as the bra, so no ⟨ϕ|O|ϕ⟩ is computed anywhere. A genuine expectation value of a two-dimensional operator additionally requires the fixed point on the other side of the network, which is a separate leading_boundary run, contracted through mixed environments.

For a single line the bra and the ket coincide and the two notions agree, which is why expectation_value forwards here for InfiniteMPS/InfiniteMPO.

See also

leading_boundary, expectation_value

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.

  • IDMRG2 fills truncation_errors with what its final sweep discarded at each bond, i.e. what the returned state is still throwing away rather than what the early, unconverged sweeps did.

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

  See AlgorithmInfo for the full list and Ground state accuracy 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

linear_solve

Solving a linear system      for an MPS , with an MPO and an MPS. The shift     follows the convention of KrylovKit.linsolve and is applied implicitly, so no shifted operator is ever formed. The sweeps are DMRG-like: DMRGSolve keeps the bond dimension of the initial guess fixed, while DMRGSolve2 updates two sites at a time and is therefore rank-adaptive. Each is parameterized by how the local subproblem is posed — Galerkin solves the effective system directly, LeastSquares minimizes       through the normal equations.

The resolvent is the canonical application:   is linear_solve(ψ₀, H, ψ₀; a₀ = z, a₁ = -1).

MPSKit.linear_solve Function

linear_solve(x₀, A, b, [algorithm], [environments]; kwargs...) -> (x, environments, info)
linear_solve!(x₀, A, b, [algorithm], [environments]; kwargs...) -> (x, environments, info)

Solve the linear system (a₀ + a₁·A)·x = b for the MPS x, where A is an MPO / MPOHamiltonian and b is an MPS. The a₀ + a₁·A shift follows KrylovKit.linsolve and is applied implicitly; the defaults a₀ = 0, a₁ = 1 give the plain system A·x = b. linear_solve! solves in-place, overwriting x₀.

Currently only finite MPS are supported.

Arguments

• x₀::AbstractMPS: initial guess (and, for linear_solve!, the state overwritten with the solution)

• A: operator of the linear system (MPO / MPOHamiltonian)

• b::AbstractMPS: right-hand side

• algorithm: linear-solve algorithm, see DMRGSolve and DMRGSolve2

• [environments]: MPS environment manager for the operator sandwich ⟨x|A|x⟩

Keyword Arguments

The keyword-based call (no explicit algorithm) selects an algorithm from the structure flags:

• tol::Float64: tolerance for convergence criterium

• maxiter::Int: maximum amount of iterations

• verbosity::Int: display progress information

• a₀, a₁: shift/scale scalars of the system (a₀ + a₁·A)·x = b

• isposdef::Bool: declare (a₀ + a₁·A) positive-definite, so the local solver is CG instead of GMRES

• trunc: if supplied, a truncated two-site sweep (DMRGSolve2) is prepended to adapt the bond dimension before the single-site algorithm polishes the result

Returns

• x::AbstractMPS: the (bond-dimension-limited) solution

• environments: operator environments corresponding to x

• info::AlgorithmInfo: how the algorithm terminated. info.galerkin is the Galerkin error, the largest local residual ‖(a₀ + a₁·A)·x − b‖ of the last sweep relative to ‖b‖, and info.converged whether it met the stopping test. DMRGSolve2 also reports info.truncation_errors. See AlgorithmInfo.

source

MPSKit.linear_solve! Function

linear_solve(x₀, A, b, [algorithm], [environments]; kwargs...) -> (x, environments, info)
linear_solve!(x₀, A, b, [algorithm], [environments]; kwargs...) -> (x, environments, info)

Solve the linear system (a₀ + a₁·A)·x = b for the MPS x, where A is an MPO / MPOHamiltonian and b is an MPS. The a₀ + a₁·A shift follows KrylovKit.linsolve and is applied implicitly; the defaults a₀ = 0, a₁ = 1 give the plain system A·x = b. linear_solve! solves in-place, overwriting x₀.

Currently only finite MPS are supported.

Arguments

• x₀::AbstractMPS: initial guess (and, for linear_solve!, the state overwritten with the solution)

• A: operator of the linear system (MPO / MPOHamiltonian)

• b::AbstractMPS: right-hand side

• algorithm: linear-solve algorithm, see DMRGSolve and DMRGSolve2

• [environments]: MPS environment manager for the operator sandwich ⟨x|A|x⟩

Keyword Arguments

The keyword-based call (no explicit algorithm) selects an algorithm from the structure flags:

• tol::Float64: tolerance for convergence criterium

• maxiter::Int: maximum amount of iterations

• verbosity::Int: display progress information

• a₀, a₁: shift/scale scalars of the system (a₀ + a₁·A)·x = b

• isposdef::Bool: declare (a₀ + a₁·A) positive-definite, so the local solver is CG instead of GMRES

• trunc: if supplied, a truncated two-site sweep (DMRGSolve2) is prepended to adapt the bond dimension before the single-site algorithm polishes the result

Returns

• x::AbstractMPS: the (bond-dimension-limited) solution

• environments: operator environments corresponding to x

• info::AlgorithmInfo: how the algorithm terminated. info.galerkin is the Galerkin error, the largest local residual ‖(a₀ + a₁·A)·x − b‖ of the last sweep relative to ‖b‖, and info.converged whether it met the stopping test. DMRGSolve2 also reports info.truncation_errors. See AlgorithmInfo.

source

MPSKit.DMRGSolve Type

struct DMRGSolve{L<:MPSKit.LinsolveFormulation, S, F, B} <: MPSKit.Algorithm

Single-site DMRG-style sweeping algorithm for the MPS linear solver linear_solve. It keeps the bond dimension of the initial guess fixed and, at each site, solves the local linear problem selected by formulation with the local Krylov solver.

Fields

• formulation::MPSKit.LinsolveFormulation: formulation of the local subproblem, either Galerkin or LeastSquares

• solver::Any: local linear solver; a plain KrylovKit solver, or one wrapped in DynamicTol for per-bond adaptive tolerances (the default)

• tol::Float64: tolerance for convergence criterium

• maxiter::Int64: maximal amount of iterations

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

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

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

source

MPSKit.DMRGSolve2 Type

struct DMRGSolve2{L<:MPSKit.LinsolveFormulation, S, G, F, B} <: MPSKit.Algorithm

Two-site DMRG-style sweeping algorithm for the MPS linear solver linear_solve. Each bond update solves the local two-site linear problem selected by formulation and truncates the enlarged bond back down with alg_gauge (a truncated SVD built from trunc), making the bond dimension adaptive.

Fields

• formulation::MPSKit.LinsolveFormulation: formulation of the local subproblem, either Galerkin or LeastSquares

• solver::Any: local linear solver; a plain KrylovKit solver, or one wrapped in DynamicTol for per-bond adaptive tolerances

• tol::Float64: tolerance for convergence criterium

• maxiter::Int64: maximal amount of iterations

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

• alg_gauge::Any: truncated SVD used for the post-update gauge

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

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

source

MPSKit.Galerkin Type

struct Galerkin <: MPSKit.LinsolveFormulation

Local formulation imposing the Galerkin condition: the residual (a₀ + a₁·A)·x − b is projected orthogonal to the local tangent space, yielding the effective system (a₀ + a₁·A_eff)·x = b_eff which is solved with the configured Krylov solver. Use a CG solver for hermitian positive-definite A, and GMRES/BiCGStab for general (non-hermitian or indefinite) A.

source

MPSKit.LeastSquares Type

struct LeastSquares <: MPSKit.LinsolveFormulation

Local formulation that minimizes ‖(a₀ + a₁·A)·x − b‖², i.e. solves the normal equations M†M·x = M†b with M = a₀ + a₁·A, built from squared operator environments. The normal operator is unconditionally positive-definite (solve it with CG) at the cost of squaring the condition number.

Currently assumes A is hermitian (the resolvent / dynamical-DMRG case), for which M†M = |a₀|² + 2·Re(ā₀·a₁)·A + |a₁|²·A² and M†b = ā₀·b + ā₁·A·b. This is a generalization of the Jeckelmann dynamical-DMRG functional.

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  . This is a thin wrapper around linear_solve, see there for the sweep, the convergence criterion and the local solver options.

MPSKit.propagator Function

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

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

Returns

• g: approximation of the propagator matrix element   

• ψ: for NaiveInvert, the MPS approximation of  ; for Jeckelmann, the vector its functional optimizes,    , i.e. the imaginary part of the propagator, from which g is reconstructed.

init is used as the initial guess and is left untouched. The underlying variational problem is solved with linear_solve; for full control over the sweep, call that directly, e.g. linear_solve(ψ₀, H, ψ₀; a₀ = z, a₁ = -1) for the NaiveInvert flavour.

source

MPSKit.DynamicalDMRG Type

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

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

This is a thin wrapper around linear_solve: the sweep, the convergence criterion and the adaptive local tolerances are those of DMRGSolve / DMRGSolve2.

Fields

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

• solver::Any: local linear solver; a plain KrylovKit solver, or one wrapped in DynamicTol for per-bond adaptive tolerances (the default)

• tol::Float64: tolerance for convergence criterium, measured as the relative residual of the linear system

• maxiter::Int64: maximal amount of iterations

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

• trunc::Any: if supplied, a truncated two-site sweep (DMRGSolve2) is prepended to adapt the bond dimension

• 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

The solution of that equation is the imaginary part of the propagator; together with equation (11) from that same paper it determines the full   . Because of that reconstruction step this flavour requires   .

Returns the approximation of    and the minimizer of the functional above (not   itself).

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