# Bit manipulation to find a monochromatic clique on k vertices knowing all on k-1 vertices

I have a very hot segment of code in a large project. I've extracted the relevant segment and a toy example to illustrate it.

The project involves an unavoidable combinatorial explosion related to graph theory, so there's no limit to how fast is 'good enough' -- but rather, each significant speed increase lets me explore slightly larger graphs. If there is a tradeoff between speed and clarity or safety or anything else but correctness, this needs speed.

I'm not only new to Rust, but am a mathematician rather than coder, so don't have expertise with any language. I'm vaguely aware that for code that needs to be fast, we might care about cache misses, CPU branch prediction, inlining, SIMD, and probably more things that I haven't even heard of.

What this code does is use u64 to represent sets. I.e., bit i is true if element i is present. I will have a possibly large Vec of these u64-sets, each of which has k members. Then, I'm looking for collections of k+1 of these that have k+1 set bits in their union.

This line: "if !(j_elements ^ candidate_elements).is_power_of_two()" is almost always true, and triggers a continue.

In addition to help with making this code faster, are there profiling tools people would recommend beyond flamegraph? E.g., some way to see if the compiled code or CPUs execution is non-optimal?

    use rustc_hash::FxHashSet;

fn main() {
let mut sets = vec![
TestSet { elements: 0b0011 },
TestSet { elements: 0b0101 },
TestSet { elements: 0b1001 },
TestSet { elements: 0b0110 },
];

let set_count = sets.len();
let subset_size = 3;
let mut explored_subsets = FxHashSet::<u64>::default();
let mut results = Vec::new();

for i in 0..=(set_count - subset_size) {
let i_elements = sets[i].elements;
for j in (i + 1)..=(set_count - subset_size + 1) {
let j_elements = sets[j].elements;
let candidate_elements = i_elements | j_elements;
if !(j_elements ^ candidate_elements).is_power_of_two() || explored_subsets.contains(&candidate_elements) {
continue;
}
explored_subsets.insert(candidate_elements);
let mut compatible_sets_needed = subset_size - 2;
for k in (j + 1)..set_count {
let k_elements = sets[k].elements;
if k_elements & candidate_elements == k_elements {
compatible_sets_needed -= 1;
if compatible_sets_needed == 0 {
results.push(candidate_elements);
break;
}
}
}
}
}

println!("Results: {:?}", results);
}

#[derive(Clone)]
struct TestSet {
elements: u64,
}


For context, while the toy code above only executes the nested loops once, in my full code they are executed many many times. As I said earlier, no speed is 'enough' but more is always better.

The full context is:

• I'm working with edge colorings of complete r-uniform hypergraphs on n vertices.
• An r-uniform hypergraph is an abstraction of a graph that redefines an edge to be a set of r vertices. I.e., an ordinary graph is a 2-uniform hypergraph.
• So let's take the simplest hypergraph that isn't an ordinary graph: The 3-uniform hypergraph on n edges. That is, r=3. Then there are choose(n, 3) edges. If we use merely 2 colors, there are 2^(choose(n, 3)) possible colorings (multi-counting isomorphisms)
• I'm doing simulated annealing with gradient-free optimization, in trying to find graphs where the count of maximal monochromatic cliques minimal. I.e., I randomly color the edges, use the provided code to find maximal cliques on k=4 vertices, then on k=5, then on k=6, and so on until there aren't enough monochromatic cliques for some k to enable k+1 to be possible.
• To find a monochromatic clique on k vertices, having already found all monochromatic cliques on k-1 vertices, I first look for a pair which have exactly k set bits in their union. Then, I need k-2 more cliques of the same color each of which has k-1 of these specific k bits set (and no other). I.e., the xor of the new clique and these k set bits is a single set bit, or a power of 2.
• The gradient descent step is flipping one edge and seeing if our count improved.

I'm not sure any of that context is important to my question, but I'm happy to provide it. What I'm asking is, assuming the nested for loops are called so often that these fast-and-simple operations are a bottleneck, how can I

1. tell which specific lines of code are taking the most time, and
2. if possible, how can I make them faster?

** update **

Here is the full method containing the above code (but without the variables renamed like they are above), and the struct it lives in:


struct HyperGraph {
edges: Vec<Clique>,
cliques: FxHashMap<u8, Vec<Clique>>,
members_of_cliques_which_should_be_deactivated: FxHashSet<u64>,
maximal_color_clique_ct: usize,
edge_order: u8,
color_ct: u8,
graph_order: u8,
graph_size: usize,
last_random_edge_index: usize,
set_bit_getter: math::SetBitGetter,
}

