A 90 min recording of a talk I gave on this post can be found here.

1 Introduction

In this post, we will develop a fast implementation of random minimizers.

The goals here are:

• a library for fast random minimizers,
• explaining the final code, and the design decisions leading to it,
• walking you through the development process, hopefully teaching some tricks and methods for how to do this yourself.

Code for this post is in the benches/ directory of github:RagnarGrootKoerkamp/minimizers. Most code snippets are directly included from the repository. I have tried to write clean and readable code that follows the order of this post. Occasionally some code snippets contain ‘future proofing’ that will only become relevant later on.

Contents¹:

• Background
• Overview of algorithms
• Analysing and optimizing the algorithms
• Rolling our own hash function
• SIMD sliding window

1.1 Results

The end result that we will be working towards here is twofold:

1. An implementation of 32-bit ntHash rolling hash that takes around 0.4ns/kmer, or around 1 CPU cycle per hash. This is 7x faster than the nthash crate by Luiz Irber.

2. An implementation of random minimizers and super-k-mers that takes around 1.6 ns/window, or 4 CPU cycles per window. This is 5x faster than Daniel Liu’s gist and 15x faster than the rust-seq/minimizer-iter crate by Igor Martayan².

   It computes all window minimizers of a 2-bit packed human genome (with non-ACGT stripped) in 5.4 seconds.

3. An implementation of canonical minimizers that runs in 1.9 ns/window, or 6.4 seconds for a human genome.

Some of the code to iterate bases is already available in the packed-seq crate in the rust-seq org. The remaining algorithms are planned to be upstreamed to rust-seq/nthash and/or rust-seq/minimizer-iter crates, after polishing the implementation, adding tests³, and designing a nice API⁴.

I welcome any (critical) feedback on both the content and styling of this post! I plan to fine-tune it a bit more based on this, and will take this into account when writing similar posts in the future. If you’re feeling adventurous, open an issue or pull request against the source of this blog.

2 Random minimizers

Many tools in bioinformatics rely on k-mers. Given a string/text/DNA sequence t: &[u8] of length \(|t|=n\), the k-mers of t are all the length-\(k\) substrings. There are \(n-k+1\) of them.

In order to speed up processing, and since consecutive k-mers are overlapping and contain redundancy, it is common to subsample the k-mers. This should be done in a deterministic way, so that similar but slightly distinct text still mostly sample the same k-mers. One such way is random minimizers, introduced in parallel by Schleimer, Wilkerson, and Aiken (2003) and Roberts et al. (2004):

1. First, a window size \(w\) is chosen, which will guarantee that at least one k-mer is sampled every \(w\) positions of \(t\), to prevent long gaps between consecutive k-mers.
2. Then, every k-mer is hashed using a pseudo-random hash function \(h\), and in each window of \(w\) consecutive k-mers (of total length \(\ell=w+k-1\)), the leftmost⁵ one with the smallest hash is selected.

Consecutive windows (that overlap in \(w-1\) k-mers) often sample the same k-mer, so that the total number of sampled k-mers is much smaller than \(n-k+1\). In fact, random minimizers have a density of \(2/(w+1)\) and thus sample approximately two k-mers per window.

My previous post on minimizers gives more background on density and other sampling schemes. One such alternative scheme is the mod-minimizer (Groot Koerkamp and Pibiri 2024), that has low density when \(k\) is large compared to \(w\). That paper also contains a review of methods and pseudocode for them. In this post, I will focus on the most widely used random minimizers.

3 Algorithms

3.1 Problem statement

Before going into algorithms to compute random minimizers, let’s first precisely state the problem itself. In fact, we’ll state three versions of the problem, that each include slightly more additional information in the output.

Problem A: Only the set of minimizers

The most basic form of computing minimizers is to simply return the list of minimizer kmers and their positions. This information is sufficient to sketch the text.

1
2
3
4

pub trait Minimizer {
    /// Problem A: The absolute positions of all minimizers in the text.
    fn minimizer_positions(&mut self, text: &[u8]) -> Vec<usize>;
}

Problem B: The minimizer of each window

When building datastructures on top of a sequence, such as SSHash (Pibiri 2022), just the set of k-mers is not sufficient. Each window should be mapped to its minimizer, but there could be multiple minimizers in a single window. Thus, in this version of the problem, we return the minimizer of each window.

Note that while in Problem A we return one result per minimizer, here we return one result per window (roughly, per character). This difference can have implications for the optimal algorithm implementation, since it’s not unlikely that the main loop of the optimal solution has exactly one iteration per returned element.

1
2
3
4
5
6

pub trait Minimizer {
    ...

    /// Problem B: For each window, the absolute position in the text of its minimizer.
    fn window_minimizers(&mut self, text: &[u8]) -> Vec<usize>;
}

Problem C: Super-k-mers

We can implement SSHash and other indexing datastructures using the output of Problem B, but in most such applications, it is probably more natural to directly iterate over the super-k-mers: the longest substrings such that all contained windows have the same minimizer. This can easily be obtained from the output of Problem B by grouping consecutive windows that share the minimal k-mer.

Note that here the result is again an iterator over minimizers, like in Problem A rather than windows as in Problem B.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

struct SuperKmer {
    /// The absolute position in the text where this super-k-mer starts.
    /// The end of the super-k-mer can be inferred from the start of the next super-k-mer.
    start_pos: usize,
    /// The absolute position in the text of the minimizer of this super-k-mer.
    minimizer_pos: usize,
}

pub trait Minimizer {
    ...

    /// Problem C: The super-k-mers of the text.
    fn super_kmers(&mut self, text: &[u8]) -> Vec<SuperKmer>;
}

Which problem to solve

In the remainder, we mostly focus on Problem B, since that is what I initially started with. However, some methods more naturally lean themselves to solve the other variants. We will benchmark each method on whichever of the problems it runs fastest, and postpone optimizing a solution for each specific problem until later.

