Skip to content

Commit

Permalink
space_partition: Simplify algorithm of merge_subspaces()
Browse files Browse the repository at this point in the history 
The new non-recursive algorithm is simpler, faster and safe against
stack overflows.

An extra variant of the `merge_subspaces()` unit test has been added to
verify stability of the function when working with larger Hilbert
spaces.
  • Loading branch information
@krivenko
krivenko committed Aug 27, 2024
1 parent 70284ba commit ce898b2
Show file tree
Hide file tree
Showing 3 changed files with 113 additions and 92 deletions.
10 changes: 6 additions & 4 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,17 +4,19 @@ All notable changes to this project will be documented in this file.

## [0.8.0] - Unreleased

- Reduce the maximum allowed number of bits in the binary representation of
a basis state index (``hilbert_space::max_n_bits``) to 63. This way
``hilbert_space::dim()`` can return a valid value of type ``sv_index_type``
even when all 63 bits are used up.
- ``n_fermion_sector_view`` and ``n_fermion_multisector_view`` are now
parametrized on the type of ranking algorithm used to map basis state indices
from a full Hilbert space to a sector. Supported ranking algorithms are
``combination_ranking`` (selected by default), ``staggered_ranking`` and
``trie_ranking``. All three algorithms are described in M. Wallerberger,
[K. Held, Phys. Rev. Research 4, 033238 (2022)](
https://doi.org/10.1103/PhysRevResearch.4.033238).
- Reduced the maximum allowed number of bits in the binary representation of
a basis state index (``hilbert_space::max_n_bits``) to 63. This way
``hilbert_space::dim()`` can return a valid value of type ``sv_index_type``
even when all 63 bits are used up.
- Improved performance and stability of ``space_partition::merge_subspaces()``
by switching to a non-recursive variant of the algorithm.
- Fixed a negative index bug in ``n_fermion_sector_view``.
Credits to Dr. Cezary Śliwa for providing the patch.
- Whenever possible, use compiler intrinsics to speed up complex bit
Expand Down
68 changes: 19 additions & 49 deletions include/libcommute/loperator/space_partition.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@
#include "loperator.hpp"
#include "sparse_state_vector.hpp"

#include <functional>
#include <map>
#include <set>
#include <stdexcept>
Expand Down Expand Up @@ -169,56 +168,27 @@ class space_partition {
set_zeros(in_state);
});

