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\),
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
The algorithm can be run either on-the-fly (matrix-free) or using a sparse matrix in the compressed-sparse-row format (see CSRMatrix).
Sources
evolve_lanczos.hpp
evolve_lanczos.cpp
evolve_lanczos.jl
Definition
On-the-fly
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.
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 thanevolve_lanczos.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)
Sparse matrix
-
Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.
template <typename idx_t, typename coeff_t> EvolveLanczosResult evolve_lanczos( CSRMatrix<idx_t, coeff_t> const &H, State psi, double tau, double precision = 1e-12, double shift = 0., bool normalize = false, int64_t max_iterations = 1000, double deflation_tol = 1e-7); template <typename idx_t, typename coeff_t> EvolveLanczosResult evolve_lanczos( CSRMatrix<idx_t, coeff_t> const &H, State psi, complex tau, double precision = 1e-12, double shift = 0., bool normalize = false, int64_t max_iterations = 1000, double deflation_tol = 1e-7);evolve_lanczos(H::CSRMatrix, psi::State, t::Float64; precision::Float64 = 1e-12, shift::Float64=0.0, normalize::Bool=false, max_iterations::Int64 = 1000, deflation_tol::Float64 = 1e-7)::EvolveLanczosResult evolve_lanczos(H::CSRMatrix, psi::State, z::ComplexF64; precision::Float64 = 1e-12, shift::Float64=0.0, normalize::Bool=false, max_iterations::Int64 = 1000, deflation_tol::Float64 = 1e-7)::EvolveLanczosResult -
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 thanevolve_lanczos.<typename idx_t, typename coeff_t> EvolveLanczosInplaceResult evolve_lanczos_inplace( CSRMatrix<idx_t, coeff_t> const &H, State &psi, double tau, double precision = 1e-12, double shift = 0., bool normalize = false, int64_t max_iterations = 1000, double deflation_tol = 1e-7); template <typename idx_t, typename coeff_t> EvolveLanczosInplaceResult evolve_lanczos_inplace( CSRMatrix<idx_t, coeff_t> const &H, State &psi, complex tau, double precision = 1e-12, double shift = 0., bool normalize = false, int64_t max_iterations = 1000, double deflation_tol = 1e-7);evolve_lanczos_inplace(H::CSRMatrix, psi::State, t::Float64; precision::Float64 = 1e-12, shift::Float64=0.0, normalize::Bool=false, max_iterations::Int64 = 1000, deflation_tol::Float64 = 1e-7)::EvolveLanczosInplaceResult evolve_lanczos_inplace(H::CSRMatrix, psi::State, z::ComplexF64; precision::Float64 = 1e-12, shift::Float64=0.0, normalize::Bool=false, max_iterations::Int64 = 1000, deflation_tol::Float64 = 1e-7)::EvolveLanczosInplaceResult
Parameters
| Name | Description | Default |
|---|---|---|
| H | OpSum or CSRMatrix 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,
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
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;
// On-the-fly
auto res = evolve_lanczos(ops, psi0, time, precision, e0, true, 500);
XDIAG_SHOW(res.alphas);
XDIAG_SHOW(res.betas);
// sparse matrix
auto spmat = csr_matrix(ops, block);
auto ress = evolve_lanczos(spmat, psi0, time, precision, e0, true, 500);
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
# on-the-fly
res = evolve_lanczos(ops, psi0, time, precision=1e-12, shift=e0, normalize=true)
@show res.alphas
@show res.betas
# sparse matrix
spmat = csr_matrix(ops, block)
res = evolve_lanczos(spmat, psi0, time, precision=1e-12, shift=e0, normalize=true)
end