Most likely though, solving Problem C will yield the most useful tool for end users. I plan to revisit this at the end of this post.

Regardless, for convenience we can implement solutions to A and C in terms of B.⁶

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

pub trait Minimizer {
    /// Problem A: The absolute positions of all minimizers in the text.
    fn minimizer_positions(&mut self, text: &[u8]) -> Vec<usize> {
        let mut minimizers = self.window_minimizers(text);
        minimizers.dedup();
        minimizers
    }

    /// Problem B: For each window, the absolute position in the text of its minimizer.
    fn window_minimizers(&mut self, _text: &[u8]) -> Vec<usize> {
        unimplemented!()
    }

    /// Problem C: The super-k-mers of the text.
    fn super_kmers(&mut self, text: &[u8]) -> Vec<SuperKmer> {
        self.window_minimizers(text)
            .into_iter()
            .enumerate()
            .group_by(|(_idx, minimizer_pos)| *minimizer_pos)
            .into_iter()
            .map(|(minimizer_pos, mut group)| SuperKmer {
                start_pos: group.next().expect("groups are non-empty").0,
                minimizer_pos,
            })
            .collect()
    }
}

Canonical k-mers

In actual bioinformatics applications, it is common to consider canonical k-mers and corresponding canonical minimizers. These are defined such that the canonical k-mer of a string and its reverse complement are the same. This will significantly impact the complexity of our code, and hence we also postpone this for later.

3.2 The naive algorithm

The naive \(O(n \cdot w)\) method simply finds the minimum of each window independently.

It looks like this:⁷

In code, it looks like this.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

pub struct NaiveMinimizer<H = FxHash> {
    pub w: usize,
    pub k: usize,
    pub hasher: H,
}

impl<H: Hasher> Minimizer for NaiveMinimizer<H>
where
    H::Out: Ord,
{
    fn window_minimizers(&mut self, text: &[u8]) -> Vec<usize> {
        // Iterate over the windows of size l=w+k-1.
        text.windows(self.w + self.k - 1)
            .enumerate()
            // For each window, starting at pos j, find the lexicographically smallest k-mer.
            .map(|(j, window)| {
                j + window
                    .windows(self.k)
                    .enumerate()
                    .min_by_key(|(_idx, kmer)| self.hasher.hash(kmer))
                    .expect("w > 0")
                    .0
            })
            .collect()
    }
}

Note that we assume that hashing each k-mer takes constant time. This is reasonable in practice, since most hash functions hash many characters at once.

Performance characteristics

I was planning to show some preliminary performance numbers here, to give some intuition about the relative performance of this and upcoming methods. Unfortunately though (and not completely surprising), performance depends heavily on small details of how the code is written. I could provide a comparison between the methods as-is, but the numbers would be hard to interpret. Instead, let’s consider some high level expectations of this algorithm for now, and postpone the quantitive comparison until later.

Since there are no if statements, the code should be highly predictable. Thus, we expect:

• few branch misses,
• high instructions per cycle,
• and generally running times linear in \(w\).

3.3 Rephrasing as sliding window minimum

After hashing all k-mers, we basically have a sequence of \(n-k+1\) pseudo-random integers. We would like to find the position of the leftmost minimum in each window of \(w\) of those integers.

Thus, the problem of finding random minimizers can be reformulated in terms of the sliding window minima. We can model this problem using the following trait:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

/// A value at an absolute position.
/// When comparing, ties between value are broken in favour of small position.
#[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Debug)]
pub struct Elem<V> {
    pub val: V,
    pub pos: usize,
}

pub trait SlidingMin<V> {
    /// Take an iterator over values of type V.
    /// Return an iterator over the minima of windows of size w, and their positions.
    fn sliding_min(&self, w: usize, it: impl Iterator<Item = V>) -> impl Iterator<Item = Elem<V>>;
}

A first implementation that mirrors the Naive implementation above could simply store the last \(w\) hashes in a RingBuffer, and in each step recompute the minimum. The boiler plate for the ring buffer is simple:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

pub struct RingBuf<V> {
    w: usize,
    idx: usize,
    data: Vec<V>,
}

impl<V: Clone> RingBuf<V> {
    #[inline(always)]
    pub fn new(w: usize, v: V) -> Self {
        let data = vec![v; w];
        RingBuf { w, idx: 0, data }
    }

    #[inline(always)]
    pub fn idx(&self) -> usize {
        self.idx
    }

    #[inline(always)]
    pub fn push(&mut self, v: V) {
        self.data[self.idx] = v;
        self.idx += 1;
        if self.idx == self.w {
            self.idx = 0;
        }
    }
}

/// A RingBuf can be used as a slice.
impl<V> std::ops::Deref for RingBuf<V> {
    type Target = [V];

    #[inline(always)]
    fn deref(&self) -> &[V] {
        &self.data
    }
}

And then a sliding window implementation also becomes trivial:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

pub struct Buffered;

impl<V: Copy + Max + Ord> SlidingMin<V> for Buffered {
    fn sliding_min(&self, w: usize, it: impl Iterator<Item = V>) -> impl Iterator<Item = Elem<V>> {
        // A ring buffer that holds the w last elements.
        let mut ring_buf = RingBuf::new(w, Elem { val: V::MAX, pos: 0 });
        // Iterate over items and their positions.
        it.enumerate()
            .map(move |(pos, val)| {
                ring_buf.push(Elem { val, pos });
                // Return the minimum element in the buffer.
                // Ties between values are broken by the position.
                *ring_buf.iter().min().expect("w > 0")
            })
            // The first w-1 elements do not complete a window.
            .skip(w - 1)
    }
}

