Practical Parallel Block Tree Construction

There exists a wide variety of different measures of compressibility. We refer to the excellent survey by Navarro [54] for a comprehensive overview. The empirical entropy [42] is based on the distribution of the characters (or substrings) of $T$. The $0$-th order entropy is defined as $H_0(T)={\sum }_{\alpha \in \mathrm{\Sigma }}\frac{n_{\alpha }}{n}\log \frac{n}{n_{\alpha }},$ where $n_{\alpha }={\mathrm{𝑟𝑎𝑛𝑘}}_{\alpha }(n+1)$. For $k>0$, the $k$-th order entropy is given by $H_k(T)={\sum }_{s\in {\mathrm{\Sigma }}^k}\frac{|T_s|}{n}{H}_0(T_s),$ where $T_s$ denotes the subsequence of characters occurring immediately after the substring $s\in {[1,\sigma ]}^k$ in $T$. For example, for $T=abcabd$, we have $T_{ab}=cd$. Overall, it holds that $\log \sigma \ge H_0(T)\ge H_1(T)\ge \mathrm{\cdots }\ge H_n(T).$ A major drawback of such statistical compression measures is that they fail to capture long-range repetitions, as they consider the text on a character-by-character basis.

Another class of measures is based on the size of a Lempel-Ziv factorization of a text. Such a factorization consists of $z$ factors, each of which can be stored as a reference to a previous occurrence in the text. Among the most prominent techniques in this class are Lempel-Ziv-based compression algorithms. The Lempel-Ziv 77 (LZ77) factorization [68] parses a text $T$ into factors $f_1,\mathrm{\dots },f_z\in {\mathrm{\Sigma }}^{+}$ such that $T=f_1\mathrm{\dots }{f}_z$ and, for each factor $f_i$ with $i\in [1,z]$, there exists a position $j\in [1,|f_1\mathrm{\dots }{f}_{i-1}|]$ in the already parsed prefix of the text such that $f_i=T[j\mathrm{\dots }j+|f_i|-1]$.

Strings can also be represented using context-free grammars [38]. Determining the size of the smallest grammar $g$ that generates only the given text is NP-complete. However, grammars of size $O\left(z\log \frac{n}{z}\right)$ can be computed in linear time [63]. The currently strongest known measure is the string complexity $\delta$ [39], defined as $\delta ={\max }_{k\in [1,n]}\frac{d_k}{k}$, where $d_k=|\{T[i..i+k-1]|i\in [0,n-k]\}|$. Interesting properties of $\delta$ include that it can be computed in linear time, that $z\in O(\delta \log \frac{n}{\delta })$, and that there exist families of length-$n$ strings that require $\mathrm{\Omega }(\delta \log \frac{n}{\delta })$ words to be encoded for every $n$ and every $\delta \in [2,n]$ [39]. Recently, it has been shown that, at least in theory, the most important string queries can be answered using $O(\delta \frac{n\log \sigma }{\delta \log n})$ words of space [37].

We analyze our algorithms in the PRAM model. In this model, multiple processors (PEs, processing elements) share a common memory. The PRAM model is synchronous, i.e., in each time step all PEs execute exactly one instruction. Several variants of the model exist, which differ in whether memory locations can be accessed exclusively or concurrently during a single time step. We further differentiate between read and write access.

The weakest variant is the EREW PRAM, which allows only exclusive read and exclusive write access to memory cells. A slightly stronger variant is the CREW PRAM model, where multiple PEs may concurrently read from the same memory cell, but only a single PE is allowed to write to a given memory cell during a time step. The strongest model considered here allows both concurrent read and concurrent write access (CRCW). We consider two variants of the CRCW model. In the Common-CRCW model, multiple PEs may write to the same memory cell simultaneously only if they all write the same value. In the Arbitrary-CRCW model, multiple PEs may write to the same memory cell without any restriction on the written values. In this case, it is unspecified which of the written values is stored.

When analyzing PRAM algorithms, we are primarily interested in their time and work. The time is the number of synchronous time steps executed by the algorithm, while the work is the total number of operations performed, including arithmetic on local data as well as memory reads and writes. The work thus corresponds to the running time of the algorithm when executed on a single PE. The differences between weaker (EREW) and stronger (CRCW) models become apparent when considering the time complexity of basic parallel primitives.

In addition to the All-Prefix operation, we frequently require parallel sorting as a basic subroutine.

We first describe the original block tree construction algorithm [6]. Note that there also exists a highly engineered construction algorithm based on fundamental toolbox data structures for compression [40]. However, the underlying algorithmic idea is the same and, at a high level, it simulates the procedure described below.

