Bounding the Average Move Structure Query for Faster and Smaller RLBWT Permutations

An array $A$ of $|A|=n$ elements is represented as $A[0..n-1]$. $\{A\}$ is the set derived from $A$ of size $|\{A\}|$. For $x\in \mathbb{N}$, let $[0,x]$ represent the array $[0,1,\mathrm{\dots },x]$. A predecessor query on a set $B$ with element $x$ is defined inclusively such that $B.pred(x)=\max \{y\in B|y\le x\}$. When dealing with a set of distinct elements, $B.rank(x)$ refers to the index of $B.pred(x)$ in sorted order. A string $S[0..n-1]$ is an array of symbols drawn from an ordered alphabet $\mathrm{\Sigma }$ of size $\sigma$. We use $S[i..j]$ to represent a substring $S[i]\mathrm{\dots }S[j]$. For strings $S_1$, $S_2$ let $S_1{S}_2$ represent their concatenation. A prefix of $S$ is some substring $S[0..j-1]$, where a suffix is some substring $S[i..n-1]$. The length of the longest common prefix ($lcp$) of $S_1$,$S_2$ is defined as $lcp(S_1,S_2)=\max \{j|S_1[0..j-1]=S_2[0..j-1]\}$. A text $T[0..n-1]$ is assumed to be terminated by the special symbol $\$\notin \mathrm{\Sigma }$ of least order so that suffix comparisons are well defined. We use the RAM word model, assuming machines words of size $𝒲=\mathrm{\Theta }(\log n)$ with basic arithmetic and logical bit operations in $O(1)$-time.

The suffix array ($\mathrm{SA}$) [8] of a text $T[0..n-1]$ is a permutation $\mathrm{SA}[0..n-1]$ of $[0,n-1]$ such that $\mathrm{SA}[i]$ is the starting position of the $i$th lexicographically smallest suffix of $T$. Let a collection of documents $𝒟=\{T_1,T_2,\mathrm{\dots },T_d\}$ be represented by their concatenation $T[0..n-1]=T_1{T}_2\mathrm{\dots }{T}_d\$$. The document array ($\mathrm{DA}$) [9] is defined as $\mathrm{DA}[0..n-1]$ where $\mathrm{DA}[i]$ stores which document $T[\mathrm{SA}[i]..n]$ begins in. We use $d$ to refer to the number of documents in a collection. The Burrows-Wheeler Transform ($\mathrm{BWT}$) [6] is a permutation $\mathrm{BWT}[0..n-1]$ of $T[0..n-1]$ such that $\mathrm{BWT}[i]=T[\mathrm{SA}[i]-1]$ if $\mathrm{SA}[i]>0$, otherwise $\mathrm{BWT}[i]=T[n-1]=\$$. Let $r$ be the number of maximal equal character runs of the $\mathrm{BWT}$. The run-length encoded BWT ($\mathrm{RLBWT}$) is an array $\mathrm{RLBWT}[0..r-1]$ of tuples where $\mathrm{RLBWT}[i].c$ is the character of the $i$th $\mathrm{BWT}$ run and $\mathrm{RLBWT}.\mathrm{\ell }$ its length.

The $\mathrm{BWT}$ is reversible by using the last-to-first ($\mathrm{LF}$) mapping $\mathrm{LF}(i)$, a permutation over $[0,n-1]$ satisfying $\mathrm{SA}[\mathrm{LF}(i)]=(\mathrm{SA}[i]-1)modn$. The first-to-last ($\mathrm{FL}$) mapping is its inverse such that $\mathrm{FL}(i)$ satisfies $\mathrm{SA}[\mathrm{FL}(i)]=(\mathrm{SA}[i]+1)modn$. The permutation $\varphi$ over $[0,n-1]$ is defined such that $\varphi (\mathrm{SA}[i])=\mathrm{SA}[(i-1)modn]$, returning the $\mathrm{SA}$ value for the suffix of preceding lexicographic rank. Its inverse ${\varphi }^{-1}$ is defined symmetrically such that ${\varphi }^{-1}(\mathrm{SA}[i])=\mathrm{SA}[(i+1)modn]$. The longest common prefix array ($\mathrm{LCP}$) stores the $lcp$ between lexicographically adjacent suffixes in $T[0..n-1]$, i.e., $\mathrm{LCP}[0..n-1]$ such that $\mathrm{LCP}[0]=0$ and $\mathrm{LCP}[i]=lcp(T[\mathrm{SA}[i-1]..n-1],T[\mathrm{SA}[i]..n-1])$ for $i>0$.