Using this trait, we can clean up⁸ the implementation of the Naive minimizer from before⁹. In fact, we can now generically implement a minimizer scheme using any implementation of the sliding window algorithm.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

/// A minimizer implementation based on a sliding window minimum algorithm.
/// Also takes a custom hash function, that defaults to FxHash.
pub struct SlidingWindowMinimizer<SlidingMinAlg, H = FxHash> {
    pub w: usize,
    pub k: usize,
    pub alg: SlidingMinAlg,
    pub hasher: H,
}

impl<H: Hasher, SlidingMinAlg: SlidingMin<H::Out>> Minimizer
    for SlidingWindowMinimizer<SlidingMinAlg, H>
{
    fn window_minimizers(&mut self, text: &[u8]) -> Vec<usize> {
        self.hasher
            // Iterate over k-mers and hash them.
            .hash_kmers(self.k, text)
            // (See source for the iterator 'extension trait'.)
            .sliding_min(self.w, &self.alg)
            .map(|elem| elem.pos)
            .collect()
    }
}

3.4 The queue

Let’s try to improve this somewhat inefficient \(O(n\cdot w)\) algorithm. A first idea is to simply keep a rolling prefix-minimum that tracks the lowest value seen so far. For example, when \(w=3\) and the hashes are \([10,9,8,7,6,…]\), the rolling prefix minimum is exactly also the minimum of each window. But alas, this won’t work for increasing sequences: when the hashes are \([1,2,3,4,5,…]\) and the window shifts from \([1,2,3]\) to \([2,3,4]\), the prefix minimum is \(1\) (at index \(0\)), but it is not a minimum of the new window anymore, which instead goes up to \(2\).

Nevertheless, as the window slides right, each time we see a new value that’s smaller than everything in the window, we can basically ‘forget’ about all existing values and just keep the new smaller value in memory. Otherwise, we can still forget about all values in the window larger than the value that is shifted in.

This is formalized by using a queue of non-decreasing values in the window. More precisely, each queue element will be smaller than all values in the window coming after it. At each step, the minimum of the window is the value at the front of the queue. Let’s look at our example again.⁷

We see that there are two reasons elements can be removed from the queue:

1. a smaller element is pushed,
2. the window shifts so far that the leftmost/smallest element falls outside it.

To handle this second case, we don’t just store values in this queue, but also their position in the original text, so that we can pop them when needed.

In code, it looks like this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

pub struct Queue;

impl<V: Ord + Copy> SlidingMin<V> for Queue {
    fn sliding_min(&self, w: usize, it: impl Iterator<Item = V>) -> impl Iterator<Item = Elem<V>> {
        // A double ended queue that is non-decreasing in both pos and val, so that the
        // smallest value is always at the front.
        let mut q = std::collections::VecDeque::<Elem<V>>::new();

        // Iterate over the items and their positions.
        it.enumerate().map(move |(pos, val)| {
            // Strictly larger preceding values are removed, so that the queue remains
            // non-decreasing.
            while q.back().is_some_and(|back| back.val >= val) {
                q.pop_back();
            }
            q.push_back(Elem { pos, val });
            // Check if the front falls out of the window.
            let front = q.front().expect("We just pushed");
            if pos - front.pos >= w {
                q.pop_front();
            }
            *q.front().expect("w > 0")
        })
    }
}

Performance characteristics

Since each element is pushed once and popped once, each call to push takes amortized constant \(O(1)\) time, so that the total runtime is \(O(n)\). We may expect a \(w=11\) times speedup over the buffered method, but most likely this will not be the case. This is due to the high number of unpredictable branches, which hurts pipelining/instruction level parallelism (ILP). Thus, while the complexity of the Queue algorithm is good, its practical performance may not be optimal.

3.5 Jumping: Away with the queue

While the queue algorithm is nice in theory, in practice it is not ideal due to the branch-misses. Thus, we would like to find some middle ground between keeping this queue and the naive approach of rescanning every window. A first idea is this: Once we find the minimizer of some window, we can jump to the first window after that position to find the next minimizer. This finds all distinct minimizers, although it does not compute the minimizer for each window individually. Thus, it solves Problem A instead of Problem B.⁷

Since the expected distance between random minimizers is \((w+1)/2\), we expect to scan each position just under two times.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

pub struct JumpingMinimizer<H = FxHash> {
    pub w: usize,
    pub k: usize,
    pub hasher: H,
}

impl<H: Hasher> Minimizer for JumpingMinimizer<H>
where
    H::Out: Ord,
{
    fn minimizer_positions(&mut self, text: &[u8]) -> Vec<usize> {
        let mut minimizer_positions = Vec::new();

        // Precompute hashes of all k-mers.
        let hashes = self.hasher.hash_kmers(self.k, text).collect::<Vec<_>>();

        let mut start = 0;
        while start < hashes.len() - self.w {
            // Position_min returns the position of the leftmost minimal hash.
            let min_pos = start
                + hashes[start..start + self.w]
                    .iter()
                    .position_min()
                    .expect("w > 0");
            minimizer_positions.push(min_pos);
            start = min_pos + 1;
        }
        // Possibly add one last minimizer.
        let start = hashes.len() - self.w;
        let min_pos = start + hashes[start..].iter().position_min().expect("w > 0");
        if minimizer_positions.last() != Some(&min_pos) {
            minimizer_positions.push(min_pos);
        }
        minimizer_positions
    }
}

Performance characteristics

One benefit of this method is that it is completely branchless. Each iteration of the loop corresponds to exactly one output. The only unpredictable step is which element of a window is the minimizer, since that determines which is the next window to scan. Thus, even though it is branchless, consecutive iterations depend on each other and ILP (instruction level parallelism, aka pipelining) may only have a limited effect.

3.6 Re-scan

