Finding invariant subspaces of a linear operator

This page explains how to use space_partition, an efficient tool to performs the following two tasks.

Partition of a Hilbert space into a set of disjoint subspaces invariant under action of a given Hermitian operator \(\hat H\) (Hamiltonian).
Merge some of the invariant subspaces together to ensure that a given operator \(\hat O\) and its Hermitian conjugate \(\hat O^\dagger\) generate only one-to-one connections between the subspaces.

Partition of the Hilbert space is a preparatory step that can immensely reduce complexity of an exact diagonalization problem by reducing it to a series of smaller problem.

Merging some of the invariant subspaces increases computational costs of diagonalization but simplifies computation of many-body correlation functions of operators \(\hat O\), \(\hat O^\dagger\).

For a detailed description of the algorithm see [SKFP16].

template<typename ScalarType> using matrix_elements_map = std::map<std::pair<sv_index_type, sv_index_type>, ScalarType>;: Sparse storage type for matrix elements of a linear operator \(\hat O\). Pairs of indices \((i,j)\) are mapped to \(\langle i| \hat O |j\rangle\).

class space_partition

Partition of a Hilbert space into disjoint subspaces.

template<typename HSType, typename LOpScalarType, int... LOpAlgebraIDs> space_partition(loperator<LOpScalarType, LOpAlgebraIDs...> const &h, HSType const &hs): Partition Hilbert space hs into invariant subspaces of Hermitian linear operator h (Hamiltonian).

template<typename HSType, typename LOpScalarType, int... LOpAlgebraIDs> space_partition(loperator<LOpScalarType, LOpAlgebraIDs...> const &h, HSType const &hs, loperator_melem_t<LOpScalarType, LOpAlgebraIDs...> &me): Partition Hilbert space hs into invariant subspaces of Hermitian linear operator h (Hamiltonian). This constructor reveals all matrix elements of h and writes them into me.

template<typename HSType, typename LOpScalarType, int... LOpAlgebraIDs> auto merge_subspaces(loperator<LOpScalarType, LOpAlgebraIDs...> const &Od, loperator<LOpScalarType, LOpAlgebraIDs...> const &O, HSType const &hs, bool store_matrix_elements = true) -> std::pair<loperator_melem_t<LOpScalarType, LOpAlgebraIDs...>, loperator_melem_t<LOpScalarType, LOpAlgebraIDs...>>

Merge some of the found invariant subspaces to ensure that linear operators Od and O generate only one-to-one connections between the subspaces.

Matrix elements of Od and O will be returned as a pair of matrix_elements_map objects if store_matrix_elements == true. Otherwise, the returned maps will be empty.

This method can be called multiple types to make the partition simultaneously fulfil the one-to-one conditions for many \(\hat O_i\), \(\hat O^\dagger_i\) pairs.

sv_index_type dim() const: Dimension of the Hilbert space used to construct this partition.

sv_index_type n_subspaces() const

Number of disjoint subspaces in this partition.

Calls to merge_subspaces() can decrease this number.

sv_index_type operator[](sv_index_type index) const

Return serial number of the disjoint subspace the basis state with index belongs to.

The disjoint subspaces can be re-enumerated after a call to merge_subspaces().

template<typename HSType, typename LOpScalarType, int... LOpAlgebraIDs> auto find_connections(loperator<LOpScalarType, LOpAlgebraIDs...> const &op, HSType const &hs) -> std::set<std::pair<sv_index_type, sv_index_type>>: Analyze connections between subspaces generated by operator op. The connections are returned as a set of pairs of subspace serial numbers, (source subspace, destination subspace).

std::vector<sv_index_type> subspace_basis(sv_index_type index) const: Build and return a list of indices of all basis states spanning a given subspace index.

std::vector<std::vector<sv_index_type>> subspace_bases() const: Build and return lists of indices of all basis states spanning all subspaces in the partition. The returned lists are disjoint and their union spans the entire Hilbert space.

template<typename F> friend void foreach(space_partition const &sp, F &&f): Apply functor f to all basis states in a given space partition. The functor must take two arguments, index of the basis state, and serial number of the subspace this basis state belongs to.

Space partition example

#include <libcommute/expression/factories.hpp>
#include <libcommute/loperator/mapped_basis_view.hpp>
#include <libcommute/loperator/space_partition.hpp>

#include <iostream>
#include <vector>

using namespace libcommute;