For a permutation $\pi$ over $[0,n-1]$, let $S$ be the set of $r$ positions $i$ such that either $i=0$ or $\pi (i-1)\ne \pi (i)-1$. Then given any $i\in [0,n-1]$, we can evaluate $\pi (i)$ using $i^{\prime }=S.pred(i)$ as

and can do so in $O(r)$-space and $O(\log \log n)$-time by storing $S_{\pi }=\{\pi (i)\mid \forall i\in S\}$ alongside a compressed sparse bitvector representing $S$ supporting predecessor queries. Intuitively, for any $j\in [0,r-1]$, there exists a disjoint interval from $S[j]$ with length $\mathrm{\ell }=S[j+1]-S[j]$ (or $\mathrm{\ell }=n-S[j]$ if $j=r-1$) such that the $j^{\text{th}}$ interval corresponds to $[S[j],S[j]+\mathrm{\ell })$ within which all positions permute contiguously.

Nishimoto and Tabei’s move structure [10] is an alternative formulation in $O(r)$-space avoiding predecessor queries which, assuming we are given the rank of the predecessor for position $i$ in $S$, achieves constant time per step when computing consecutive permutations ${\pi }^k(i)$ for $k\ge 1$ . Assume, for some $i$, that we know the rank $j^{\prime }$ of its predecessor $i^{\prime }=S.pred(i)$ in $S$. Then we compute $\pi (i)$ as in (1) by accessing $i^{\prime }=S[j^{\prime }]$ and $\pi (i^{\prime })=S_{\pi }[j^{\prime }]$. However, to compute ${\pi }^2(i)=\pi \left(\pi (i)\right)$ we would need a new predecessor query to find the rank $j^{\prime \prime }$ of its predecessor $i^{\prime \prime }=S.pred(\pi (i))$ in $S$. If we also store the position of predecessors of $S_{\pi }$ in $S$, $S_{rank(\pi )}=\{S.rank(\pi (i))\mid \forall i\in S\}$, then we could scan down starting from $S_{rank(\pi )}[j^{\prime }]$ in $S$ until we find the largest rank $j^{\prime \prime }$ such that $S[j^{\prime \prime }]\le \pi (i)$. Then, we compute ${\pi }^2(i)$ as we did before using accesses with $j^{\prime \prime }$, and continue similarly for any ${\pi }^k(i)$ with $k>1$. A single iteration of this procedure is called a move query. Figure 1 shows the relationship between a permutation and its move structure.

The move query time is proportional to what Zakeri et al. [14] call fast forwards, the distance from $S_{rank(\pi )}[j^{\prime }]$ to $j^{\prime \prime }$. As described, there could be $O(r)$ fast forwards in the worst case. Brown, Gagie, and Rossi [5] showed that even over all distinct move queries, ${\pi }^n(i)$ for any $i$ when formed from a single cycle, the worst case total number of fast forwards is $\mathrm{\Theta }(n\cdot r)$ and hence average $\mathrm{\Theta }(r)$-time per distinct query. However, Nishimoto and Tabei [10] proved in their original result that we can balance the intervals by arbitrarily splitting some of them. Let $S^{\prime }\supseteq S$ also contain the positions of the new intervals introduced by their balancing, with $|S^{\prime }|=r^{\prime }$. Then the number of fast forwards for any move query is $O(1)$, and $r^{\prime }\in O(r)$. Brown, Gagie, and Rossi [5] showed how to parameterize balancing for a fixed parameter $\alpha$ into at most $r^{\prime }=r+\frac{r}{\alpha -1}$ intervals and less than $2\alpha$ fast forwards for any step. Bertram, Fischer, and Nalbach [1] gave a practical $O(r\log r)$-time and $O(r)$-space balancing algorithm for parameter $\alpha$. We call this version with theoretical guarantees a balanced move structure, in contrast to the original unbalanced move structure.