I learned of an improved method via Daniel Liu’s gist for winnowing, but I believe the method is folklore¹⁰. This is a slightly improved variant of the method above. Above, we do a new scan for each new minimizer, but it turns out this can sometimes be avoided. Now, we will only re-scan when a previous minimum falls outside the window, and the minimum of the window goes up.⁷

Thus, we keep a ring buffer of the last \(w\) values, and scan it as needed. One point of attention is that we need the leftmost minimal value, but since the buffer is cyclic, the leftmost element in the buffer is not necessarily the leftmost one in the original text. Thus, we don’t only store elements but also their positions, and break ties by position.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

pub struct Rescan;

impl<V: Ord + Copy + Max> SlidingMin<V> for Rescan {
    fn sliding_min(&self, w: usize, it: impl Iterator<Item = V>) -> impl Iterator<Item = Elem<V>> {
        let mut min = Elem { val: V::MAX, pos: 0 };
        let mut ring_buf = RingBuf::new(w, min);

        it.enumerate()
            .map(move |(pos, val)| {
                let elem = Elem { val, pos };
                ring_buf.push(elem);
                min = min.min(elem);
                // If the minimum falls out of the window, rescan to find the new minimum.
                if pos - min.pos == w {
                    min = *ring_buf.iter().min().expect("w > 0");
                }
                min
            })
            .skip(w - 1)
    }
}

Performance characteristics

This method should have much fewer branch-misses than the queue, since we rescan a full window of size \(w\) at a time. Nevertheless, these re-scans occur at somewhat random moments, so some branch-misses around this are still to be expected.

3.7 Split windows

In an ideal world, we would use a fully predictable algorithm, i.e., an algorithm without any data-dependent branches/if-statements.¹¹

With some Googling¹² I found this codeforces blog on computing sliding window minima (see also this stackoverflow answer and an implementation in the moving_min_max crate). In the post, the idea is presented using two stacks. My explanation below goes in a slightly different direction, but in the end it comes down to the same.

Aside: Accumulating non-invertible functions. So, as we saw, most methods keep some form of rolling prefix minimum that is occasionally ‘reset’, since ‘old’ small elements should not affect the current window. This all has to do with the fact that the min operation is not invertible: knowing only the minimum of a set \(S\) is not sufficient to know the minimum of \(S - \{x\}\). This is unlike for example addition: sliding window sums are easy by just adding the new value to the sum and subtracting the value that falls out of the window.

This distinction is also exactly the difference between segment trees and Fenwick trees¹³: a Fenwick tree computes some function \(f\) on a range \([l, r)\) by removing \(f(a_0, \dots, a_{l-1})\) from \(f(a_0, \dots, a_{r-1})\), by decomposing both ranges into intervals with power-of-\(2\) lengths corresponding to the binary representation of their length. But this only works if \(f\) is invertible. If that is not the case, a segment tree must be used, which partitions the interval \([l, r)\) itself, without going ‘outside’ it.

The algorithm. Based on this, we would also like to partition our windows such that we can compute their minimum as a minimum over the different parts. It turns out that this is possible, and in a very clean way as well!

First, we conceptually chunk the input text into parts with length \(w\), as shown in the top row of Figure 5 for \(w=4\). Then, each window of length \(w\) either exactly matches such a chunk, or overlaps with two consecutive chunks. In the latter case, it covers a suffix of chunk \(i\) and a prefix of chunk \(i+1\). Thus, to compute the minimum for all windows, we need to compute the prefix and suffix minima of all chunks. We could do that as a separate pass over the data in a pre-processing step, but it can also be done on the fly:

1. Each time a new kmer is processed, its hash is written to a size-\(w\) buffer. (The orange values in the top-right of Figure 5.)
2. After processing the window corresponding to a chunk, the buffer contains exactly the hashes of the chunk. (Fourth row.)
3. Then, we compute the suffix-minima of the chunk by scanning the buffer backwards. (Fifth row.)
4. For the next \(w\) windows, we intialise a new rolling prefix minimum in the new chunk (green outline).
5. The minimum of each window is the minimum between the rolling prefix minimum in chunk \(i+1\) (green outline), and the suffix minimum in chunk \(i\) that we computed in the buffer (blue outline).⁷

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

pub struct Split;

impl<V: Ord + Copy + Max> SlidingMin<V> for Split {
    fn sliding_min(&self, w: usize, it: impl Iterator<Item = V>) -> impl Iterator<Item = Elem<V>> {
        let mut prefix_min = Elem { val: V::MAX, pos: 0 };
        let mut ring_buf = RingBuf::new(w, prefix_min);

        it.enumerate()
            .map(move |(pos, val)| {
                let elem = Elem { val, pos };
                ring_buf.push(elem);
                prefix_min = prefix_min.min(elem);
                // After a chunk has been filled, compute suffix minima.
                if ring_buf.idx() == 0 {
                    let mut suffix_min = ring_buf[w - 1];
                    for elem in ring_buf[0..w - 1].iter_mut().rev() {
                        suffix_min = suffix_min.min(*elem);
                        *elem = suffix_min;
                    }
                    prefix_min = Elem { val: V::MAX, pos: 0 };
                }
                prefix_min.min(ring_buf[ring_buf.idx()])
            })
            .skip(w - 1)
    }
}

Performance characteristics

The big benefit of this algorithm is that its execution is completely deterministic, similar to the jump algorithm: every \(w\) iterations we compute the suffix-minima, and then we continue again with the rolling prefix minimum. Since this ‘side loop’ is taken exactly once in every \(w\) iterations, the corresponding branch is easy to predict and does not cause branch-misses. Thus, we expect nearly \(0\%\) of branch-misses, and up to \(4\) instructions per cycle.

4 Analysing what we have so far

4.1 Counting comparisons

