Algorithms

EvolutionAlgorithm

the supertype of the methods of evolve and thermalize, Tdvp and ApproxW. An algorithm of one's own is a subtype with a method of TensorMixedStates.evolve_step, which makes one step, see evolve_step.

Tdvp(; krylov = Krylov(), expand_period = 0, hermitianize_period = 0)

time evolution by tdvp, for evolve and thermalize. The bond dimension is enlarged by a global Krylov expansion before any step that starts from a state whose largest bond dimension is below min(5, maxdim), as a product state: tdvp is held to the bonds of the state it starts from, and would leave the error of the first steps, 5e-5 instead of 1e-7 on a Heisenberg chain. A single site, which tdvp does not take, is exponentiated whole.

• krylov: the Krylov exponentiation of each local step, see Krylov
• expand_period: also enlarge it before the first step and then every expand_period steps (default 0, never)
• hermitianize_period: make a mixed state hermitian every hermitianize_period steps (default 0, never)

Examples

Tdvp()
Tdvp(expand_period = 5)
Tdvp(krylov = Krylov(tol = 1e-10))

ApproxW(; order, w = 2, apply_algo = "naive", hermitianize_period = 0)

time evolution by the WI or WII approximations of the exponential of the evolver, combined into an approximation of the given order, for evolve and thermalize.

• order: the order of the approximation, from 1 to 4. In imaginary time, as for thermalize, two of the seven sub steps of order 4 go backwards, which amplifies the high energies by about $e^{0.03\,d\beta\,\Delta E}$, ΔE the width of the spectrum: unstable unless dβ ΔE is small, where orders 1 to 3 only go forwards
• w: 1 or 2 for WI or WII (default 2)
• apply_algo: the algorithm of the product of the state by each MPO, see apply: "naive" (default), exact then truncated, or "densitymatrix", faster at large bond dimensions without conservation (2.4x at 128 on a Heisenberg chain) but with the cutoff eps() it needs, which stops the error near 1e-7 and lets it grow with the number of steps
• hermitianize_period: make a mixed state hermitian every hermitianize_period steps (default 0, never): from order 2, the complex fractions of the step keep a density matrix hermitian only up to the order of the approximation

Examples

ApproxW(order = 2)
ApproxW(order = 4, w = 1)

SearchAlgorithm

the supertype of the methods of ground_state and steady_state, Dmrg. An algorithm of one's own is a subtype with a method of TensorMixedStates.search, see search.

Dmrg(; nsweeps, noise = 0, krylov = Krylov(), tol = 0)

the search of an eigenvector by dmrg, for ground_state and steady_state.

• nsweeps: the number of sweeps, at most
• noise: the noise of each sweep, a number or one value per sweep, continued with its last
• krylov: the Krylov search of each local step, see Krylov
• tol: the search stops when the energy changes by less than tol from one sweep to the next (default 0, never)

Examples

Dmrg(nsweeps = 10)
Dmrg(nsweeps = 20, noise = [1e-6, 1e-7, 0], tol = 1e-10)

ground_state(hamiltonian, state; algo, limits = Limits(), mpo_limits = Limits(),
             orthogonal_to = [], weight = nothing)

the ground state of the Hamiltonian, searched from the pure state by the method algo, see SearchAlgorithm, as the named tuple (energy, state). The Hamiltonian may be given as its PreMPO(state, hamiltonian), prepared once for several searches and keeping its MPO. limits, see Limits, may give one value per sweep, with a default cutoff of 1e-26, and mpo_limits truncates the MPO of the Hamiltonian, see PreMPO (default none). A mixed state is refused: search the pure state and mix it; the steady state of a Lindbladian is given by steady_state.

Given orthogonal_to, eigenstates of the Hamiltonian found before, the search gives the lowest eigenstate orthogonal to them, an excited state, weight times their projectors being added to the Hamiltonian; it must exceed the gap between the state searched and each of them, and is by default 1.1 times a bound of the width of the spectrum, sure but slow to converge, a weight just above the gap converging faster. Its precision is that of the states given, and the penalty is exact only for eigenstates. A warning says when the state found keeps an overlap with one of them. The energy is that of the Hamiltonian alone.

Examples

energy, state = ground_state(H, state; algo = Dmrg(nsweeps = 10),
                             limits = Limits(maxdim = [10, 20, 50]))
e1, excited = ground_state(H, state; algo = Dmrg(nsweeps = 20), orthogonal_to = [ground])

steady_state(lindbladian, state; algo, limits = Limits(), mpo_limits = Limits())