In practice, move structures often do not store $S$, $S_{\pi }$ and $S_{rank(\pi )}$ exactly as described, but in some $O(r)$-space representation of the components which support access for an arbitrary rank $j$. Brown, Gagie, and Rossi [5] suggested a practical improvement by switching to relative positions, replacing $S$ with $S_{\mathrm{\ell }}$, where $S_{\mathrm{\ell }}[j]$ is the length of the $j^{\text{th}}$ interval. Using relative positions, a move query now returns the interval and offset within that interval rather than the absolute value; however, representing $S_{\mathrm{\ell }}$ uses less space than representing $S$ naively in practice, since $S_{\mathrm{\ell }}$ is the delta encoding of $S$. Similarly, they noted that, whether using absolute or relative positions, you could replace $S_{\pi }$ with $S_{\mathrm{\Delta }}$, where $S_{\mathrm{\Delta }}[j]$ is the offset of $S_{\pi }[j]$ from its predecessor in $S$.

Brown, Gagie, and Rossi [5] then introduced splitting intervals with respect to a maximum length to practically reduce the number of fast forwards, experimenting with values as a function of the average interval length $\frac{n}{r}$. Bertram, Fischer, and Nalbach [1] (and later, Zakeri et al. [15]) leveraged length splitting using an arbitrary constant to instead reduce the size of representing any component in $S_{\mathrm{\ell }}$ and/or $S_{\mathrm{\Delta }}$. We show that a special type of length splitting, which we call length capping, guarantees additional properties with respect to move query time.

Consider length capping with a fixed parameter $c$, which splits intervals of a move structure such that the maximum length of any interval is $L=c\cdot \frac{n}{r}$, a constant factor of the average interval length. Let the amortized complexity of a move structure refer to supporting $n$ move queries consecutively from an arbitrary starting position of permutations formed from a single cycle, which is equivalent to performing each of the $n$ distinct move queries. Length splitting in this way gives an average of at most $c+1$ fast forwards across distinct move queries while introducing at most $\frac{r}{c}$ new intervals:

A move structure maintaining $O(1)$-time move queries can be represented by $S$ as a sparse bitvector using $r\log \frac{n}{r}+O(r)$-bits, $S_{\pi }$ in $r\log n$ bits, and $S_{rank(\pi )}$ in $r\log r$ bits. Overall, this representation results in $O(r\log n)$-bits. Although Brown et al. [4] showed that move structures for $\mathrm{LF}$ and $\mathrm{FL}$ can be represented in $O(r\log \frac{n}{r})$-bits, we do not know how to remove the $r\log n$ bit components for general permutations, even when switching $S_{\pi }$ for $S_{\mathrm{\Delta }}$.

However, after length capping, since any element in $S_{\mathrm{\Delta }}$ is the offset within an interval, it is upper bounded by the maximum length $L\in O(\frac{n}{r})$. Thus, with length splitting, we can support both absolute and relative positions by representing $S_{\mathrm{\Delta }}$ in $O(r\log \frac{n}{r})$-bits. If using relative positions, since any element in $S_{\mathrm{\ell }}$ is also bounded by $L$, we can represent it in the same $O(r\log \frac{n}{r})$-bits as representing $S$ using a sparse bitvector. This may be useful in practice since array access is cache friendly compared to sparse bitvector select queries.

Overall, this results in $O(r\log r+r\log \frac{n}{r})=O(r\log n)$-bits, so the asymptotic bound is unchanged. However, no single component requires $O(r\log n)$-bits. Let $\widehat{r}=\mathrm{\Theta }(r)$ to suppress constant factors of $r$ introduced by balancing and/or length capping. Then the total space in bits by applying length capping changes from

to

effectively saving $O(r\log r)$-bits of space. This result can also be applied by length capping alongside balancing.

