Exact Diagonalization

The module Rimu.ExactDiagonalization provides a framework for exact diagonalization of quantum many-body systems defined by an AbstractHamiltonian type.

The main usage is through defining an ExactDiagonalizationProblem and solving it with the solve function. The module provides a unified interface for accessing different solver algorithms, which make use of solvers provided by external packages.

Exports

• ExactDiagonalizationProblem

• BasisSetRepresentation

• build_basis

• KrylovKitSolver

• LinearAlgebraSolver

• ArpackSolver

• LOBPCGSolver

ExactDiagonalizationProblem(hamiltonian::AbstractHamiltonian, [v0]; kwargs...)

Defines an exact diagonalization problem with an AbstractHamiltonian hamiltonian. Optionally, a starting vector of type AbstractDVec, or a single address or a collection of addresses can be passed as v0.

ExactDiagonalizationProblems can be solved with solve.

Keyword arguments

• algorithm=LinearAlgebraSolver(): The algorithm to use for solving the problem. The algorithm can also be specified as the second positional argument in the init function.
• Optional keyword arguments will be passed on to the init and solve functions.

Algorithms

• LinearAlgebraSolver(): An algorithm for solving the problem using the dense-matrix eigensolver from the LinearAlgebra standard library (eventually using LAPACK). Only suitable for small matrices.
• KrylovKitSolver(matrix_free=true): An algorithm for finding a few eigenvalues and vectors. With matrix_free=true the problem is solved without instatiating a matrix. This is suitable for large dimensions. With matrix_free=false the problem is solved after instantiating a sparse matrix. This is faster if sufficient memory is available. Requires using KrylovKit.
• ArpackSolver(): An algorithm for solving the problem after instantiating a sparse matrix and using the Arpack Fortran library. Requires using Arpack.
• LOBPCGSolver(): An algorithm for solving the problem after instantiating a sparse matrix using the LOBPCG method. Requires using IterativeSolvers.

Keyword arguments for matrix-based algorithms (also accepted by init)

See BasisSetRepresentation for more information.

• sizelim: The maximum size of the basis set representation. The default is 10^6 for sparse matrices and 10^5 for dense matrices.
• cutoff: A cutoff value for the basis set representation.
• filter: A filter function for the basis set representation.
• max_depth = Inf: Limit the depth when building the matrix.
• minimum_size = Inf: Stop building the matrix after this size is reached.
• nnzs = 0: A hint for the number of non-zero elements in the basis set representation. Setting a non-zero value can speed up the computation.
• col_hint = 0: A hint for the number of columns in the basis set representation.
• sort = false: Whether to sort the basis set representation.

Keyword arguments for iterative algorithms (also accepted by solve)

• verbose = false: Whether to print additional information.
• abstol = nothing: The absolute tolerance for the solver. If nothing, the solver chooses a default value.
• howmany = 1: The minimum number of eigenvalues to compute.
• which = :SR: Whether to compute the largest or smallest eigenvalues.
• maxiters = nothing: The maximum number of iterations for the solver. If nothing, the solver chooses a default value.

Solving an ExactDiagonalizationProblem

The solve function can be called directly on an ExactDiagonalizationProblem to solve it. Alternatively, the init function can be used to initialize a solver, which can then be solved with solve. The solve function returns a result type with the eigenvalues, eigenvectors, and convergence information.

Result type

The result type for the solve function is determined by the algorithm used. It has the following fields:

• values::Vector: The eigenvalues.
• vectors::Vector{<:AbstractDVec}: The eigenvectors.
• success::Bool: A boolean flag indicating whether the solver was successful.
• info: Convergence information.
• algorithm: The algorithm used for the computation.
• problem: The ExactDiagonalizationProblem that was solved.
• Additional fields may be present depending on the algorithm used.

Iterating the result type will yield the eigenvalues, eigenvectors, and a boolean flag success in that order.

Examples

julia> p = ExactDiagonalizationProblem(HubbardReal1D(BoseFS(1,1,1)))
ExactDiagonalizationProblem(
  HubbardReal1D(fs"|1 1 1⟩"; u=1.0, t=1.0),
  nothing;
  NamedTuple()...
)