the steady state of the Lindbladian $L$, -im * H plus dissipators, the lowest eigenvector of $L^\dagger L$ searched from the mixed state by the method algo, see SearchAlgorithm, as the named tuple (residual, state), the state of trace one and residual the eigenvalue reached, zero for a steady state. The MPO of $L^\dagger L$, built from that of $L$, has as few channels as its rank across each link. The Lindbladian may be given as its PreMPO(state, lindbladian), prepared once and keeping the MPO of $L^\dagger L$.

• limits: the limits of the sweeps, see Limits, with a default cutoff of 1e-26
• mpo_limits: the truncation of the MPO of $L^\dagger L$, see PreMPO (default none)

The Krylov search of Dmrg takes dim = 8 and maxiter = 3 by default, which the spectrum of $L^\dagger L$, crowded near zero, needs. A Lindbladian may have several steady states, one in each sector of a quantity it conserves for instance, and dmrg returns any combination of them, in general not a density matrix; declaring the quantity strong keeps the search in the sector of the starting state. A warning is given when the state found is not hermitian, its HermiticityDeviation above 1e-6: the search has not converged, or the steady state is not unique.

Examples

value, rho = steady_state(-im * H + D, rho; algo = Dmrg(nsweeps = 10),
                          limits = Limits(maxdim = [10, 20, 50]))

evolve(evolver, duration, state; nsteps | time_step, algo = Tdvp(), limits = Limits(),
       mpo_limits = Limits(), time = 0, step = 1)

the state evolved for duration under evolver, in nsteps steps of duration / nsteps, or in steps of time_step adjusted to the nearest one dividing duration into a whole number of them (none when duration is shorter than half a step), one of the two being given, by the method algo, see EvolutionAlgorithm. The evolver is -im * H for a Hamiltonian H, plus dissipators for a mixed state, or its PreMPO(state, evolver), prepared once for several calls and keeping the MPOs they build. One that depends on time is a vector of terms, each an operator or an operator with its real function of time, [A, B => g] for $A + g(t) B$, $t$ starting at time. step is the number of the first step, for the periods of the methods and for limits, see Limits, which may give one value per step, step k taking the k-th. Its default cutoff is 1e-26 for tdvp and apply_algo = "naive", eps() for "densitymatrix". mpo_limits truncates the MPO of the evolver, see PreMPO (default none).

Examples

evolve(-im * H, 1., state; nsteps = 10, limits = Limits(cutoff = 1e-12, maxdim = 50))
evolve([-im * H0, -im * H1 => t -> cos(t)], 1., state; time_step = 0.05,
       algo = ApproxW(order = 4))

thermalize(hamiltonian, beta, state; nsteps | beta_step, algo = Tdvp(), limits = Limits(),
           mpo_limits = Limits(), step = 1)

the mixed state $\rho$, normalized to trace one, taken to $e^{-\beta H/2} \rho \, e^{-\beta H/2}$ by an evolution in imaginary time in nsteps steps of beta / nsteps, or in steps of beta_step adjusted as the time_step of evolve, one of the two being given, by the method algo, see EvolutionAlgorithm (ApproxW has its limits there), the state normalized to trace one at each step. It is returned as the named tuple (log_trace, state), log_trace being $\log \mathrm{tr}(e^{-\beta H/2} \rho \, e^{-\beta H/2})$ for $\rho$ of trace one, so that the log_trace of successive calls add up.

From "FullyMixed", it gives the thermal state $e^{-\beta H}/Z$, and log_trace is $\log Z - \sum_i \log d_i$, $d_i$ being the dimension of site i. A state commuting with $H$ gives the thermal state restricted to what it describes: fully_mixed(system, N => m) the canonical state of m particles, a product state of populations $e^{\beta\mu n_i}$ the grand canonical state at chemical potential $\mu$. The Hamiltonian may be given as its PreMPO on the mixed state, prepared once and keeping its MPOs, and mpo_limits truncates its MPO, see PreMPO (default none). step is the number of the first step, for the periods of the method and for limits, which may give one value per step; from a later step, the state is taken as a computation left it, of trace one, and not normalized again. A pure state is refused.

Examples

log_trace, rho = thermalize(H, 2.0, State{Mixed}(system, "FullyMixed"); nsteps = 20,
                            limits = Limits(cutoff = 1e-12, maxdim = 64))

Krylov(; dim = nothing, maxiter = nothing, tol = nothing)

the parameters of the Krylov method of each local step: the exponentiation of tdvp, the eigenvector search of dmrg. A field left nothing takes the default of the method, field by field, that of KrylovKit for tdvp (30, 100, 1e-12) and of ITensorMPS for dmrg (3, 1, 1e-14), steady_state taking dim = 8 and maxiter = 3.