The block tree is constructed in two phases, followed by an optional pruning phase. In the first phase, we mark all leftmost occurrences of consecutive blocks $B_i\cdot B_{i+1}$. In the second phase, we compute the pointers (and offsets) from unmarked blocks to their leftmost occurrences on the same level. During the pruning phase, we might reduce the number of nodes in the tree by adding additional pointers (and offsets). This phase, however, does not improve the asymptotic size of the block tree.

First Phase:

In the first phase, we use a hash table to identify the leftmost occurrence of each pair of consecutive blocks. If a pair is already present in the hash table, it cannot correspond to the leftmost occurrence. Since we are not only interested in matches at block boundaries, we additionally scan the text to identify positions where such a pair may occur across block boundaries. To efficiently compare blocks, we use Karp-Rabin fingerprints [36]. Although matching fingerprints must still be verified to ensure that the corresponding substrings are equal, each substring needs to be verified only once. This results in a running time linear in the size of the level.

Second Phase:

In the second phase, we again use hash tables and Karp-Rabin fingerprints, this time focusing on unmarked blocks. For each such block, we compute the pointer (and offset) to its leftmost occurrence, which is contained within a pair of consecutive marked blocks. To this end, we store the fingerprints of all unmarked blocks in a hash table. We then scan the current level to locate the leftmost occurrences of these blocks using their fingerprints. As before, this can be done in time linear in the size of the level, since each block is considered only once.

Optional Pruning Phase:

Finally, there is an optional pruning step that can further reduce the size of the block tree in practice. However, this step does not yield any asymptotic improvement. Recall that, in the first phase, we mark all leftmost occurrences of consecutive pairs of blocks. This ensures that the leftmost occurrence of any unmarked block exists on the same level, possibly spanning one or two marked blocks. However, not all marked blocks are the target of a pointer. During the pruning step, these marked but untargeted blocks are removed. We do not consider this step in our parallel algorithm described below, as it does not improve the asymptotic space bounds. Nevertheless, it is included in our implementation; see Section 5.4.

The running time depends on the sizes of the levels of the block tree. Let $b_k$ denote the block size on the $k$-th level and let $c_k$ denote the number of blocks on that level. Thus, $b_k\cdot c_k$ characters must be processed on level $k$. According to Lemma 5, there are at most $3z\tau$ blocks per level.

For levels starting from $\mathrm{\ell }\ge 1+{\log }_{\tau }\frac{3z}{s}$, we can bound the total number of characters processed across all levels from $\mathrm{\ell }$ to the last level $h$ as ${\sum }_{m=\mathrm{\ell }}^h{b}_m\cdot c_m\le {\sum }_{m=\mathrm{\ell }}^h\frac{n}{{\tau }^{m-\mathrm{\ell }}}=n{\sum }_{m=\mathrm{\ell }}^h\frac{1}{{\tau }^{m-\mathrm{\ell }}}\le \frac{n}{1-\frac{1}{\tau }}=O(n)$. All preceding levels also contain at most $O(n)$ characters each. Since each character can be processed in constant expected time with high probability, we obtain the following running time bound.

Our first parallel block tree construction algorithm replaces hash tables with sorting. Whenever fingerprints are required, we compute them in parallel.

Since we require the fingerprints of all pairs of blocks in the first phase and of all substrings of length $b_{\mathrm{\ell }}$ in the second phase on each level $\mathrm{\ell }$, we apply Lemma 7.

The space required to store these fingerprints constitutes the main memory bottleneck of the construction algorithm, as this storage is needed on every level (not at the same time).

We can now efficiently compute the fingerprints, but we still need to validate all comparisons. This verification can be performed in $O(\log n)$ time in the EREW and CREW models and in $O(1)$ time in the CRCW models. In both cases, the verification requires $O(n)$ work and no additional space.

Next, we sort tuples consisting of a fingerprint and the corresponding block position. In the sorted order, equal fingerprints appear consecutively. Using Lemmas 3 and 4, this sorting step can be performed in $O(\log n)$ time in the EREW model and $O\left(\frac{\log n}{\log \log \log n}\right)$ time in the CRCW model. Once the fingerprints are sorted (first by fingerprint and then by position), the leftmost occurrences can be identified in $O(1)$ time.

The second phase is slightly more involved. Here, we again use sorting to identify, for each unmarked block, the block and position to which it must point. To this end, we combine the fingerprints obtained during the scan with the fingerprints of the unmarked blocks. We sort these fingerprints as tuples consisting of the fingerprint value and the corresponding text position. For unmarked blocks, the text position additionally encodes that the entry corresponds to an unmarked block. After sorting, equal fingerprints are grouped together and ordered by text position, with the unmarked blocks placed at the end of each group.