julia> result = solve(p) # convert to dense matrix and solve with LinearAlgebra.eigen
EDResult for algorithm LinearAlgebraSolver() with 10 eigenvalue(s),
  values = [-5.09593, -1.51882, -1.51882, 1.55611, 1.6093, 1.6093, 4.0, 4.53982, 4.90952, 4.90952]
  Convergence info: "Dense matrix eigensolver solution from `LinearAlgebra.eigen`", with howmany = 10 eigenvalues requested.
  success = true.

julia> using KrylovKit # an external package has to be installed and loaded

julia> s = init(p; algorithm = KrylovKitSolver(true)) # solve without building a matrix
IterativeEDSolver
 with algorithm KrylovKitSolver(matrix_free = true,) for hamiltonian = HubbardReal1D(fs"|1 1 1⟩"; u=1.0, t=1.0),
  kwargs = NamedTuple()
)

julia> values, vectors, success = solve(s);

julia> result.values[1] ≈ values[1]
true

See also solve(::ExactDiagonalizationProblem), init(::ExactDiagonalizationProblem), KrylovKitSolver, ArpackSolver, LinearAlgebraSolver.

solve(p::ExactDiagonalizationProblem, [algorithm]; kwargs...)

Solve an ExactDiagonalizationProblem p directly. Optionally specify an algorithm. Returns a result type with the eigenvalues, eigenvectors, and convergence information.

For a description of the keyword arguments, see the documentation for ExactDiagonalizationProblem.

See also solve(::ProjectorMonteCarloProblem).

init(p::ExactDiagonalizationProblem, [algorithm]; kwargs...)

Initialize a solver for an ExactDiagonalizationProblem p with an optional algorithm. Returns a solver instance that can be solved with solve.

For a description of the keyword arguments, see the documentation for ExactDiagonalizationProblem.

KrylovKitSolver(matrix_free::Bool; kwargs...)
KrylovKitSolver(; matrix_free = false, kwargs...)

Algorithm for solving a large ExactDiagonalizationProblem to find a few eigenvalues and vectors using the KrylovKit.jl package. The Lanczos method is used for hermitian matrices, and the Arnoldi method is used for non-hermitian matrices.

Arguments

• matrix_free = false: Whether to use a matrix-free algorithm. If false, a sparse matrix will be instantiated. This is typically faster and recommended for small matrices, but requires more memory. If true, the matrix is not instantiated, which is useful for large matrices that would not fit into memory. The calculation will parallelise using threading if available by making use of LinearMap.
• kwargs: Additional keyword arguments are passed on to the function KrylovKit.eigsolve().

See also ExactDiagonalizationProblem, solve, ArpackSolver, LOBPCGSolver.

LinearAlgebraSolver(; kwargs...)

Algorithm for solving an ExactDiagonalizationProblem using the dense-matrix eigensolver from the LinearAlgebra standard library. This is only suitable for small matrices.

The kwargs are passed on to function LinearAlgebra.eigen.

Keyword arguments

• permute = true: Whether to permute the matrix before diagonalization.
• scale = true: Whether to scale the matrix before diagonalization.
• sortby: The sorting order for the eigenvalues.

See also ExactDiagonalizationProblem, solve.

ArpackSolver(; kwargs...)

Algorithm for solving an ExactDiagonalizationProblem after instantiating a sparse matrix. It uses the Lanzcos method for hermitian problems, and the Arnoldi method for non-hermitian problems, using the Arpack Fortran library. This is faster than KrylovKitSolver(; matrix_free=true), but it requires more memory and will only be useful if the matrix fits into memory.

Arguments

• matrix_free = false: Whether to use a matrix-free algorithm. If false, a sparse matrix will be instantiated. This is typically faster and recommended for small matrices, but requires more memory. If true, the matrix is not instantiated, which is useful for large matrices that would not fit into memory. The calculation will parallelise using threading if available by making use of LinearMap.
• Additional kwargs are passed on to the function Arpack.eigs().

See also ExactDiagonalizationProblem, solve, KrylovKitSolver, LOBPCGSolver.

LOBPCGSolver(; kwargs...)

The Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG). Algorithm for solving an ExactDiagonalizationProblem after instantiating a sparse matrix.

LOBPCG is not suitable for non-hermitian eigenvalue problems.

Arguments

• matrix_free = false: Whether to use a matrix-free algorithm. If false, a sparse matrix will be instantiated. This is typically faster and recommended for small matrices, but requires more memory. If true, the matrix is not instantiated, which is useful for large matrices that would not fit into memory. The calculation will parallelise using threading if available by making use of LinearMap.
• Additional kwargs are passed on to the function

IterativeSolvers.lobpcg().

See also ExactDiagonalizationProblem, solve, KrylovKitSolver, ArpackSolver.

BasisSetRepresentation(
    hamiltonian::AbstractOperator, addr=starting_address(hamiltonian);
    sizelim=10^7, cutoff, filter, max_depth, minimum_size, sort=false, kwargs...
)
BasisSetRepresentation(hamiltonian::AbstractOperator, addresses::AbstractVector; kwargs...)

Eagerly construct the basis set representation of the operator hamiltonian with all addresses reachable from addr. Instead of a single address, a vector of addresses can be passed. The second argument is mandatory when the operator is not an AbstractHamiltonian.

An ArgumentError is thrown if dimension(hamiltonian) > sizelim in order to prevent memory overflow. Set sizelim = Inf in order to disable this behaviour.

Providing the number nnzs of expected calculated matrix elements and col_hint for the estimated number of nonzero off-diagonal matrix elements in each matrix column may improve performance.

Providing an energy cutoff will skip the columns and rows with diagonal elements greater than cutoff. Alternatively, an arbitrary filter function can be used instead. Addresses passed as arguments are not filtered. To generate the matrix truncated to the subspace spanned by the addresses, use filter = Returns(false).

Providing a max_depth will limit the size of the matrix and basis by only visiting addresses that are connected to the starting_address through max_depth hops through the Hamiltonian. Similarly, providing minimum_size will stop the bulding process after the basis reaches a length of at least minimum_size.

Setting sort to true will sort the matrix rows and columns. This is useful when the order of the columns matters, e.g. when comparing matrices. Any additional keyword arguments are passed on to Base.sortperm.

The order of the returned basis and matrix rows and columns is arbitrary and
non-deterministic. Use `sort=true` if the ordering matters.

Fields

• sparse_matrix: sparse matrix representing hamiltonian in the basis basis
• basis: vector of addresses
• hamiltonian: the Hamiltonian hamiltonian

Example

julia> hamiltonian = HubbardReal1D(BoseFS(1,0,0));

julia> bsr = BasisSetRepresentation(hamiltonian)
BasisSetRepresentation(HubbardReal1D(fs"|1 0 0⟩"; u=1.0, t=1.0)) with dimension 3 and 6 stored entries:3×3 SparseArrays.SparseMatrixCSC{Float64, Int32} with 6 stored entries:
   ⋅   -1.0  -1.0
 -1.0    ⋅   -1.0
 -1.0  -1.0    ⋅

julia> BasisSetRepresentation(hamiltonian, bsr.basis[1:2]; filter = Returns(false)) # passing addresses and truncating
BasisSetRepresentation(HubbardReal1D(fs"|1 0 0⟩"; u=1.0, t=1.0)) with dimension 2 and 2 stored entries:2×2 SparseArrays.SparseMatrixCSC{Float64, Int32} with 2 stored entries:
   ⋅   -1.0
 -1.0    ⋅

julia> using LinearAlgebra; round.(eigvals(Matrix(bsr)); digits = 4) # eigenvalues
3-element Vector{Float64}:
 -2.0
  1.0
  1.0

julia> ev = eigvecs(Matrix(bsr))[:,1]; ev = ev .* sign(ev[1]) # ground state eigenvector
3-element Vector{Float64}:
 0.5773502691896257
 0.5773502691896255
 0.5773502691896257

julia> dv = DVec(zip(bsr.basis, ev)) # ground state as DVec
DVec{BoseFS{1, 3, BitString{3, 1, UInt8}},Float64} with 3 entries, style = IsDeterministic{Float64}()
  fs"|1 0 0⟩" => 0.57735
  fs"|0 1 0⟩" => 0.57735
  fs"|0 0 1⟩" => 0.57735

