Diagonalization

XDiag features several iterative algorithms to perform diagonalization and time evolution, which do not require dense matrix storage of the involved operators. Instead, applications of operators are implemented on-the-fly (i.e. matrix-free) to minimize memory requirements. All routines described here therefore take the operator directly as an OpSum together with the block it acts on, rather than a precomputed matrix.

Ground state

A fundamental property of a quantum system is its ground state energy, which in XDiag can be easily computed using the eigval0 function.

JuliaC++

e0 = eigval0(ops, block)

double e0 = eigval0(ops, block);

Similarly, the ground state itself can be computed using the eig0 function.

JuliaC++

e0, psi0 = eig0(ops, block)

auto [e0, psi0] = eig0(ops, block);

Here, e0 is a double precision real number and psi0 is a State object. These two functions are convenient shortcuts for the underlying Lanczos implementation described next: they set up a Lanczos iteration, converge to the extremal (lowest) eigenvalue, and hide the additional bookkeeping.

The Lanczos interface

The Lanczos algorithm builds up a Krylov space and projects the operator onto a small tridiagonal matrix, whose extremal eigenvalues rapidly approximate those of the full operator. XDiag exposes this machinery directly through the eigvals_lanczos and eigs_lanczos functions, the latter also returning the eigenvectors.

JuliaC++

res = eigs_lanczos(ops, block)
eigenvalues = res.eigenvalues   # Ritz eigenvalue estimates
eigenvectors = res.eigenvectors
alphas = res.alphas             # diagonal of the tridiagonal matrix
betas = res.betas               # off-diagonal of the tridiagonal matrix

auto res = eigs_lanczos(ops, block);
arma::vec eigenvalues = res.eigenvalues;   // Ritz eigenvalue estimates
State eigenvectors = res.eigenvectors;
arma::vec alphas = res.alphas;             // diagonal of the tridiagonal matrix
arma::vec betas = res.betas;               // off-diagonal of the tridiagonal matrix

These routines return a result object containing more than just the eigenvalues: the coefficients alphas and betas of the tridiagonal Lanczos matrix, the number of iterations niterations, and the convergence criterion that was met. They also offer more control over the convergence properties, for example through the precision, max_iterations, and deflation_tol parameters, and they can be started from a user-provided initial State instead of a random one. The eig0 and eigs_lanczos functions perform the Lanczos iteration twice: first to compute the tridiagonal matrix, and in a second run to build the eigenvectors, in order to minimize memory requirements.

A limitation of the Lanczos algorithm is that, converging essentially one vector at a time, it is not well suited to resolve degeneracies: eigenvalues with a multiplicity greater than one typically show up only once, and nearly degenerate excited states can be difficult to separate reliably.

Excited states

To reliably compute several of the lowest eigenstates, including their degeneracies, XDiag provides the eigvals and eigs functions. eigvals computes the neigs algebraically smallest eigenvalues, and eigs additionally returns the corresponding eigenvectors.

JuliaC++

neigs = 3
eigenvalues = eigvals(ops, block, neigs)

int64_t neigs = 3;
arma::vec eigenvalues = eigvals(ops, block, neigs);

---

JuliaC++

eigenvalues, eigenvectors = eigs(ops, block, neigs)

auto [eigenvalues, eigenvectors] = eigs(ops, block, neigs);

Here eigenvalues is an arma::vec (Julia: Vector{Float64}) holding the neigs lowest eigenvalues in ascending order, while eigenvectors is a single State object holding the neigs eigenvectors as its columns.

These functions are based on the LOBPCG algorithm (Locally Optimal Block Preconditioned Conjugate Gradient). In contrast to Lanczos, LOBPCG is a block eigensolver that iterates a whole set of trial vectors simultaneously. This requires more memory, since several vectors have to be kept in storage at once, but in return it reliably resolves excited states and, in particular, their degeneracies.

The LOBPCG interface

For finer control, the eigs_lobpcg function exposes the algorithm directly.

JuliaC++

res = eigs_lobpcg(ops, block, neigs)
eigenvalues = res.eigenvalues     # the neigs lowest eigenvalues
eigenvectors = res.eigenvectors   # corresponding eigenvectors
residuals = res.residual_norms    # final residual norm of each

auto res = eigs_lobpcg(ops, block, neigs);
arma::vec eigenvalues = res.eigenvalues;    // the neigs lowest eigenvalues
State eigenvectors = res.eigenvectors;      // corresponding eigenvectors
arma::vec residuals = res.residual_norms;   // final residual norm of each

It returns a result object collecting, besides the eigenvalues and eigenvectors, additional information about the run: the final residual_norms of every eigenvector, the number of iterations niterations, the convergence criterion that was met, and the full eigenvalue_history and residual_norms_history across iterations, which are useful for monitoring convergence. Beyond neigs, the tolerance tol and the maximal number of iterations max_iterations can be set. The guard parameter enlarges the iterated block to neigs + guard vectors, so that a degenerate multiplet sitting exactly at the neigs-th eigenvalue is still captured with the correct multiplicity.