We then construct an array $\mathrm{𝑜𝑐𝑐}$ that helps us compute the pointers and offsets. It has an entry for every fingerprint we just sorted. Then, for each fingerprint, we write a $0$ to $\mathrm{𝑜𝑐𝑐}$ if it is equal to the fingerprint immediately to its left, and otherwise write its text position. For the rightmost occurrence of a fingerprint, we instead write the negated text position of the first occurrence to $\mathrm{𝑜𝑐𝑐}$.²²2Negated, i.e., $-\mathrm{𝑝𝑜𝑠}$ for a position $\mathrm{𝑝𝑜𝑠}$ and not bit-wise negated, where we flip each bit individually. After computing the prefix sum over $\mathrm{𝑜𝑐𝑐}$, the array contains the correct text positions for all but the rightmost occurrence of each fingerprint group (for which the correct value is found to its left). The prefix sum can be computed in $O(\log n)$ time in the EREW model and in $O\left(\frac{\log n}{\log \log n}\right)$ time in the Common-CRCW model.

This procedure allows us to compute every level of the block tree using only sorting and prefix sums. In the description above, we implicitly assume that we operate on $O(n)$ objects throughout. However, on each level of the block tree, some blocks may be replaced by references.

A major drawback of the algorithm described above is its high working-memory requirement. We address this issue by employing domain decomposition to introduce a trade-off between scalability and memory usage. The idea of domain decomposition is to partition the input and process each partition as independently as possible.

For block trees, however, each level depends heavily on the previous one. Therefore, we have to partition the input on each level, after combining the results of the previous level. The algorithm from Lemma 7 remains largely unchanged, with the key difference that we now apply a local pre-filtering step.

Our algorithm uses $K$ partitions. For each partition, we maintain a hash table that acts as a local filter. Only the locally leftmost occurrences of block pairs (and, in the second phase, substrings encountered during the scan) are forwarded to the sorting step. This significantly reduces the memory overhead, as we now operate on only a $K$-approximation of the blocks required for the construction.

We accelerate the computation of Karp-Rabin fingerprints using SIMD techniques. Our implementation was initially developed for the Intel AVX2 instruction set extension and later ported to AArch64 ASIMD. We use the parameters $q=2^{32}$, $r=33$, and $\sigma =256$, similar to the well-known djb2 hash function.

Since AVX2 vectors are 32 bytes wide, they can store eight 32-bit hash values. In each iteration of our function, 16 new hashes $\varphi (s..e),\varphi (s+1,e+1),\mathrm{\dots },\varphi (s+15,e+15)$, stored in two vectors, are computed from the previous 16 hashes according to the scheme

The 16-tuple hashes $\varphi (e-15..e)$ to $\varphi (e..e+15)$ are computed in two phases. First, hashes of 4-tuples are computed as shown in Figure 2. A chunk of input (a) is loaded into a SIMD register and permuted to form sliding windows (b) of four characters each. Using a data-parallel pairwise multiply–add instruction, pairs of characters $(T[i],T[i+1])$ are multiplied with $(r,1)$ to produce pair hashes (c) according to $\varphi (i,i+1)=T[i]\cdot r+T[i+1]$. This phase is then repeated by multiply–adding pairs of pair hashes with $(r^2,1)$, yielding the 4-tuple hashes (d). Since , the pair hashes fit into 16 bits, preserving the vector width throughout steps (a) to (d).

These 4-tuple hashes are then combined into 8-tuple hashes and finally into 16-tuple hashes as depicted in Figure 3, using element-wise operations. Some intermediate values are carried over from the previous iteration to reduce the number of required operations. Vectors of 4-tuple hashes offset by 4, which are needed to compute the 8-tuple hashes, are formed by combining the rear half of the previous vector with the front half of the next vector.

The 16-character hashes are cached in a ring buffer so that they can be subtracted once the start of the sliding window reaches the current end. For windows of 17–32 characters, a variant of the above algorithm is used in which the ring buffer remains in vector registers to avoid costly store-reload latency. For hashes of up to 16 characters, specialized routines directly compute each hash from the corresponding characters.

The total computational effort per iteration consists of 4 pairwise multiply-add operations, 8 multiplications, 8 additions, and 6 permutations. However, due to the high degree of data parallelism and the latency of 32-bit SIMD multiplications on current Intel microarchitectures, throughput is nevertheless latency-bound at approximately 11 cycles per 16 characters (10 cycles to compute the product $\varphi (s-16..e-16)\cdot r^{16}$ and 1 cycle to add the remaining terms) [1]. This compares favorably with the $\sim 2.5$ cycles per character required to compute Karp–Rabin fingerprints using the conventional recurrence [36]

and yields an overall performance improvement of approximately $5\%$ across all tested configurations of our algorithm.

We implemented our block tree in C++, following the same interface as previous block tree implementations [6, 40]. A practical block tree implementation requires only one bit vector per level (to mark inner nodes), along with pointers and offsets to earlier occurrences for inner nodes. We implemented the memory-efficient algorithm based on domain decomposition presented in Section 5.2, since the number of available PEs is determined by the hardware and is typically much smaller than the input size. Preliminary experiments showed that using a number of partitions different from the number of PEs negatively impacts performance.