A major benefit of length capping is the $O(r)$-time construction in comparison to balancing with $O(r\log r)$-time. However, a technicality of Theorem 2 is that it specified being given an unbalanced move structure, required since finding the move structure itself takes superlinear time with respect to $r$. That is, given $S$ and $S_{\pi }$, finding $S_{rank(\pi )}$ in $O(r)$-space requires computing an $O(r\log r)$-time sorting of $S_{\pi }$ or computing the required predecessor queries in $O(r\log \log n)$-time using a sparse bitvector. Based on known lower bounds for comparison based sorting and predecessor queries, it would appear that this problem cannot be solved in $O(r)$-time and space for an arbitrary permutation.

However, various permutations defined by the $\mathrm{RLBWT}$ do benefit from the improved time complexity to construct a length capped move structure. Both $\mathrm{LF}$ and $\varphi$ as well as their respective inverses $\mathrm{FL}$ and ${\varphi }^{-1}$ can be represented as $O(r)$-space move structures [10], where $r$ is the number of runs in the $\mathrm{BWT}$. For $\mathrm{LF}$, it easily follows that $\mathrm{BWT}[i]=\mathrm{BWT}[i+1]$ implies $\mathrm{LF}(i+1)=\mathrm{LF}(i)+1$ and thus its permutation has $r$ contiguously permuted runs. For $\varphi$, given the properties of $\mathrm{LF}$, it follows that

and thus there are at most $r$ contiguously permuted runs in the permutation of $\varphi$ [7]. Further, these $\mathrm{RLBWT}$ based permutations are formed from a single cycle, allowing us to apply the amortized results of Corollary 3.

The properties of $\mathrm{RLBWT}$ based permutations also permit alternative constructions of their corresponding move structure. For $\mathrm{LF}$ and $\mathrm{FL}$ it is possible to leverage the properties of the $\mathrm{BWT}$ to find $S_{rank(\pi )}$ in $O(r)$-time and space. This is since the $\mathrm{LF}$ permutation is determined by the lexicographic order and individual rank of its $\mathrm{BWT}$ characters: once the $\mathrm{RLBWT}$ has been scanned and $S_{\pi }$ computed, we can obtain $S_{rank(\pi )}$ by iterating over the intervals by order of their corresponding output as computed by lexicographic rank of its $\mathrm{RLBWT}$ character [3]. This allows us to apply the result of Theorem 2 to derive the construction time of a length capped move structure for these permutations.

For $\varphi$ and ${\varphi }^{-1}$, the best bound in terms of $r$ is still $O(r\log r)$-time and $O(r)$-space, which is the complexity of finding the unbalanced move structure of an arbitrary permutation. However, Sanaullah et al. [11] noted that you could iterate over the text using $\mathrm{LF}$ to sample $\mathrm{SA}$ positions using $O(n)$ consecutive move queries to derive $S_{rank(\pi )}$ in $O(n)$-time and $O(r)$-space. Although $O(r\log r)$-time may often be preferred to $O(n)$-time in practice, the latter is useful when constructing a length capped move structure for the purpose of performing $n$ consecutive queries as described in Corollary 3 so that the final complexity depends only on $n$. Their method assumes a balanced move structure and thus overall they require $O(n+r\log r)$-time when including time to balance. By applying Lemma 6, we achieve new bounds for construction of an unbalanced move structure and a length capped move structure for $\varphi$ and ${\varphi }^{-1}$.

To further illustrate the benefit of length capped move structures for scenarios where query time is amortized $O(1)$-time over $n$ consecutive steps, we apply the results of Section 4.1 to known $\mathrm{RLBWT}$ algorithms. The most basic query involving a $\mathrm{BWT}$ is inversion, which uses $\mathrm{LF}$ to iterate over every character in $\mathrm{BWT}[0..n-1]$ to return the original text $T[0..n-1]$. Such queries involving traversal using $\mathrm{LF}$ or $\mathrm{FL}$ could be considered iterating in text order. However, it is not known how to do this in optimal $O(n)$-time and $O(r)$-space since obtaining $O(1)$-time $\mathrm{LF}$ requires $O(r\log r)$-time balancing (Theorem 1). Using a length capped move structure, we obtain the optimal-time for inversion using only $O(r)$-space in addition to the output $T[0..n-1]$.