• dim: the largest dimension of a Krylov space
• maxiter: the number of Krylov spaces built in turn
• tol: the tolerance, for tdvp relative to the norm of the state, whose step is made on the state of norm one; for tdvp, it sets the number of vectors, dim only bounding it

Examples

Krylov(tol = 1e-10)            # tdvp: fewer vectors, a lower precision
Krylov(dim = 8, maxiter = 3)   # dmrg: a more accurate local step

apply(op, state; limits = Limits())
apply(mpo::MPO, state; limits = Limits(), apply_algo = "naive")

the state with the operator op, placed on sites, applied, or the MPO mpo from make_mpo. A product acts as the operator it denotes, its rightmost factor first, factor by factor and never developed: a factor on at most two sites, or given by its matrix on its sites, is applied as a dense block there, any other through its MPO, which costs less for a sum or a product on more sites, as for a factor too large or laying Jordan-Wigner strings on other sites. On a pure state, a channel of a single Kraus operator K, as Gate(K), is applied as K. On a mixed state, a channel is applied factor by factor in the same way, Gate(K) as K on the kets and its conjugate on the bras, an operator K being Gate(K); a generator L gives L(ρ), a dephasing of more than two sites through the MPO of the gates of its projectors when they have few terms. The state is not normalised.

• limits: the truncation of the bonds a gate of several sites crosses, or of the whole state for an MPO, cutoff = 1e-26 by default, eps() with "densitymatrix", whose density matrices give the squares of the singular values only to eps() of the largest
• apply_algo: for an MPO, the algorithm of its product with the state: "naive", the exact product truncated, or "densitymatrix", faster when it truncates much, but whose error may stay near 1e-8, see Camaño, Epperly and Tropp, arXiv:2504.06475

Examples

apply(controlled(Z)(1, 3) * H(2) * controlled(X)(3, 4), state)
apply(exp(-0.1im * (X ⊗ X + Y ⊗ Y))(2, 5), state)
apply(depolarizing_channel(0.01, 2)(1, 2) * controlled(X)(1, 2), ρ)

set_threading(threading::ThreadingState)
set_threading(mode)
set_threading(system)
set_threading(state)

set how ITensors threads the tensor contractions, and return the ThreadingState it replaces. The threading is given as a ThreadingState, or as one of two modes:

• :dense: no block sparse multithreading, BLAS on its threads, Strided on a single one, as ITensors recommends: started with Julia, Strided has as many threads as Julia, which compete with those of BLAS and can slow down dense contractions considerably. BLAS keeps its threads, except after :blocks, where it gets back those Julia starts it with. runTMS applies this mode by default, and a program calling the functions of TMS directly should start with set_threading(:dense);
• :blocks: block sparse multithreading on the threads of Julia, BLAS and Strided on a single thread, for a system that conserves something. Julia has to be started with several threads, julia -t N: on a single one, everything runs on one core, and ITensors warns.

Given a system or a state, the mode is :blocks when the system conserves something and Julia has several threads, :dense otherwise; which is faster depends on the calculation, see Threads and performance. These settings of the process last until they are changed; the threads of Julia and of its garbage collector, fixed when Julia starts, are not among them.

Examples

set_threading(:dense)
old = set_threading(System(20, Fermion(conserve = N)))   # :blocks, given several threads
set_threading(old)                                       # back to what it was
set_threading(ThreadingState(blas = 2))

ThreadingState(; blas, strided, blocksparse)

how the tensor contractions are threaded, in what a program can change, see threading_settings:

• blas: the threads of BLAS, which runs the products of dense matrices
• strided: the threads of Strided, which ITensors uses for the dense permutations
• blocksparse: whether ITensors multithreads the contractions of block sparse tensors

A field left out takes the value in force, so that ThreadingState() is the threading in force. set_threading applies one, returning the one it replaces.

Examples

ThreadingState()                          # the threading in force
set_threading(ThreadingState(blas = 2))   # BLAS on two threads, the rest as it is

threading_settings()

how the tensor contractions are threaded, as a named tuple:

• julia: the threads of Julia, on which block sparse multithreading runs
• gc: the threads of the garbage collector
• blas, blas_library: the threads of BLAS, which runs the products of dense matrices, and its library, OpenBLAS unless MKL was loaded
• strided: the threads of Strided, which ITensors uses for the dense permutations
• blocksparse: whether ITensors multithreads the contractions of block sparse tensors

julia and gc are fixed when Julia starts, by julia --threads=N --gcthreads=M; the others can be changed at any time, see ThreadingState and set_threading.