Has methods for dimension, sparse, Matrix, starting_address.

Part of the AbstractHamiltonian interface. See also build_basis, AbstractOperator.

build_basis(
    ham, address=starting_address(ham);
    cutoff, filter, sizelim, sort=false, kwargs...
) -> basis
build_basis(ham, addresses::AbstractVector; kwargs...)

Get all basis element of a linear operator ham that are reachable (via non-zero matrix elements) from the address address, returned as a vector. Instead of a single address, a vector of addresses can be passed. Does not return the matrix, for that purpose use BasisSetRepresentation.

Providing an energy cutoff will skip addresses with diagonal elements greater than cutoff. Alternatively, an arbitrary filter function can be used instead. Addresses passed as arguments are not filtered.

Providing a max_depth will limit the size of the basis by only visiting addresses that are connected to the starting_address through max_depth hops through the Hamiltonian. Similarly, providing minimum_size will stop the bulding process after the basis reaches a length of at least minimum_size.

A maximum basis size sizelim can be set which will throw an error if the expected dimension of ham is larger than sizelim. This may be useful when memory may be a concern. These options are disabled by default.

The order the basis is returned in is arbitrary and non-deterministic. Use
`sort=true` if the ordering matters.

build_basis(addr::AbstractFockAddress)
build_basis(::Type{<:AbstractFockAddress}) -> basis

Return all possible Fock states of a given type as a vector. This method is much faster than build_basis(::AbstractHamiltonian, ...), but does not take matrix blocking into account. This version of build_basis accepts no additional arguments.

All address types except OccupationNumberFS are supported.

Returns a sorted vector of length equal to the dimension of addr.

See also AbstractFockAddress.

Matrix(
    hamiltonian::AbstractHamiltonian, addr=starting_address(hamiltonian);
    sizelim=10^4, kwargs...
)
Matrix(bsr::BasisSetRepresentation)

Return a dense matrix representation of hamiltonian or bsr. kwargs are passed to BasisSetRepresentation.

See BasisSetRepresentation.

sparse(hamiltonian::AbstractHamiltonian, addr=starting_address(hamiltonian); kwargs...)
sparse(bsr::BasisSetRepresentation)

Return a sparse matrix representation of hamiltonian or bsr. kwargs are passed to BasisSetRepresentation.

See BasisSetRepresentation.

LinearMap(::AbstractOperator{T}, basis)
LinearMap(::AbstractOperator{T}; basis)
LinearMap(
    op::AbstractOperator{T};
    starting_address=starting_address(op),
    full_basis=false
)

Wrapper for an AbstractOperator and a basis that allows multiplying regular Julia vectors with the operator without storing the matrix representation of the operator in memory.

If an op::AbstractOperator with no basis is passed, the basis is constructed automatically from the starting_address (note that function is only defined for AbstractHamiltonians). In that case, when full_basis=true the entire basis is constructed from an address as build_basis(starting_address), otherwise it is constructed as build_basis(op, starting_address). You may want to set full_basis=false when dealing with Hamiltonians that block, such as HubbardMom1D, otherwise setting full_basis=true is more efficient.

Implements the LinearMaps.jl interface, and can be used in Base.:*, mul! and the three-argument dot. These functions are multithreaded.

Example

julia> H = HubbardReal1D(BoseFS(1, 1, 1, 1));

julia> bsr = BasisSetRepresentation(H);

julia> v = ones(length(bsr.basis));

julia> w1 = bsr.sparse_matrix * v;

julia> op = LinearMap(H, bsr.basis);

julia> w2 = op * v;

julia> w1 ≈ w2
true

julia> dot(w1, bsr.sparse_matrix, v) ≈ dot(w1, op, v)
true

OperatorAsMap <: LinearMap

Wrapper over AbstractOperator that allows it to be used as a LinearMap from LinearMaps.jl.

See LinearMap for usage.

BasisSetRep(args...; kwargs...)

BasisSetRep is deprecated. Use BasisSetRepresentation instead.