We explore length capping through the Orbit move structure library featuring the results of Theorem 2, available at https://github.com/drnatebrown/orbit. Inspired by the optimizations of Bertram et al.’s Move-r [1] and Depuydt et al.’s b-move, our move structure is built over a bit packed matrix. That is, we use bitpacking to support the components of the move structure as columns which use the minimum possible fixed bitwidth, while keeping rows contiguous in memory to support cache efficiency. This is essential to length capping, since it gives the maximum benefit to the bounded size of any elements in $S_{\mathrm{\ell }}$ or $S_{\mathrm{\Delta }}$. Most importantly, we accept the length capping parameter $c$ as input to construction.

Beyond a proof of concept for length capping, Orbit is designed as a header only library that can be inserted into existing software. As such, it supports various representations of a move structure. For absolute positions, the components $S$, $S_{\mathrm{\Delta }}$ and $S_{rank(\pi )}$ are stored in a bitpacked matrix using roughly $r\log r+2r\log n$ bits (we do not use sparse bitvectors, optimizing instead for speed). Alternatively, relative positions are supported by storing the components $S_{\mathrm{\ell }}$, $S_{\mathrm{\Delta }}$, and $S_{rank(\pi )}$ in roughly $r\log r+2r\log \frac{n}{r}$ bits as described in Section 3.3. Move queries support using either a traditional linear scan over fast forwards or exponential search as in Corollary 5.

Within the library are specialized move structures for supporting the $\mathrm{RLBWT}$ based permutations $\mathrm{LF}$, $\mathrm{FL}$, $\varphi$, and ${\varphi }^{-1}$. These form the basis needed to perform queries such as count and locate as originally described by Nishimoto and Tabei [10] when using move structures to construct a full text index for pattern matching. We do not provide an implementation for text indexing, but have designed the implementation such that additional columns can be added to the move structure, such as the $\mathrm{BWT}$ character for $\mathrm{LF}$. These can either be bitpacked alongside move structure columns for cache efficiency or stored separated in their own dedicated memory. This makes Orbit well suited to be used in various scenarios as a flexible data structure library to construct larger and more intricate applications requiring a move structure.

To evaluate the efficiency of length capping and its implementation in practice, we evaluate it against approaches which use balancing for the $\mathrm{RLBWT}$ based permutations $\mathrm{LF}$ and $\varphi$. To our knowledge, the only known implementation of a generic move structure library supporting balancing is Move-r [1]. However, their approach performs some fixed length splitting that cannot be turned off. For that reason, we also use the r-permute tool¹¹1https://github.com/drnatebrown/r-permute to obtain balanced intervals for $\mathrm{LF}$ without any length splitting. Unfortunately, it does not support arbitrary balancing and thus we cannot do the same for $\varphi$.

To generate the permutations for $\mathrm{LF}$ and $\varphi$, we computed the $\mathrm{RLBWT}$ of $16,32,64,128$, $256$ and $1000$ concatenated haplotypes of human chromosome-19 (chr19) used in the experiments of Boucher et al. [2]. The specifics of these collections are described in Table 1. For $\mathrm{LF}$, we perform the move queries involved in $\mathrm{BWT}$ inversion, and use relative positions (i.e., using $S_{\mathrm{\ell }}$). For $\varphi$, we enumerate the suffix array, which requires absolute positions (i.e., using $S$). Both of these processes involve exactly $n$ distinct move queries and thus benefit from length capping. Construction and query time were measured using GNU time on a server with an Intel(R) Xeon(R) Gold 6248R CPU running at 3.00 GHz with 48 cores and 1.5TB DDR4 memory. We report build times and peak memory usage for Orbit across length capping parameters and Move-r across balancing parameters; we do not measure for other modes since they involve interfacing between implementations, whereas these modes are encapsulated by the listed tool.

We use $16$ copies of chr19 to explore the effect of various parameters for length capping and balancing. Timings in this section are averaged over 5 iterations. To evaluate $\mathrm{LF}$, we compared metrics using one thread for construction and query using:

• $\mathrm{■}$

  Orbit on an unbalanced move structure.