fn find_cliques_from_scratch(&mut self) {
self.cliques.clear();
self.maximal_color_clique_ct = 0;
let mut smaller_cliques_count: usize;
let mut order_usize: usize;
let mut compatible_cliques_needed: usize;
let mut key_small_cliques: u8;  // key for order-1 cliques
let mut explored_subsets = FxHashSet::<u64>::default();

for color in 0..self.color_ct {
// println!("\n  color: {}", color);

// Create all trivial color-cliques (those having order equal to EDGE_ORDER - 1)
Clique::generate_all_cliques(color, self.edge_order - 1, self.graph_order),
self.get_key(color, self.edge_order - 1));

// Create all color-edge-cliques
let mut edge_cliques = Vec::<Clique>::new();
for edge in &mut self.edges {
if edge.color == color {
edge_cliques.push(edge.clone());
}
}
edge_cliques,
self.get_key(color, self.edge_order));

// Mark fully contained trivial color-cliques as inactive (aka non-maximal)
self.mark_nonmaximal_cliques_inactive(color, self.edge_order);

// Create all color-order-cliques for orders > edge_order

// 'smaller_cliques' refers to all cliques that match the color
// we're examining, and have (order-1) vertices.
// The k'th bit of a cliques_index corresponds to the index in cliques[(color, order - 1)]

// These are used in a hot loop; allocating memory once here.
let mut candidate_clique_members: u64;
let mut cur_smaller_clique_members: u64;
let mut i_members: u64;
let mut j_members: u64;

for order in (self.edge_order + 1)..(self.graph_order + 1) {
// println!("    order: {}", order);

order_usize = order as usize;
key_small_cliques = self.get_key(color, order - 1);

// Confirm we have enough smaller cliques to build a bigger clique out of
let smaller_cliques = match self.cliques.get(&key_small_cliques) {
Some(cliques) => cliques,
None => break,
};

smaller_cliques_count = smaller_cliques.len();
if smaller_cliques_count < order_usize {
break;
}

explored_subsets.clear(); // Clear explored subsets

for i in 0..=(smaller_cliques_count - order_usize) {
i_members = smaller_cliques[i].members;
for j in (i + 1)..=(smaller_cliques_count - order_usize + 1) {
j_members = smaller_cliques[j].members;
candidate_clique_members = i_members | j_members;
if !(j_members ^ candidate_clique_members).is_power_of_two() || explored_subsets.contains(&candidate_clique_members) {
continue;
}
explored_subsets.insert(candidate_clique_members);
compatible_cliques_needed = order_usize - 2;
for k in (j + 1)..smaller_cliques_count {
cur_smaller_clique_members = smaller_cliques[k].members;
if cur_smaller_clique_members & candidate_clique_members == cur_smaller_clique_members {
compatible_cliques_needed -= 1;
if compatible_cliques_needed == 0 {
// We've found an order-clique!
}
continue;
}
}
}
}

while let Some(clique) = cliques_to_add.pop() {
}

// Having finished finding color-order cliques, mark any nonmaximal smaller cliques as inactive
self.mark_nonmaximal_cliques_inactive(color, order);
}
}
}
$$$$

• Also, is there a way to get flamegraph to tell me which line number it's spending time on? I know almost all of my time is being spent somewhere in this nested loop, but don't know which line(s) are particularly expensive.
– Dave
Commented Aug 11 at 1:48
• So do the sets all fit in u64, or are you using multiple u64 per set? I'm not familiar with Rust. Are you sure there isn't a faster approach with dynamic programming?
– qwr
Commented Aug 11 at 2:30
• @qwr It's a single u64 per set, and I could get away with a u32, though not u16. In my tests, u64 ran as fast as, possibly slightly faster than, u32 so I stuck with u64.
– Dave
Commented Aug 11 at 2:32
• 64-bit is almost certainly your processor's native word size, so it makes sense it would be fastest. Although 32-bit might fit twice as much into caches.
– qwr
Commented Aug 11 at 2:38
• Also flipping one edge and seeing if count improved is not really gradient descent, unless you really have a meaningful gradient from finite differences. I would call it gradient-free optimization. Genetic algorithms or particle swarms would be an interesting approach.
– qwr
Commented Aug 11 at 4:35

## Initial Reasoning

Let V be the order of the graph. Let Q = subset_size. Let a Q-set be a u64 with Q set bits. We can get a speedup of O(choose(V, Q)^2 / (Q * V)) times by improving the algorithm.

Start with the innermost loop, where the iterator k traverses (j+1)..set_count. For worst-case input, this means O(set_count) = O(choose(V, Q)) iterations on average for each i, j.

Instead of checking all remaining (Q-1)-sets for Q-2 other subsets of the candidate Q-set, construct all Q subsets of the candidate, which are (Q-1)-sets, and check their inclusion.

You can do so by initializing a FxHashSet::<u64> included from sets such that included[x] records whether x is a member of sets. Once you've computed candidate_elements, you can take out one element at a time and test the resulting subset with included.contains, which runs in O(1) amortized time. Since there are Q elements to test in candidate_elements, this modification replaces the O(choose(V, Q)) innermost loop with an O(Q) loop.

Now let's move to the middle loop with iterator j. As you observed, the continue statement runs very frequently. This suggests that reformulating the loop could improve performance.

Instead of constructing candidate Q-sets by the union of a pair of (Q-1)-sets, take each (Q-1)-set and try adding each remaining vertex in isolation.

