Skip to content

eigvals_lanczos

Performs an iterative eigenvalue calculation using the Lanczos algorithm. Returns the tridiagonal matrix, eigenvalues, number of iterations and the stopping criterion.

Sources
eigvals_lanczos.hpp
eigvals_lanczos.cpp
eigvals_lanczos.jl

Definition

The Lanczos algorithm can be run in thre distinct ways:

  1. A random intial state \(|\psi_0\rangle = |r\rangle\) with normal distributed entries is used.

    EigvalsLanczosResult
eigvals_lanczos(OpSum const &ops, Block const &block, int64_t neigvals = 1,
                double precision = 1e-12, int64_t max_iterations = 1000,
                double deflation_tol = 1e-7, int64_t random_seed = 42);

    eigvals_lanczos(ops::OpSum, block::Block; neigvals::Int64 = 1, 
                precision::Float64 = 1e-12, max_iterations::Int64 = 1000,
                deflation_tol::Float64 = 1e-7, random_seed::Int64 = 42)

  2. The initial state \(|\psi_0\rangle\) is explicitly specified. 

    EigvalsLanczosResult 
eigvals_lanczos(OpSum const &ops, State psi0, int64_t neigvals = 1,
                double precision = 1e-12, int64_t max_iterations = 1000,
                double deflation_tol = 1e-7);

    eigvals_lanczos(ops::OpSum, psi0::State; neigvals::Int64 = 1, 
                precision::Float64 = 1e-12, max_iterations::Int64 = 1000,
                deflation_tol::Float64 = 1e-7)

    Notice this version copies the initial state, which requires memory but keeps the orginal state intact.

  3. The initial state \(|\psi_0\rangle\) is explicitly specified and overwritten in the process. This version can save memory, but the initial state \(|\psi_0\rangle\) cannot be used later.

    EigvalsLanczosResult 
eigvals_lanczos_inplace(OpSum const &ops, State &psi0, int64_t neigvals = 1,
                        double precision = 1e-12, int64_t max_iterations = 1000,
                        double deflation_tol = 1e-7);

    eigvals_lanczos_inplace(ops::OpSum, psi0::State; neigvals::Int64 = 1, 
                        precision::Float64 = 1e-12, max_iterations::Int64 = 1000,
                        deflation_tol::Float64 = 1e-7)

Parameters

Name Description Default
ops OpSum defining the bonds of the operator
block block on which the operator is defined
psi0 Initial State from which the Lanczos algorithm is started
neigvals number \(k\) of eigenvalue to converge 1
precision accuracy of the computed ground state 1e-12
max_iterations maximum number of iterations 1000
deflation_tol tolerance for deflation, i.e. breakdown of Lanczos due to Krylow space exhaustion 1e-7
random_seed random seed for setting up the initial vector 42

Returns

A struct of type EigvalsLanczosResult with the following entries.

Entry Description
alphas diagonal elements of the tridiagonal matrix
betas off-diagonal elements of the tridiagonal matrix
eigenvalues the computed Ritz eigenvalues of the tridiagonal matrix
niterations number of iterations performed
criterion string denoting the reason why the algorithm stopped

Convergence criterion

The algorithm terminates if the \(k\)-th (\(k\) is the argument neigvals) approximate eigenvalue changes only by a fraction smaller than \(\epsilon\) (\(k\) is the argument precision), i.e.

\[ (\tilde{e}_k^{(n)} - \tilde{e}_k^{(n-1)}) / \tilde{e}_k^{(n)} < \epsilon.\]

Here, \(\tilde{e}_k^{(n)}\) denotes the Lanczos approximation to the \(k\)-th eigenvalue after \(n\) iterations.

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 += "J" * Op("SdotS", {i, (i+1) % N});
}
ops["J"] = 1.0;

// With random intial state
auto res = eigvals_lanczos(ops, block);
XDIAG_SHOW(res.alphas);
XDIAG_SHOW(res.betas);
XDIAG_SHOW(res.eigenvalues);

// With specific initial state
auto psi0 = product_state(block, {"Up", "Dn", "Up", "Dn", "Up", "Dn", "Up", "Dn"});
auto res2 = eigvals_lanczos(ops, psi0);
XDIAG_SHOW(res.alphas);
XDIAG_SHOW(res.betas);
XDIAG_SHOW(res.eigenvalues);

let
    N = 8
    nup = N รท 2
    block = Spinhalf(N, nup)

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

    # With random intial state
    res = eigvals_lanczos(ops, block)
    @show res.alphas
    @show res.betas
    @show res.eigenvalues

    # With specific initial state
    psi0 = product_state(block, ["Up", "Dn", "Up", "Dn", "Up", "Dn", "Up", "Dn"])
    res2 = eigvals_lanczos(ops, psi0)
    @show res.alphas
    @show res.betas
    @show res.eigenvalues
end