Before we look at runtime performance, let’s have a more theoretical look at the number of comparisons each method makes. We can do this by replacing the u64 hash type by a wrapper type that increments a counter each time .cmp() or .partial_cmp() is called on it. I exclude comparisons with the MAX sentinel value since those could possibly be optimized away anyway. The code is not very interesting, but you can find it here.

1
2
3
4
5
6

           w=10    w=20     formula
Buffered:  9.000  19.000    w-1
Queue:     1.818   1.905    2-2/(w+1)
Jumping:   1.457   1.595    ?
Rescan:    1.818   1.905    2-2/(w+1)
Split:     2.700   2.850    3-3/w

Observations:

• Buffered Does \(w-1\) comparisons for each window to find the max of \(w\) elements.
• Queue and Rescan both do exactly \(2-2/(w+1)\) comparisons per element.
• For Queue, each pop_back is preceded by a comparison, and each push has one additional comparison with the last element in the queu that is smaller than the new element. Only when the new element is minimal (probability \(1/(w+1)\)), the queue becomes empty and that comparison doesn’t happen. Also when the smallest element is popped from the front (probability \(1/(w+1)\)), that means it won’t be popped from the back and hence also saves a comparison.
• Unsurprisingly, the Jumping algorithm uses the fewest comparisons, since it also returns strictly less information than the other methods.
• The Split method, which seems very promising because of its data-independence, actually uses around \(50\%\) more comparisons (\(3-3/w\), to be precise), because for each window, a minimum must be taken between the suffix and prefix.

Open problem: Theoretical lower bounds

It would be cool to have a formal analysis showing that we cannot do better than \(2-2/(w+1)\) expected comparisons per element to compute the random minimizer of all windows, and similarly it would be nice to have a lower bound on the minimal number of comparisons needed to only compute the positions of the minimizers.

4.2 Setting up benchmarking

Adding criterion

Before we start our implementation, let’s set up a benchmark so we can easily compare them. We will use criterion, together with the cargo-criterion helper. First we add criterion as a dev-dependency to Cargo.toml.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

[dev-dependencies]
criterion = "*"

# Do not build the library as a benchmarking target.
[lib]
bench = false

[[bench]]
name = "bench"
harness = false

We add benches/bench.rs which looks like this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

        .warm_up_time(Duration::from_millis(500))
        .measurement_time(Duration::from_millis(2000))
        .sample_size(10);
    targets = initial_runtime_comparison,
        blog::counting::count_comparisons_bench,
        optimized, ext_nthash, buffered, local_nthash,
        simd_minimizer, human_genome
);
criterion_main!(group);