• $\mathrm{■}$

  Orbit balanced using $\alpha \in \{2,4,8,16\}$ with r-permute. Values larger than $16$ are identical to an unbalanced move structure.

• $\mathrm{■}$

  Orbit length capped using $c\in \{0.5,1,2,4,8,16,32\}$.

• $\mathrm{■}$

  Move-r balanced using $\alpha \in \{2,4,8,16\}$.We note that Move-r length splits first through binning into preset bitwidths.

• $\mathrm{■}$

  Orbit using pairs $(c,\alpha )\in \{(1,2),(2,4),(4,8),\mathrm{\dots },(32,64)\}$ by length capping first with Orbit, then passing this result to Move-r to balance the length capped intervals. This result is then loaded back into Orbit.

• $\mathrm{■}$

  Move-r using the same procedure as above using $(c,\alpha )$ pairs, but this time loaded into Move-r.

Note that we chose to use Move-r balancing instead of r-permute for the combined case of length capping and balancing so that we could use the same underlying move structure. This allows a direct comparison of Orbit and Move-r to evaluate the benefit of bitpacking when using both optimizations.

The results of the experiment are shown in Table 2. We observe that balancing alone has little effect compared to an unbalanced move structure; this result has been observed in other studies for $\mathrm{LF}$, which seems to have good practical query time even unbalanced [5]. However, length capping results in the fastest and smallest approaches, achieving a space reduction of $\sim 46\%$ compared to the unbalanced version. When directly compared to Move-r when using both balancing and length capping, Orbit is faster and smaller, although speed improvements are slight.

Although it is obvious from our theory why length capping results should be smaller, speed improvements in length capping may also be due to the smaller representation causing more rows of the move structure to fit in cache. Thus, our experiments are too limited to imply that Orbit itself is a faster implementation than Move-r, but does imply the benefit of length capping itself empirically. Construction times are improved over balancing, showing the theoretical improvement gained by length capping. Construction memory is also smaller for reasonable choices of $c$ and $\alpha$, due to the simplicity of splitting by length when compared with balancing. Overall, the best query times combine both length capping and balancing.

The equivalent experiment using the $\varphi$ permutation was performed with slight changes. First, we cannot obtain a balanced move structure without the implicit length splitting of Move-r, since this was only available for $\mathrm{LF}$. Second, the $\varphi$ permutation permits higher $\alpha$ values before stagnation than $\mathrm{LF}$. The results of the experiment are shown in Table 3. We see that unbalanced move structure queries for $\varphi$ performs $\sim 10$ times slower or worse when compared to splitting approaches. Further, the best query approaches involve some form of balancing, with the very best approaches being a combination of length capping and balancing. This implies that, in contrast to $\mathrm{LF}$, in practice $\varphi$ does not have good query time when left unbalanced, meaning even just the length capping guarantees result in a large performance improvement.

When directly compared to Move-r when balancing and length capping, query speed improvements of Orbit can be quite small, seemingly correlating with a slight space improvement. However, both the fastest and smallest approaches for queries use Orbit formulations. Note that only Orbit can make proper use of length capping, since its data structures bitpack to fit. Thus, the smaller space improvement could be attributed to requiring absolute positions rather than relative, which limits the benefit of bounding interval lengths. Regardless, length capping is still a clear query speed improvement when compared to an unbalanced move structure. We also again see that construction times for length capping are faster than balancing, with memory also smaller for reasonable choice of $c$ and $\alpha$.

To evaluate splitting approaches for larger collections and as $n/r$ grows, we perform $\mathrm{LF}$ inversion move query steps using Orbit. Overall, we compared an unbalanced move structure, length capping with $c=4$, balancing with $\alpha =8$, and the combination with $c=4$ and $d=8$. These parameters were selected due to their strong performance in both experiments of Section 5.3.

Results for splitting at scale are shown in Figure 2. Methods which perform length capping are consistently faster than those that do not, with only a slight improvement from balancing when compared to unbalanced. The time per query for length capping and the combination of splitting approaches is comparable with no decisive best approach. With respect to space, length capping approaches are at least $\sim 40\%$ smaller than the unbalanced or only balanced move structures. This trend is consistent across varying collections of chr19 sequences.