int main() {

  using namespace static_indices; // For c(), c_dag() and n()

  //
  // Build Hamiltonian of the 3-orbital Hubbard-Kanamori atom
  //

  int const n_orbs = 3;
  double const mu = 0.7; // Chemical
  double const U = 3.0;  // Coulomb repulsion
  double const J = 0.3;  // Hund interaction

  // Hamiltonian
  expression<double, std::string, int> H;

  // Chemical potential terms
  for(int o = 0; o < n_orbs; ++o) {
    H += -mu * (n("up", o) + n("dn", o));
  }

  // Intraorbital interactions
  for(int o = 0; o < n_orbs; ++o) {
    H += U * n("up", o) * n("dn", o);
  }
  // Interorbital interactions, different spins
  for(int o1 = 0; o1 < n_orbs; ++o1) {
    for(int o2 = 0; o2 < n_orbs; ++o2) {
      if(o1 == o2) continue;
      H += (U - 2 * J) * n("up", o1) * n("dn", o2);
    }
  }
  // Interorbital interactions, equal spins
  for(int o1 = 0; o1 < n_orbs; ++o1) {
    for(int o2 = 0; o2 < n_orbs; ++o2) {
      if(o2 >= o1) continue;
      H += (U - 3 * J) * n("up", o1) * n("up", o2);
      H += (U - 3 * J) * n("dn", o1) * n("dn", o2);
    }
  }
  // Spin flip and pair hopping terms
  for(int o1 = 0; o1 < n_orbs; ++o1) {
    for(int o2 = 0; o2 < n_orbs; ++o2) {
      if(o1 == o2) continue;
      H += -J * c_dag("up", o1) * c_dag("dn", o1) * c("up", o2) * c("dn", o2);
      H += -J * c_dag("up", o1) * c_dag("dn", o2) * c("up", o2) * c("dn", o1);
    }
  }

  //
  // Hilbert space of the problem
  //

  auto hs = make_hilbert_space(H);

  //
  // Linear operator form of the Hamiltonian
  //

  auto Hop = make_loperator(H, hs);

  //
  // Reveal invariant subspaces of H
  //

  auto sp1 = space_partition(Hop, hs);

  std::cout << "Total dimension of the Hilbert space is " << sp1.dim() << '\n';

  std::cout << "H has " << sp1.n_subspaces() << " invariant subspaces\n";

  //
  // Once again, this time saving the matrix elements of H
  //

  matrix_elements_map<double> H_elements;
  auto sp2 = space_partition(Hop, hs, H_elements);

  std::cout << "H has " << H_elements.size() << " non-vanishing matrix elements"
            << '\n';

  //
  // Now merge some invariant subspaces to make sure that all electron
  // creation and annihilation operators connect one subspace to one subspace.
  //

  for(std::string spin : {"up", "dn"}) {
    for(int o = 0; o < n_orbs; ++o) {
      auto Cdagop = make_loperator(c_dag(spin, o), hs);
      auto Cop = make_loperator(c(spin, o), hs);

      auto matrix_elements = sp2.merge_subspaces(Cdagop, Cop, hs);

      std::cout << c_dag(spin, o) << " has " << matrix_elements.first.size()
                << " non-vanishing matrix elements\n";
      std::cout << c(spin, o) << " has " << matrix_elements.first.size()
                << " non-vanishing matrix elements\n";
    }
  }

  std::cout << "Number of disjoint subspaces after merging is "
            << sp2.n_subspaces() << '\n';

  //
  // foreach()
  //

  std::cout << "Distribution of basis states over the subspaces:\n";
  foreach(sp2, [](sv_index_type index, sv_index_type subspace) {
    std::cout << index << " => " << subspace << '\n';
  });

  //
  // Bonus: Using mapped_basis_view to act on a vector in one of the invariant
  // subspaces.
  //

  // Collect all basis state indices from one subspace
  std::vector<sv_index_type> basis_states_in_subspace_24 =
      sp2.subspace_basis(24);

  auto sp24_dim = basis_states_in_subspace_24.size();
  std::cout << "Subspace 24 is " << sp24_dim << "-dimensional\n";

  // Make a mapper object
  basis_mapper mapper(basis_states_in_subspace_24);

  // 'psi' and 'phi' are small vectors entirely lying in the 24-th subspace
  std::vector<double> psi(sp24_dim), phi(sp24_dim);

  //
  // These views have dimension of the full Hilbert space!
  //
  auto psi_view = mapper.make_const_view(psi);
  auto phi_view = mapper.make_view(phi);

  //
  // Act on 'psi_view' as if it was a vector in the full space
  // No exception here because 'Hop' connects subspace 24 only to itself.
  //
  Hop(psi_view, phi_view);

  return 0;
}

[SKFP16]

“TRIQS/CTHYB: A continuous-time quantum Monte Carlo hybridisation expansion solver for quantum impurity problems”, P. Seth, I. Krivenko, M. Ferrero and O. Parcollet, Comp. Phys. Comm. 200, March 2016, 274-284, http://dx.doi.org/10.1016/j.cpc.2015.10.023 (section 4.2)