This replaces the O(choose(V, Q)) iterations of the middle loop with O(V) iterations. You can go even faster by observing that you only need to try remaining vertices that are larger than the vertices in the clique. In other words, aim to find each Q-set from the (Q-1)-set produced by removing the vertex corresponding to the most significant bit of the Q-set. Notably, this removes the need for explored_subsets because each candidate_elements will be unique. This approach has optimal O(set_count) space complexity:

use std::collections::HashSet;

// Is setdiff(cand, {x}) in 'included' for all x in cand?
fn can_expand(included: &HashSet::<u64>, cand: u64) -> bool
{
let mut unit = 1; // single-element set {x}

while unit < cand {
if unit & cand == unit {
let setdiff = cand - unit;
if !included.contains(&setdiff) {
return false;
}
}

unit <<= 1; // from {x} to {x+1}
}

return true;
}

// Get all expansions of 'sets'.
fn get_expansions(sets: &Vec::<u64>) -> Vec<u64>
{
// handle sets = []
if sets.len() == 0 {
return vec![];
}

// handle sets = [0]
if sets[0] == 0 {
return vec![1];
}

let max_set = *(sets.iter().max().unwrap());
let included = sets.iter().cloned().collect();
let mut results = Vec::new();

for &set in sets {
let mut unit = 1; // to be added

// only add elements larger than max(set)
while unit <= set {
unit <<= 1;
}

// add the other elements one-by-one
while unit <= max_set {
let candidate = set + unit;

if can_expand(&included, candidate) {
results.push(candidate);
}

unit <<= 1;
}
}

return results;
}


## Major Update

Here's another idea that uses O(set_count * V) space. You can get a speedup of O(choose(V, Q)^2 / V) times. This approach has optimal O(set_count * V) time complexity:

use std::collections::HashMap;

// Gets the number of elements in 'set'.
fn get_set_bits(set: u64) -> u64
{
let mut set_bits = 0;
let mut unit = 1;

while unit <= set {
if unit & set == unit {
set_bits += 1;
}

unit <<= 1;
}

return set_bits;
}

// Get all expansions of 'sets'.
fn get_expansions(sets: &Vec::<u64>) -> Vec<u64>
{
// handle sets = []
if sets.len() == 0 {
return vec![];
}

// handle sets = [0]
if sets[0] == 0 {
return vec![1];
}

// let 2^(V-1) <= max_set, Q = min_bits + 1
let max_set = *(sets.iter().max().unwrap());
let min_bits = get_set_bits(sets[0]);

// votes[cand] is the tally of sets for cand

// each (Q-1)-set votes for candidate Q-sets
for &set in sets {
let mut unit = 1;

while unit <= max_set {
if unit & set == 0 {
let cand = set + unit;
*tally.or_insert(0) += 1;
}

unit <<= 1;
}
}

// return the candidates with Q tallied votes
let mut results = Vec::new();

for (cand, tally) in votes {
if tally > min_bits {
results.push(cand);
}
}

return results;
}


Imagine a population sets, each a (Q-1)-set, voting for Q-set candidates in results. As an example, 0b1010 would vote for 0b1110 and 0b1011. We use a HashMap to tally the votes for each candidate. Then simply return the candidates with Q or more votes. This has more similarity to your initial approach, but with one key distinction: instead of using explored_subsets to avoid duplicates, we use duplicates to count votes for that candidate.

## Test Cases

fn test_expand(sets: &Vec::<u64>, want: &Vec::<u64>) {
assert!(&get_expansions(sets) == want);
}

fn main() {
test_expand(
&vec![],
&vec![]
);

test_expand(
&vec![0b0],
&vec![0b1],
);

test_expand(
&vec![0b1],
&vec![],
);

test_expand(
&vec![0b10, 0b01],
&vec![0b11],
);

test_expand(
&vec![0b100000, 0b000100],
&vec![0b100100],
);

test_expand(
&vec![0b0011, 0b0101, 0b1001, 0b0110],
&vec![0b0111],
);
test_expand(
&vec![0b01011, 0b10011, 0b11001, 0b11010],
&vec![0b11011],
);
}


Hope this helps! Your demonstration was well-contained; I was able to get it running despite having never tried Rust before. I used standard libraries due to installation issues.

• Great stuff, thanks!
– Dave
Commented Aug 19 at 18:33
• If you like optimization, you might have fun with the algorithms section on S/O (stackoverflow.com/questions/tagged/algorithm)
– Dave
Commented Aug 19 at 22:25
• Thanks for the link! Would love to hear if the changes I proposed translate into runtime improvements for your research problem. It sounds very interesting. Commented Aug 20 at 6:43
• The speedup depends on the inputs, but call it roughly 5x, which is helpful. For my problem, not producing 0b11 as described doesn't come up, but I may reuse this in another context someday where it matters.
– Dave
Commented Aug 20 at 10:27
• I'm on vacation; will check it out in a wk. Thanks!
– Dave
Commented Aug 24 at 1:30

insert() returns whether the set has grown.
            if !(j_elements ^ candidate_elements).is_power_of_two()
`