fn initial_runtime_comparison(c: &mut Criterion) {
    // Create a random string of length 1Mbp.
    let text = &(0..1000000)
        .map(|_| b"ACGT"[rand::random::<u8>() as usize % 4])
        .collect::<Vec<_>>();

    let hasher = FxHash;

    let w = 11;
    let k = 21;

    #[rustfmt::skip]
    let minimizers: &mut [(&str, &mut dyn Minimizer, bool)] = &mut [
        ("naive", &mut NaiveMinimizer { w, k, hasher }, true),
        ("buffered", &mut SlidingWindowMinimizer { w, k, alg: Buffered, hasher }, true),
        ("queue", &mut SlidingWindowMinimizer { w, k, alg: Queue, hasher }, true),
        ("jumping", &mut JumpingMinimizer { w, k, hasher }, false),
        ("rescan", &mut SlidingWindowMinimizer { w, k, alg: Rescan, hasher }, true),
        ("split", &mut SlidingWindowMinimizer { w, k, alg: Split, hasher }, true),
        ("queue_igor", &mut QueueIgor { w, k }, false),
        ("rescan_daniel", &mut RescanDaniel { w, k }, true),
    ];

    let mut g = c.benchmark_group("g");

Making criterion fast

By default, the benchmarks are quite slow. That’s for a couple of reasons:

• Before each benchmark, a 3 second warmup is done.
• Each benchmark is 5 seconds.
• After each benchmark, some plots are generated.
• At the end, some HTML reports are generated.

I’m impatient, and all this waiting really impacts my iteration time, so let’s make it faster:

• Reduce warmup time to 0.5s.
• Reduce benchmark time to 2s, and only take 10 samples.
• Disable plots and html generation.

The first two are done from code:

1
2
3
4
5
6
7
8
9

criterion_group!(
    name = group;
    config = Criterion::default()
        // Make sure that benchmarks are fast.
        .warm_up_time(Duration::from_millis(500))
        .measurement_time(Duration::from_millis(2000))
        .sample_size(10);
    targets = initial_runtime_comparison, count_comparisons_bench
);

The last is best done by adding the --plotting-backend disabled flag. For convenience, we add this rule to the justfile so we can just do just bench. I’m also adding quiet to hide the comparison between runs to simplify presentation.

1
2

bench:
    cargo criterion --plotting-backend disabled --output-format quiet

When we now do just bench, we see this:

A note on CPU frequency

Most consumer CPUs support turboboost to increase the clock frequency for short periods of time. That’s nice, but not good for stable measurements. Thus, I always pin the frequency of my i7-10750H to the default 2.6GHz:

1

sudo cpupower frequency-set --governor powersave -d 2.6GHz -u 2.6GHz

This usually results in very stable measurements.

Similarly, I have hyper threading disabled, so that each program has a full core to itself, without interference with other programs.

4.3 Runtime comparison with other implementations

Let’s compare the runtime of our method so far with two other Rust implementations:

• minimizer-iter is an implementation written by Igor Martayan. It returns an iterator over all distinct minimizers and uses the queue algorithm.
• Daniel Liu’s implementation of the rescan algorithm.

1
2
3
4
5
6
7
8

g/naive                 time:   [64.399 ms 64.401 ms 64.404 ms]
g/buffered              time:   [55.154 ms 55.160 ms 55.173 ms]
g/queue                 time:   [22.305 ms 22.322 ms 22.335 ms]
g/jumping               time:   [7.0913 ms 7.0969 ms 7.1042 ms]
g/rescan                time:   [14.067 ms 14.071 ms 14.078 ms]
g/split                 time:   [14.021 ms 14.023 ms 14.027 ms]
g/queue_igor            time:   [26.680 ms 26.689 ms 26.700 ms]
g/rescan_daniel         time:   [9.0335 ms 9.0340 ms 9.0349 ms]

Some initial observations:

• Measurements between samples are very stable, varying less than \(1\%\).
• Each method takes 7-55ms to process 1 million characters. That’s 7-55ns per character, or 7s to 55s for 1Gbp.
• Our queue implementation performs similar to the one in minimizer-iter.
• Daniel’s implementation of Rescan is around twice as fast as ours! Time to start looking deeper into our code.

Looking back at our expectations:

• As expected, Buffered is the slowest method by far, but even for \(w=11\), it is only \(2\times\) slower than the next-slowest method.
• Split is about as fast as Rescan, and even a bit slower. We were expecting a big speedup here, so something is weird.
• RescanDaniel is even much faster than our Split, so clearly there must be a lot of headroom.

4.4 Deeper inspection using perf stat

Let’s look some more at performance statistics using perf stat. We can reuse the benchmarks for this, and call criterion with the --profile-time 2 flag, which simple runs a benchmark for \(2\) seconds.

1
2
3

stat test='':
    cargo build --profile bench --benches
    perf stat -d cargo criterion -- --profile-time 2 {{test}}

The metrics that are relevant to us now are summarized in the following table.

Observe:

• Naive has an instructions per cycle (IPC) of over 4! I have never actually seen IPC that high. It’s not measurement noise, so most likely it comes from the fact that some instructions such as compare followed by conditional jump are fused.
• Naive also has pretty much no branch misses.
• Buffered only has half the IPC of Naive but is still slightly faster since it only hashes each kmer once. For some reason, there are a lot of branch misses, which we’ll have to investigate further.
• The two Queue variants also have relatively few instructions per cycle, and also relatively a lot of branch misses. Here, the branch-misses are expected, since the queue pushing and popping are inherently unpredictable.
• In fact, I argued that the last element of the Queue is popped with probability \(50\%\). But even then, the branch predictor is able to get down to \(4\%\) of branch misses. Most likely there are correlations that make these events somewhat predictable, and of course, there are also other easier branches that count towards the total.
• Jumping has an very good IPC of nearly \(4\), as we somewhat hoped for, and correspondingly has almost no branch misses.
• The remaining methods all have a decent IPC just above \(3\).
• RescanDaniel has much fewer branch misses than Rescan. We should investigate this.
• We argued that Split is completely predictable, but we see a good number of branch misses anyway. Another thing to investigate.

4.5 A first optimization pass

It’s finally time to do some actual optimizations! We’ll go over each of the algorithms we saw and optimize them. This step is ‘just’ so we can make a fair benchmark between them, and in the end we will only optimize one of them further, so I will try to keep it short. The final optimized version of each algorithm can be found in the repository, but intermediate diffs are only shown here.

Optimizing Buffered: reducing branch misses

Let’s have a look at perf record:

1
2
3
4
5
6

 1.16 │       cmp    %r8,%rbp          ; compare two values
 0.92 │       setne  %r11b
 0.70 │    ┌──jb     976               ; conditionally jump down
*6.15*│    │  mov    %r11b,%r13b       ; overwrite the minimum (I think)
 0.02 │    │  mov    %r13d,%r12d       ; overwrite the position (I think)
*5.44*│976:└─→mov    0x8(%rdx),%r10    ; continue

Most of these branches come from the fact that computing the minimum of two Elem { val: usize, pos: usize } is very inefficient, since two comparisons are needed, and then the result is conditionally overwritten.

The only reason we are storing an Elem that includes the position is that we must return the leftmost minimum. But since e.g. position_min() returns the leftmost minimum anyway, we can try to avoid storing the position and always iterate over values from left to right. We’ll add a method to our RingBuf for this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

        &mut self.data
    }
}

/// Extensions for ringbuf
impl<V> RingBuf<V> {
    pub fn forward_slices(&self) -> [&[V]; 2] {
        let (a, b) = self.data.split_at(self.idx);
        [b, a]
    }

    pub fn forward_min(&self) -> Elem<V>
    where
        V: Copy + Ord + Max,
    {
        let mut min = Elem { val: V::MAX, pos: 0 };
        for (idx, &v) in self
            .forward_slices()
            .into_iter()
            .flat_map(|part| part)
            .enumerate()
        {
            if v < min.val {
                min = Elem { val: v, pos: idx };
            }
        }
        min
    }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

pub struct BufferedOpt;

impl<V: Copy + Max + Ord> SlidingMin<V> for BufferedOpt {
    fn sliding_min(&self, w: usize, it: impl Iterator<Item = V>) -> impl Iterator<Item = Elem<V>> {
        // Note: this only stores V now, not Elem<V>.
        let mut ring_buf = RingBuf::new(w, V::MAX);
        it.enumerate()
            .map(move |(pos, val)| {
                ring_buf.push(val);
                let mut min = ring_buf.forward_min();
                min.pos += pos - w + 1;
                min
            })
            .skip(w - 1)
    }
}

The hot loop of the code looks much, much better now, using cmov instead of branches to overwrite the minimum.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