In our implementation, we do not perform sorting. Although sorting simplifies the theoretical analysis, it is unnecessary in practice. Instead, each PE is responsible for a local partition of the current level. Blocks are initially marked only with respect to occurrences within the local partition; for this purpose, we use a hash table (ankerl::unordered_dense::map) as a filter to identify locally leftmost occurrences.

These marked blocks are then communicated to the responsible PE. To this end, we use a wait-free multiple-producer single-consumer (MPSC) queue [2]. Since fingerprints are uniformly distributed, they are used to assign blocks to PEs via a modulo operation. In practice, this results in a well-balanced workload without noticeable skew. Once all blocks have been marked globally, we proceed with the second phase in the same manner.

Another observation is that, beyond a certain threshold, the size of the MPSC queue does not influence performance. In our experimental evaluation (Section 6), we use a queue capacity of 512 fingerprints. Smaller queues led to load imbalances between PEs, while larger queues did not yield any measurable performance improvement.

We first evaluate the scalability of all implementations. To this end, we conducted a strong-scaling experiment in which each PE is mapped to a single physical CPU core. The throughput, measured as processed input in MiB per second, is shown in Figure 4.

Overall, all tested algorithms exhibit similar behavior across inputs, largely independent of input size, alphabet size, and average Lempel-Ziv factor length. We include Par-LZ to illustrate the inherent difficulty of efficiently computing the LZ77 factorization in parallel. Consequently, Par-LZ is the slowest of all evaluated algorithms. Interestingly, it scales similarly to Par-BT, despite Par-BT computing an LZ77-style approximation together with additional information required for query support rather than the LZ77 factorization itself.

Among the general-purpose compression tools, pigz is the fastest and best-scaling algorithm overall. However, this performance comes at the cost of a weaker compression ratio compared to pxz, the parallel LZMA compressor included in our benchmark. Overall, pxz achieves performance comparable to Par-BT(on some inputs it is faster, on others similar, and on some slower). Notably, it does not scale as well as pigz. We emphasize again that these compressors cannot be used for block tree construction, nor does the data they produce support access, rank, or select queries.

Currently, LPF-BT is the state-of-the-art block tree construction algorithm in terms of construction speed. Its authors also provide a naive parallelization, which, to the best of our knowledge, is the only other parallel block tree construction approach. This parallelization relies on parallel algorithms for toolbox compression data structures. Since only the construction of these toolbox data structures is parallelized, rather than the entire algorithm, LPF-BT exhibits limited scalability, with little to no speedup beyond 8 PEs.

In contrast, our algorithm Par-BT achieves speedups of up to 15.5 (on einstein.de) and an average speedup of 11.1 when using 64 PEs. While this is far from linear speedup, it is consistent with expectations given the difficulty of the underlying problem, the LZ77-like processing with its data dependencies and the scalability behavior observed for Par-LZ.

Finally, even when using a single PE, the achievable memory throughput is significantly higher than the block tree construction throughput. Specifically, memory throughput exceeds 1 GiB/s on a single thread, which is several orders of magnitude higher than the construction throughput. This indicates that additional memory bandwidth from multiple threads provides little benefit and that the algorithms are CPU-bound.

In addition to construction speed, we are interested in the amount of memory required during construction. Excessive working-memory demands can severely limit the applicability of these algorithms in memory-constrained environments, such as HPC clusters or shared hardware. We report the memory allocated on the heap as a percentage of the input size in Figure 5. To the best of our knowledge, none of the evaluated implementations allocate significant memory on the stack, so considering heap usage alone is sufficient.

Although very fast, LPF-BT requires a significant amount memory, as it relies on several auxiliary data structures that need $O(n\log n)$ bits for an input of length $n$. A similar observation holds for Par-LZ. Interestingly, the working-memory requirement of LPF-BT is independent of the number of PEs used. As expected, Par-BT requires significantly less space than LPF-BT; its space overhead is roughly an order of magnitude smaller. Among all evaluated tools, Par-LZ exhibits the highest memory consumption during construction.

Both pigz and pxz are comparatively memory-efficient, but their memory usage increases proportionally with the number of PEs. Notably, for small PE counts, pigz does not require any heap allocations. However, as discussed earlier, this also results in weaker compression ratios compared to pxz. We observe slight variations in the space overhead of Par-BT depending on the number of PEs. This behavior can be explained by the local hash tables, whose sizes depend on the distribution of fingerprints. Since the used hash function is only pseudo-random, perfect load balancing across partitions cannot be guaranteed. Still, the observed distribution is sufficiently balanced to keep hash table sizes comparable across PEs.