Spin-1 Heisenberg ring: Compressed storage of state vectors

In this example we explore the concept of a sparse Hilbert space, and illustrate how using the compressed_state_view objects can lead to a dramatic reduction in memory footprint. To this end, we consider a ring of \(N\) spins \(S=1\) interacting via the antiferromagnetic Heisenberg exchange coupling,

\[\begin{split}\hat H = \sum_{i=0}^{N-1} \mathbf{S}_i \cdot \mathbf{S}_{(i+1)\bmod N} = \sum_{\substack{i=0,\\j = (i+1)\bmod N}}^{N-1} S_{z,i} S_{z,j} + \frac{1}{2} \left[ S_{+,i} S_{-,j} + S_{-,i} S_{+,j} \right].\end{split}\]

The state space of each spin is spanned by the \(S_z\)-triplet \(|0\rangle, |\pm1\rangle\), and the dimension of the model’s Hilbert space is thus \(3^N\). From the computing standpoint, however, it takes 2 bits to specify one of the three basis states, with the bit combination 0b11 not representing any physical state. According to the convention adopted in libcommute, one then needs a bit string of length \(2N\) to encode a basis state of the entire system. Although not all such strings correspond to physical basis states (most of the strings are, in fact, unphysical for \(N\geq 3\)), naively one would have to allocate an array of size \(4^N\) to store a state vector compatible with linear operators acting in this Hilbert space. This situation is an example of a sparse Hilbert space, i.e. of a space whose basis state indices do not form a continuous integer range. Sparse Hilbert spaces are characterized by non-power-of-two dimensions.

It is a valid (and faster) option to allocate std::vector (Eigen::Vector, etc) objects of the size \(4^N\) to store state vectors, even though some of their elements will never be referenced. Another option offered by libcommute is to allocate containers of the exponentially smaller size \(3^N\), and to employ the compressed_state_view adaptor that – at the cost of a certain speed loss – translates indices of the physical basis states onto the continuous range \([0;3^N-1]\).

// 'state_vector_eigen3.hpp' must be included before 'libcommute.hpp'.
// This way libcommute will know how to deal with Eigen's vector types.
#include <libcommute/loperator/state_vector_eigen3.hpp>

#include <libcommute/libcommute.hpp>

#include <Eigen/Eigenvalues>

using namespace libcommute;

int main() {
  int const N = 6; // Number of spins

  // Expression of the Heisenberg Hamiltonian to be constructed.
  // Every creation and annihilation operator found in the expression carries
  // one integer index (ring site).
  static_indices::expr_real<int> H;

  // The following functions return operators with statically typed indices.
  using static_indices::S_z; // Spin projection operator S_z
  using static_indices::S_p; // Spin-raising operator S_+
  using static_indices::S_m; // Spin-lowering operator S_-

  // Heisenberg interaction of the S=1 spins
  for(int i = 0; i < N; ++i) {
    int j = (i + 1) % N;
    H += S_z<3>(i) * S_z<3>(j) +
         0.5 * (S_p<3>(i) * S_m<3>(j) + S_m<3>(i) * S_p<3>(j));
  }

  // Automatically analyse structure of 'H' and construct a Hilbert space.
  auto hs = make_hilbert_space(H);

  // The Hilbert space 'hs' is sparse: The minimal required container size for
  // state vectors exceeds the dimension of the space.
  std::cout << "Is the Hilbert space sparse? " << hs.is_sparse() << '\n';
  std::cout << "Hilbert space dimension: " << hs.dim() << '\n';
  std::cout << "Minimal container size for state vectors: " << hs.vec_size()
            << '\n';

  auto const dim = hs.dim(); // 3^N

  // Construct an 'loperator' object that represents action of the expression
  // 'H' on state vectors in the Hilbert space 'hs'.
  auto Hop = make_loperator(H, hs);

  // Preallocate an input state vector of the length 3^N.
  // This container is compressed, i.e. it stores quantum amplitudes of
  // the *physical* basis states only!
  Eigen::VectorXd psi(dim);

  // Preallocate a Hamiltonian matrix to be filled and diagonalized.
  // This matrix is also compressed.
  Eigen::MatrixXd Hmat(dim, dim);

  // Create a view of the compressed state 'psi'.
  auto psi_view = make_const_comp_state_view(psi, hs);

  // Fill 'Hmat'
  // Important: k takes on values corresponding to the *physical* basis states
  // only.
  foreach(hs, [dim, &hs, &Hop, &psi, &psi_view, &Hmat](sv_index_type k) {
    // Translate k into the compressed (continuous) range [0; 3^N-1]
    auto const k_comp = psi_view.map_index(k);

    // Reset vector |\psi> to zero.
    psi = Eigen::VectorXd::Zero((int)dim);

    // Now, 'psi' is the basis vector corresponding to 'k'.
    psi((int)k_comp) = 1;

    // Create a view of the column of 'Hmat' that corresponds to 'k'.
    auto HMat_col_view = make_comp_state_view(Hmat.col((int)k_comp), hs);

    // Act with 'Hop' on the basis state 'k' and write the result into
    // the column of 'Hmat' that also corresponds to 'k'.
    Hop(psi_view, HMat_col_view);
  });

  // Diagonalize 'Hmat'.
  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> diag(Hmat,
                                                      Eigen::EigenvaluesOnly);

  std::cout << "10 lowest eigenvalues of H = "
            << diag.eigenvalues().head(10).transpose() << '\n';

  return 0;
}