 1.59 │898:┌─→lea        0x8(%r14),%r13
 4.85 │    │  mov        %rbx,0xe0(%rsp)
 2.72 │    │  mov        %r9,0xf8(%rsp)
 5.29 │    │  mov        %r13,0xf0(%rsp)
 1.85 │    │  lea        0x1(%rax),%rbp
 4.78 │    │  mov        %rbp,0x110(%rsp)
11.20 │    │  mov        (%r14),%r14
 7.24 │    │  cmp        %rcx,%r14       ; compare
 3.41 │    │  cmovb      %rax,%r11       ; conditionally overwrite pos
 4.04 │    │  cmovb      %r14,%rcx       ; conditionally overwrite min
 2.50 │    │  mov        %rbp,%rax
 5.19 │    │  mov        %r13,%r14
 3.01 │8d4:│  test       %r14,%r14
 0.92 │    │↓ je         8de
 0.44 │    ├──cmp        %r9,%r14
 3.90 │    └──jne        898

Also the perf stat is great now: 4.23 IPC and 0.13% branch misses, so we’ll leave it at this.

Queue is hopelessly branchy

We can confirm that indeed a lot of time is spent on branch-misses around the popping of the queue.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

 1.20 │3a0:   lea        (%rdx,%rax,1),%rdi
 0.63 │ ↑     cmp        %rcx,%rdi
 1.34 │ │     mov        $0x0,%edi
 0.97 │ │     cmovae     %rcx,%rdi
 0.46 │ │     shl        $0x4,%rdi
 0.49 │ │     mov        %rsi,%r8
 1.58 │ │     sub        %rdi,%r8
 1.52 │ │  ┌──cmp        %rbp,(%r8)    ; compare the new element with the back
 1.97 │ │  ├──jb         3d2           ; if back is smaller, break loop
*9.85*│ │  │  dec        %rax          ; pop: decrease capacity
 3.70 │ │  │  mov        %rax,0x28(%rsp)
 0.41 │ │  │  add        $0xfffffffffffffff0,%rsi
 0.08 │ │  │  test       %rax,%rax     ; queue not yet empty
 0.37 │ └──│─ jne        3a0           ; try popping another element
 0.38 │    │  xor        %eax,%eax
*9.64*│3d2:└─→cmp        %rcx,%rax

Sadly, I don’t really see much room for improvement here.

Jumping is already very efficient

If we do a perf report here, we see that only 25% of the time is spent in the minimizer_positions function:

1
2

63.69%  <map::Map<I,F> as Iterator>::fold
24.25%  <JumpingMinimizer<H> as Minimizer>::minimizer_positions

The hot loop for finding the minimum of a window is very compact again.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

2.84 │170:┌─→lea     0x1(%rcx),%rsi
1.58 │    │  mov     -0x18(%rdx,%rcx,8),%rdi
4.17 │    │  mov     -0x10(%rdx,%rcx,8),%r8
0.09 │    │  cmp     %rdi,%rax
3.60 │    │  cmovb   %rax,%rdi
9.34 │    │  cmovbe  %rbp,%rsi
     │    │  lea     0x2(%rcx),%rax
0.94 │    │  cmp     %r8,%rdi
2.59 │    │  cmovae  %r8,%rdi
9.36 │    │  cmovbe  %rsi,%rax
1.35 │    │  mov     -0x8(%rdx,%rcx,8),%rsi
0.27 │    │  lea     0x3(%rcx),%rbp
0.86 │    │  cmp     %rsi,%rdi
5.69 │    │  cmovae  %rsi,%rdi
5.94 │    │  cmovbe  %rax,%rbp
1.48 │    │  mov     (%rdx,%rcx,8),%rsi
0.72 │    │  add     $0x4,%rcx
0.32 │    │  cmp     %rsi,%rdi
4.98 │    │  mov     %rdi,%rax
1.31 │    │  cmovae  %rsi,%rax
5.76 │    │  cmova   %rcx,%rbp
     │    ├──cmp     %rbx,%rcx
0.13 │    └──jne     170

More problematic is the 63% of time spent on hashing kmers, but we will postpone introducing a faster hash function a bit longer.

Optimizing Rescan

The original Rescan implementation suffers from the same issue as the Buffered and compares positions, so let’s reuse the RingBuf::forward_min() function, and compare elements only by value.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

pub struct RescanOpt;

impl<V: Ord + Copy + Max> SlidingMin<V> for RescanOpt {
    fn sliding_min(&self, w: usize, it: impl Iterator<Item = V>) -> impl Iterator<Item = Elem<V>> {
        let mut min = Elem { val: V::MAX, pos: 0 };
        // Store V instead of Elem.
        let mut ring_buf = RingBuf::new(w, V::MAX);

        it.enumerate()
            .map(move |(pos, val)| {
                ring_buf.push(val);
                if val < min.val {
                    min = Elem { val, pos };
                }
                // If the minimum falls out of the window, rescan to find the new minimum.
                if pos - min.pos == w {
                    min = ring_buf.forward_min();
                    min.pos += pos - w + 1;
                }
                min
            })
            .skip(w - 1)
    }
}

This has 3.55 IPC up from 3.02, and 0.91% branch misses, down from 2.26%, which results in the runtime going down from 14.1ms to 10.2ms.

The hottest part of the remaining code is now the branch-miss when if pos - min.pos == w is true, taking 6% of time, since this is an unpredictable branch. This seems rather unavoidable. Otherwise most time is spent hashing kmers.

Optimizing Split

Also here we start by avoiding position comparisons. This time it is slightly more tricky though: the RingBuf does have to store positions, since after taking suffix-minima, the value at each position does not have to correspond anymore to its original position.

The tricky part here is to ensure that the compiler generates a branchless minimum, since usually it creates a branch to avoid overwriting two values with their own value:

1
2
3
4
5
6
7
8

-if ring_buf[i + 1].val <= ring_buf[i].val {
-    ring_buf[i] = ring_buf[i + 1];
-};
+ring_buf[i] = if ring_buf[i + 1].val <= ring_buf[i].val {
+    ring_buf[i + 1]
+} else {
+    ring_buf[i]
+};

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

pub struct SplitOpt;

impl<V: Ord + Copy + Max> SlidingMin<V> for SplitOpt {
    fn sliding_min(&self, w: usize, it: impl Iterator<Item = V>) -> impl Iterator<Item = Elem<V>> {
        let mut prefix_min = Elem { val: V::MAX, pos: 0 };
        let mut ring_buf = RingBuf::new(w, prefix_min);

        it.enumerate()
            .map(move |(pos, val)| {
                let elem = Elem { val, pos };
                ring_buf.push(elem);
                if val < prefix_min.val {
                    prefix_min = elem;
                }
                // After a chunk has been filled, compute suffix minima.
                if ring_buf.idx() == 0 {
                    for i in (0..w - 1).rev() {
                        ring_buf[i] = if ring_buf[i + 1].val <= ring_buf[i].val {
                            ring_buf[i + 1]
                        } else {
                            ring_buf[i]
                        };
                    }
                    prefix_min = Elem { val: V::MAX, pos: 0 };
                }
                let suffix_min = ring_buf[ring_buf.idx()];
                if suffix_min.val <= prefix_min.val {
                    suffix_min
                } else {
                    prefix_min
                }
            })
            .skip(w - 1)
    }
}

Now, it runs in 11.2ms instead of 14.0ms, and has 3.46 IPC (up from 3.37 IPC) and 0.63% branch misses, down from 1.55%. The remaining branch misses seem to be related to the kmer hashing, although this is somewhat surprising as Jumping doesn’t have them.

4.6 A new performance comparison

Now that we have our optimized versions, let’s summarize the new results.

Note how BufferedOpt is now almost as fast as Queue, even though it’s \(O(nw)\) instead of \(O(n)\). Unpredictable branches really are terrible. Also, RescanOpt is now nearly as fast as RescanDaniel, although Jumping is still quite a bit faster. For reasons I do not understand SplitOpt is still not faster than RescanOpt, but we’ll have to leave it at this.

5 Rolling¹⁴ our own hash

As we saw, Jumping spends two thirds of its time on hashing the kmers, and in fact, SplitOpt and RescanOpt spend around half of their time on hashing. Thus, it is finally time to investigate some different hash functions.

5.1 FxHash

So far, we have been using FxHash, a very simple hash function originally used by FireFox and also used in the Rust compiler. For each u64 chunk of 8 bytes, it bit-rotates the hash so far, xor’s it with the new words, and multiplies by a fixed constant. For inputs less than 8 characters, it’s only a single multiplication. Even though this method iterates over the input length, all k-mers have the same length, so that the branches involved are usually very predictable. Nevertheless, hashing each k-mer involves quite some instructions, and especially a lot of code is generated to to handle all lengths modulo \(8\) efficiently.

WyHash

WyHash is currently the fastest cryptographically secure hash function in the SMHasher benchmark, but using it instead of FxHash slows down the benchmarks by roughly 2ns per character. The reason is that FxHash itself is not secure, and hence can be simpler.

5.2 NtHash: a rolling hash

We could avoid this for loop for the hashing of each k-mer by using a rolling hash that only needs to update the hash based on the character being added and removed. The classsic rolling hash is the Karp-Rabin hash, where the hash of a kmer \(x\) is \[ h(x) = \sum_{i=0}^{k-1} x_i \cdot p^i \pmod m, \] for some modulus \(m\) (e.g. \(2^{64}\)) and coprime integer \(p\). This can be computed in an incremental way by multiplying the previous hash by \(p\) and adding/subtracting terms for the character being added and removed.

NtHash (Mohamadi et al. 2016) is a more efficient variant¹⁵ of this based on ‘buzhash’, using multiplication and addition \((\times, +)\) uses bitwise rotation and xor \((\text{rot}, \oplus)\): \[ h(x) = \bigoplus_{i=0}^{k-1} \text{rot}^i(f(x_i)), \] where \(f\) is a static table lookup that maps each character to some fixed random value. (See also this post on ntHash).

The main benefit of this method over Karp-Rabin hashing is that typically table lookups and bit rotate instructions are much faster than multiplications.

One limitation of NtHash is that it is periodic: when \(k> 64\), characters that are exactly 64 positions apart will have the same rotation applies to their \(f(x_i)\). This means that the hashes of A....A and B....B are be the same, and similarly the hashes of A....B and B....A are the same. NtHash2 (Kazemi et al. 2022) fixes this by using multiple rotation periods, but we will stick with the original NtHash. In practice, minimizers with \(k> 64\) are rarely used anyway, and even if so, then the probability of hash collisions between the \(w\) k-mers in a window should be small enough.

The nthash crate

Using the nthash crate, the implementation of the Hasher trait that we’ve implicitly been using so far is simple:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

pub trait Hasher: Clone {
    type Out;
    fn hash(&self, t: &[u8]) -> Self::Out;
    #[inline(always)]
    fn hash_kmers(&mut self, k: usize, t: &[u8]) -> impl Iterator<Item = Self::Out> {
        t.windows(k).map(|kmer| self.hash(kmer))
    }
}

#[derive(Clone, Copy, Debug)]
pub struct FxHash;
impl Hasher for FxHash {
    type Out = u64;
    fn hash(&self, t: &[u8]) -> u64 {
        fxhash::hash64(t)
    }
}

#[derive(Clone, Copy, Debug)]
pub struct ExtNtHash;
impl Hasher for ExtNtHash {
    type Out = u64;
    fn hash(&self, t: &[u8]) -> u64 {
        nthash::ntf64(t, 0, t.len())
    }
    fn hash_kmers(&mut self, k: usize, t: &[u8]) -> impl Iterator<Item = Self::Out> {
        nthash::NtHashForwardIterator::new(t, k).unwrap()
    }
}