Sparse matrices

Working with many-body quantum systems often involves matrices which only have a small number of non-zero elements, also known as sparse matrices. Since storing and handling such objects in the conventional way (i.e. element by element) is very inefficient, XDiag includes basic implementations of three common sparse-matrix types: the coordinate (COO), the compressed-sparse-row (CSR), and the compressed-sparse-column (CSC) formats.

Just as the matrix function can be used to obtain the full matrix representing a given operator ops (in the form of an OpSum) on a given Hilbert space block (see the Dense matrices section), there are functions coo_matrix, csr_matrix, and csc_matrix to obtain the same matrix in the respective sparse format.

JuliaC++

coo_mat = coo_matrix(ops, block)
csr_mat = csr_matrix(ops, block)
csc_mat = csc_matrix(ops, block)

auto coo_mat = coo_matrix(ops, block);
auto csr_mat = csr_matrix(ops, block);
auto csc_mat = csc_matrix(ops, block);

Note that the C++ implementation distinguishes between real and complex matrices, e.g., there are the csr_matrix and csr_matrixC functions.

The objects returned by these functions are "raw" in the sense that they are not an instance of a sparse matrix implementation by another library, but contain all the information to call the respective sparse-matrix constructor of your sparse-matrix library of choice. For instance, the output of csc_matrix can be used to construct the sp_mat type implemented by the C++ Armadillo library or the SparseMatrixCSC type implemented by SparseArrays in julia.

JuliaC++

using SparseArrays
csc_mat = csc_matrix(ops, block)
jmat = SparseMatrixCSC(csc_mat.nrows, csc_mat.ncols, csc_mat.colptr, csc_mat.row, csc_mat.data)

auto csc_mat = csc_matrix(ops, block);
auto colptr = arma::conv_to<arma::uvec>::from(csc_mat.colptr);
auto row = arma::conv_to<arma::uvec>::from(csc_mat.row);
auto A = arma::sp_mat(colptr, row, csc_mat.data, csc_mat.nrows, csc_mat.ncols);

While the COO, CSR, and CSC formats are supported for extracting sparse matrices from XDiag, only the CSR format is used internally because it is the only one suitable for parallelized matrix-vector multiplications. The following julia example finds the ground state of an open-boundary Heisenberg chain and compares the default (matrix-free, i.e., on-the-fly) implementation to first converting the Hamiltonian to the CSR format.

Julia

# Heisenberg chain with open boundary conditions
N = 22
block = Spinhalf(N)
ops = OpSum()
for i in 1:(N-1)
    ops += "J" * Op("SdotS", [i, i+1])
end
ops["J"] = 1.0

# perform default (matrix-free) Lanczos iterations
@time begin
    println(eigval0(ops, block))
end

# perform Lanczos iterations on sparse CSR Hamiltonian
@time begin
    csr_mat = csr_matrix(ops, block)
    println(eigval0(csr_mat, block))
end

# Example output:
#
#   -9.568075875976074
#       6.859141 seconds (7 allocations: 512 bytes)
#   -9.568075875976076
#       5.954734 seconds (376 allocations: 2.000 GiB, 16.86% gc time)

It can be observed that (at medium block sizes) the runtime of both versions is comparable, while the matrix-free variant will always require substantially less memory. Note that the computational effort and memory required to construct the CSR representation grows significantly with system size such that the matrix-free version will generally run faster on larger systems. The CSR implementations of XDiag functions are thus only recommended when the CSR matrix can be pre-computed once and needs to be reused frequently.