evolve_lanczos

Computes the exponential of a Hermitian operator \(H\) with an arbitrary real or complex prefactor \(z\) applied to a State \(\vert \psi_0\rangle\),

\[\vert \psi(z) \rangle = e^{z(H - \delta)} \vert \psi_0\rangle.\]

Here, \(\delta\) denotes a real number shifting the spectrum of \(H\). The algorithm implemented is described in the following publication.

On Krylov Subspace Approximations to the Matrix Exponential Operator
Marlis Hochbruck and Christian Lubich
SIAM Journal on Numerical Analysis, Vol. 34, Iss. 5 (1997)
DOI: 10.1137/S0036142995280572

Sources
evolve_lanczos.hpp
evolve_lanczos.cpp
evolve_lanczos.jl

Definition

The method is provided in two variants:

Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.

C++Julia

EvolveLanczosResult
evolve_lanczos(OpSum const &H, State psi, double t, double precision = 1e-12,
               double shift = 0., bool normalize = false,
               int64_t max_iterations = 1000, double deflation_tol = 1e-7);

EvolveLanczosResult
evolve_lanczos(OpSum const &H, State psi, complex z, double precision = 1e-12,
               double shift = 0., bool normalize = false,
               int64_t max_iterations = 1000, double deflation_tol = 1e-7);

evolve_lanczos(H::OpSum, psi::State, t::Float64; precision::Float64 = 1e-12,
               shift::Float64=0.0, normalize::Bool=false,
               max_iterations::Int64 = 1000, deflation_tol::Float64 = 1e-7)

evolve_lanczos(H::OpSum, psi::State, z::ComplexF64; precision::Float64 = 1e-12,
               shift::Float64=0.0, normalize::Bool=false,
               max_iterations::Int64 = 1000, deflation_tol::Float64 = 1e-7)

An inplace variant evolve_lanczos_inplace, where the input state is overwritten and contains the time evolved state upon exit. This version is more memory efficient than evolve_lanczos.

C++Julia

EvolveLanczosInplaceResult
evolve_lanczos_inplace(OpSum const &H, State &psi, double t, 
                       double precision = 1e-12, double shift = 0.,
                       bool normalize = false, int64_t max_iterations = 1000, 
                       double deflation_tol = 1e-7);

EvolveLanczosInplaceResult
evolve_lanczos_inplace(OpSum const &H, State &psi, complex z, 
                       double precision = 1e-12, double shift = 0.,
                       bool normalize = false, int64_t max_iterations = 1000, 
                       double deflation_tol = 1e-7);

evolve_lanczos_inplace(H::OpSum, psi::State, t::Float64; precision::Float64 = 1e-12,
                       shift::Float64=0.0, normalize::Bool=false,
                       max_iterations::Int64 = 1000, deflation_tol::Float64 = 1e-7)

evolve_lanczos_inplace(H::OpSum, psi::State, z::ComplexF64; precision::Float64 = 1e-12,
                       shift::Float64=0.0, normalize::Bool=false,
                       max_iterations::Int64 = 1000, deflation_tol::Float64 = 1e-7)

Parameters

Name	Description	Default
H	OpSum defining the hermitian operator \(H\) for time evolution
psi0	initial State \(\vert \psi_0 \rangle\) of the time evolution
time	time \(\tau\) until which the state is evolved
precision	accuracy of the computed time evolved state \(\vert \psi(t) \rangle\)	1e-12
shift	the offset \(\delta\) when computing \(\vert \psi(t) \rangle = e^{-(H - \delta) \tau} \vert \psi_0\rangle\)	0.0
normalize	flag whether or not the evolved state should be normalized	false
max_iterations	maximum number of Lanczos iterations performed	1000
deflation_tol	tolerance for deflation, i.e. breakdown of Lanczos due to Krylow space exhaustion	1e-7

The parameter shift can be used to turn all eigenvalues of the matrix \(H - \delta \;\textrm{Id}\) positive whenever \(\delta < E_0\), where \(E_0\) denotes the ground state energy of \(H\).

Returns

A struct with the following entries

Entry	Description
alphas	diagonal elements of the Lanczos tridiagonal matrix
betas	off-diagonal elements of the Lanczos tridiagonal matrix
eigenvalues	the computed Ritz eigenvalues of the tridiagonal matrix
niterations	number of iterations performed
criterion	string denoting the reason why the Lanczosalgorithm stopped
state	time-evolved State \(\vert \psi(t)\rangle\) (not defined for inplace variant)

Convergence criterion

The algorithm is estimating the following error,

\[ \varepsilon = \parallel \vert \tilde{\psi}(t)\rangle - e^{z(H - \delta)} \vert \psi_0\rangle \parallel_2, \]

where \(\vert \tilde{\psi}(t) \rangle\) denotes the approximation computed during the algorithm. As the exact solution is not available this error is estimated using the method described by Algorithm 2 in

Expokit: A Software Package for Computing Matrix Exponentials
Roger B. Sidje
ACM Trans. Math. Softw., 24(1):130-156, 1998. (1998)
DOI: 10.1145/285861.285868

Usage Example

C++Julia

int N = 8;
int nup = N / 2;
auto block = Spinhalf(N, nup);

// Define the nearest-neighbor Heisenberg model
auto ops = OpSum();
for (int i=0; i<N; ++i) {
  ops += Op("SdotS", {i, (i+1) % N});
}

// Compute ground state energy
double e0 = eigval0(ops, block);

auto psi0 = product_state(block, {"Up", "Dn", "Up", "Dn", "Up", "Dn", "Up", "Dn"});
double time = 1.0;
double precision = 1e-12;
auto res = evolve_lanczos(ops, psi0, time, precision, e0, true, 500);
XDIAG_SHOW(res.alphas);
XDIAG_SHOW(res.betas);

let
    N = 8
    nup = N ÷ 2
    block = Spinhalf(N, nup)

    # Define the nearest-neighbor Heisenberg model
    ops = OpSum()
    for i in 1:N
        ops += Op("SdotS", [i, mod1(i+1, N)])
    end

    # Compute ground state energy
    e0 = eigval0(ops, block)

    psi0 = product_state(block, ["Up", "Dn", "Up", "Dn", "Up", "Dn", "Up", "Dn"])
    time = 1.0
    res = evolve_lanczos(ops, psi0, time, precision=1e-12, shift=e0, normalize=true)
    @show res.alphas
    @show res.betas
end