// 'Zigzag' traversal algorithm
while(!Cd_conn.empty()) {

// Take one C^+ - connection
// C^+|lower_subspace> = |upper_subspace>

// C++17 structured bindings would be the real solution here
// NOLINTNEXTLINE(cppcoreguidelines-init-variables)
sv_index_type lower_subspace, upper_subspace;
std::tie(lower_subspace, upper_subspace) = *std::begin(Cd_conn);

// - Reveals all subspaces reachable from lower_subspace by application of
// a 'zigzag' product C^+ C C^+ C C^+ ... of any length.
// - Removes all visited connections from Cd_connections/C_connections.
// - Merges lower_subspace with all subspaces generated from
// lower_subspace by application of (C C^+)^(2*n).
// - Merges upper_subspace with all subspaces generated from
// upper_subspace by application of (C^+ C)^(2*n).
std::function<void(sv_index_type, bool)> zigzag_traversal =
[this,
lower_subspace,
upper_subspace,
&Cd_conn,
&C_conn,
&zigzag_traversal](
sv_index_type
in_subspace, // find all connections from in_subspace
bool upwards // if true, C^+ connection, otherwise C connection
) {
std::multimap<sv_index_type, sv_index_type>::iterator it;
while((it = (upwards ? Cd_conn : C_conn).find(in_subspace)) !=
(upwards ? Cd_conn : C_conn).end()) {

auto out_subspace = it->second;
(upwards ? Cd_conn : C_conn).erase(it);

if(upwards)
ds.set_union(out_subspace, upper_subspace);
else
ds.set_union(out_subspace, lower_subspace);

// Recursively apply to all found out_subspace's with
// a 'flipped' direction
zigzag_traversal(out_subspace, !upwards);
// Merge all 'out' subspaces corresponding to the same 'in' subspace
// in 'conn'.
auto merge_conn_targets =
[this](std::multimap<sv_index_type, sv_index_type> const& conn) {
if(conn.empty()) return;

auto conn_it = conn.cbegin();
sv_index_type in_subspace = conn_it->first;
sv_index_type out_subspace = conn_it->second;
++conn_it;
for(; conn_it != conn.cend(); ++conn_it) {
if(conn_it->first == in_subspace) {
ds.set_union(out_subspace, conn_it->second);
} else {
std::tie(in_subspace, out_subspace) = *conn_it;
}
};
}
};

// Apply to all C^+ connections starting from lower_subspace
zigzag_traversal(lower_subspace, true);
}
merge_conn_targets(Cd_conn);
merge_conn_targets(C_conn);

update_root_to_subspace();

Expand Down
127 changes: 88 additions & 39 deletions test/space_partition.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@

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

#include <algorithm>
#include <cmath>
Expand Down Expand Up @@ -343,54 +344,102 @@ TEST_CASE("Automatic Hilbert space partition", "[space_partition]") {

SECTION("merge_subspaces()") {

auto sp = space_partition(Hop, hs);
auto check_one_to_one = [](space_partition const& sp,
std::vector<decltype(Hop)> const& ops) {
// Calculated classification of states
std::vector<std::set<sv_index_type>> v_cl(sp.n_subspaces());
foreach(sp, [&](int i, int subspace) { v_cl[subspace].insert(i); });
std::set<std::set<sv_index_type>> cl{v_cl.cbegin(), v_cl.cend()};

sparse_state_vector<double> in_state(sp.dim());
sparse_state_vector<double> out_state(sp.dim());

for(auto const& op : ops) {
for(auto const& i_sp : cl) {
std::set<sv_index_type> f_sp;
for(auto i : i_sp) {
in_state[i] = 1.0;
op(in_state, out_state);
foreach(out_state, [&f_sp](sv_index_type f, double a) {
if(std::abs(a) < 1e-10) return;
f_sp.insert(f);
});
set_zeros(in_state);
}

// op maps i_sp onto zero
if(f_sp.size() == 0) continue;

// Check if op maps i_sp to only one subspace
auto n =
std::count_if(cl.begin(),
cl.end(),
[&f_sp](std::set<sv_index_type> const& f_sp_ref) {
return std::includes(f_sp_ref.cbegin(),
f_sp_ref.cend(),
f_sp.cbegin(),
f_sp.cend());
});
CHECK(n == 1);
}
}
};

std::vector<decltype(Hop)> ops1, ops2, all_ops;
for(int o = 0; o < n_orbs; ++o) {
ops1.emplace_back(c_dag("up", o) + c("dn", o), hs);
ops2.emplace_back(c("up", o) + c_dag("dn", o), hs);
SECTION("Hubbard-Kanamori") {
auto sp = space_partition(Hop, hs);

all_ops.emplace_back(ops1.back());
all_ops.emplace_back(ops2.back());
std::vector<decltype(Hop)> ops1, ops2, all_ops;
for(int o = 0; o < n_orbs; ++o) {
ops1.emplace_back(c_dag("up", o) + c("dn", o), hs);
ops2.emplace_back(c("up", o) + c_dag("dn", o), hs);

sp.merge_subspaces(ops1.back(), ops2.back(), hs);
}
all_ops.emplace_back(ops1.back());
all_ops.emplace_back(ops2.back());

// Calculated classification of states
std::vector<std::set<sv_index_type>> v_cl(sp.n_subspaces());
foreach(sp, [&](int i, int subspace) { v_cl[subspace].insert(i); });
std::set<std::set<sv_index_type>> cl{v_cl.cbegin(), v_cl.cend()};
sp.merge_subspaces(ops1.back(), ops2.back(), hs);
}

std::vector<double> in_state(sp.dim());
check_one_to_one(sp, all_ops);
}

SECTION("4 bath orbitals") {
int const n_bath_orbs = 4;
std::vector<double> const eps = {-0.2, -0.1, 0.1, 0.2};
double const V = 0.7;

for(auto const& op : all_ops) {
for(auto const& i_sp : cl) {
std::set<sv_index_type> f_sp;
for(auto i : i_sp) {
in_state[i] = 1.0;
auto out_state = op(in_state);
foreach(out_state, [&f_sp](sv_index_type f, double a) {
if(std::abs(a) < 1e-10) return;
f_sp.insert(f);
});
in_state[i] = 0;
for(std::string spin : {"dn", "up"}) {
for(int o = 0; o < n_bath_orbs; ++o) {
H += eps[o] * n(spin, n_orbs + o);
}
for(int o1 = 0; o1 < n_orbs; ++o1) {
for(int o2 = 0; o2 < n_bath_orbs; ++o2) {
H += V * (c_dag(spin, o1) * c(spin, n_orbs + o2) +
c_dag(spin, n_orbs + o2) * c(spin, o1));
}
}
}

// Hilbert space
auto hs_bath = make_hilbert_space(H);
// Linear operator form of the Hamiltonian
auto H_bath_op = make_loperator(H, hs_bath);

// op maps i_sp onto zero
if(f_sp.size() == 0) continue;

// Check if op maps i_sp to only one subspace
auto n =
std::count_if(cl.begin(),
cl.end(),
[&f_sp](std::set<sv_index_type> const& f_sp_ref) {
return std::includes(f_sp_ref.cbegin(),
f_sp_ref.cend(),
f_sp.cbegin(),
f_sp.cend());
});
CHECK(n == 1);
auto sp = space_partition(H_bath_op, hs_bath);

std::vector<decltype(Hop)> Cd, C, all_ops;
for(std::string spin : {"dn", "up"}) {
for(int o = 0; o < n_orbs + n_bath_orbs; ++o) {
Cd.emplace_back(c_dag(spin, o), hs_bath);
C.emplace_back(c(spin, o), hs_bath);

all_ops.emplace_back(Cd.back());
all_ops.emplace_back(C.back());

sp.merge_subspaces(Cd.back(), C.back(), hs_bath);
}
}

check_one_to_one(sp, all_ops);
}
}

Expand Down

0 comments on commit ce898b2

Please sign in to comment.