# 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)