后缀BWT与循环移位BWT及快速计算

# 后缀 BWT 与循环移位 BWT，以及快速计算 来源：https://purplesyringa.moe/blog/suffix-bwt-vs-cyclic-shift-bwt-and-fast-computation/ 2026年7月3日 > 目标读者：数据压缩爱好者。 Burrows-Wheeler 变换（https://en.wikipedia.org/wiki/Burrows%E2%80%93Wheeler_transform）接收一个字符串作为输入，并对其字符进行重排，使它们按上下文分组。该变换是可逆的，只需额外输入 O(1) 数据，这些特性使其在数据压缩和基因组比对中占有一席之地。维基百科没有告诉你的是，实际上 BWT 有两种变体，它们在性能特性和简单性上有着微妙的差异。这些差异似乎几乎没有被文档记录，因此这篇博文是我试图理清它们的一次尝试。 ### 循环移位 BWT 我先介绍我称为“循环移位 BWT”的变体。以字符串 `bcacaba` 为例。我们写下该字符串的所有循环移位（或旋转），以及它们的偏移量： ``` bcacaba (0) cacabab (1) acababc (2) cababca (3) ababcac (4) babcaca (5) abcacab (6) ``` 然后按字典序排序： ``` ababcac (4) abcacab (6) acababc (2) babcaca (5) bcacaba (0) cababca (3) cacabab (1) ``` BWT 的输出是每个循环移位开始位置之前的那个字符。例如，列表中第一个循环移位从索引 4 开始，所以我们写下输入中偏移量 3 处的字符，即 `c`。最终输出为 `cbcaaab`。（注意这等同于取上面表格的最后一列。） 排在一起的两个字符串很可能有较长的公共前缀，因此在实际数据中，它们之前的符号也往往匹配。例如，在给定文本中，如果字符串 `ender` 前面总是跟着 `g`、`b` 或空格，那么输出中就会有一个完全由这三个符号组成的长子序列。 该算法是可逆的：如果你知道 BWT 的输出，你可以恢复输入（直到一个循环移位）。具体做法如下： - 准备工作：取 BWT 的输出（`cbcaaab`）并将其字符按字典序排序。这样就得到了上面表格的第一列：`aaabbcc`。 - 现在取输出的第一个字符 `c`。它是由某个循环移位 t‖"c" 的最后一个字符产生的。我们将寻找这个 `c` 之前的字符，即 t 的最后一个字符。 - 文本中有两个 `c`，所以有两个以 `c` 结尾的循环移位：t1‖"c" 和 t2‖"c"。由于我们取的是输出的第一个 `c`，而不是第二个，它一定属于字典序较小的那个移位 t‖"c" = t1‖"c" < t2‖"c"。 - 但我们还不知道哪个移位更大。不过，这并不重要，因为我们现在正在恢复字符串，而循环移位 BWT 的逆变换只能恢复到一个循环移位。因此，我们随意选取一个，比如 t1。由于 `c` 是输出中的第一个 `c`，它对应第一列中的第一个 `c`（在 `aaabbcc` 中位置为 5）。因此，该行的第一列是 `c`，它是导致该行成为字典序最小的候选者的原因：t1‖"c" 以 `c` 开头，而第二行的第一个字符 `c` 并不能与 t2 的 `b` 竞争。 - 让我们用同样的方式处理 `b`、`c` 等等。完整算法称为“last-to-first mapping”（LF mapping），因为你将输出中特定位置的字符映射到第一列中对应字符的第 k 次出现。最终我们得到 `bcacaba`——一个循环移位。 ### 后缀 BWT 另一种变体称为“后缀 BWT”。它采用字符串的所有后缀，而不是循环移位。但后缀没有最后一个字符——除了空后缀（empty suffix）。所以我们需要一个技巧来实现“环绕”语义。对于字符串 `bcacaba`，我们显式地将空后缀放在最后，或者假定它在最前面（取决于约定）。对后缀进行排序，并输出每个后缀前面的字符： ``` [a] bcacaba [ ] (empty suffix) ``` 后缀 `bcacaba` 整体上相对于 `cacaba`、`acaba` 等进行比较。但我们也需要比较后缀的相对顺序，以确定 BWT 输出顺序。假设我们按字典序排序： ``` (empty suffix) a aba acaba ba bcacaba caba cacaba ``` 现在 BWT 输出是每个后缀前面一个字符： ``` (empty suffix): 前面没有字符——假设为某默认字符，比如 '' a: 前面是 b aba: 前面是 c acaba: 前面是 c ba: 前面是 a bcacaba: 前面是 ''? 又是边界 caba: 前面是 a cacaba: 前面是 b ``` 所以我们得到了 `bccaaab`，以及一个 undefined。边界问题可以避免：某些实现约定空后缀放在开头，或者使用终止符 `$`。但无论如何，边界处始终是个问题。 在经典实现中，你通常会看到这样的解码算法： ```c char decode(const char *bwt, size_t primary_index) { // 第一步：计算字符计数 // 第二步：LF mapping，但这次跳过该 primary_index 而不是空后缀 } ``` 循环移位 BWT 的主索引（primary index）就是给定循环移位在排序列表中的索引。在后缀 BWT 中，它是空后缀被放置的索引。如果使用空后缀（或 `$`），编解码器会跳过它。 ### 细微差别 我用来强调边界问题的例子是 `"a"` 是否出现在文本中。在 `"bcacaba"` 中，`"a"` 出现了，所以后缀列表包含 `"a"`。但是，在循环移位 BWT 中，循环移位 `"a"` 并不存在，因为 `"a"` 不是循环移位。所以循环移位 BWT 处理 `"a"` 与 `"aa"` 的方式可能不同。 让我们更仔细地比较这两种排序。对于后缀： ``` a # 前面是 b aba # 前面是 c acaba # 前面是 c ba # 前面是 a bcacaba # 前面是 ? (边界) caba # 前面是 a cacaba # 前面是 b ``` 对于循环移位（从偏移量 4 开始）： ``` ababcac # 前面是 c (来自偏移量 3) abcacab # 前面是 b acababc # 前面是 c babcaca # 前面是 a bcacaba # 前面是 ? (边界) cababca # 前面是 a cacabab # 前面是 b ``` 注意 `a` 作为后缀是有的，但作为循环移位却没有——它只作为 `ababcac` 的开头出现。这导致了 BWT 输出的不同。 如果我们天真地试图用后缀 BWT 来模拟循环移位 BWT，通过丢弃空后缀并假设最后一个字符（在 full suffix 前面）是某个值，就会遇到问题。假设我们在后缀列表的顶部插入空后缀，并说 full suffix 前面是某个字符（例如 `a`，因为它应该被环绕）。那么排序就变成了： ``` (empty): 前面是 ? -> 假设为 "a" a: 前面是 b aba: 前面是 c acaba: 前面是 c ba: 前面是 a bcacaba: 前面是 ? 边界，但这里不关心 caba: 前面是 a cacaba: 前面是 b ``` 我们想要的是：在字典序上，`"a"` 作为后缀应该小于空后缀？或者空后缀小于所有？循环移位的排序是 `"a"`（实际上不存在）应该介于空后缀和 `"aba"` 之间？但事实上 `"a"` 不是循环移位，所以这是不可能实现的。 这种不一致会导致解码失败。具体来说，当你用后缀 BWT 输出 `bccaaab` 并告诉解码器 full suffix 在索引 4 时，解码器会错误地将倒数第二个字符（即 `a`）视为 full suffix 前面的字符。但实际上那个位置对应的 suffix 是 `ba`，而不是 `a`。这是因为编码时我们丢失了 `"a"` 作为循环移位的信息。 ### 比较 鉴于这些细微差别，你可能会合理地认为每个人都使用循环移位 BWT，而后缀 BWT 是一个遗留物。但事实并非如此。虽然后缀 BWT 解码稍复杂且稍慢，但编码快得多。有多种快速方法用于排序后缀（https://en.wikipedia.org/wiki/Suffix_array），最实用的是线性时间的 SA-IS（https://zork.net/~st/jottings/sais.html）算法，但很少有用于循环移位的算法。新手第一个 O(n log n) 的后缀数组实现（https://cp-algorithms.com/string/suffix-array.html#on-log-n-approach）确实适用于循环移位，但我找不到文献中任何专注于循环移位的线性时间算法。所以大多数人这样做：导入 libsais（https://github.com/IlyaGrebnov/libsais），构建后缀数组，丢弃全后缀，插入空后缀，然后通过 `s[sa[i] - 1]` 构造后缀 BWT。或者直接调用 `libsais_bwt` 函数做同样的事情。在很多情况下，这完全没问题。你可以经常在字符串前加上一个比所有其他字符都小的 `^` 字符，这样全后缀就会在索引 0 处，你不需要调整 `i`。或者附加一个 `$` 字符，这样后缀 BWT 和循环移位 BWT 在 `$` 的位置上重合，效果几乎相同。 ### 欺骗 SA-IS 但如果你确实需要循环移位 BWT 呢？也许你的字母表是满的，实现真正的后缀 BWT 解码器太昂贵，或者你的问题明确要求循环移位 BWT。在这种情况下，我们仍然可以通过一些调整使用 SA-IS。 实现环绕语义的最简单方法是加倍字符串。字符串 s‖s 中起始偏移量 i, j < |s| 的两个后缀与 s 的对应循环移位比较方式完全相同。所以你可以将 s 加倍，计算后缀数组，丢弃索引 ≥ |s| 的元素，然后从该列表中计算循环 BWT 和主索引。 如果你不想付出双倍代价，还有另一种方法。取 s 的字典序最小的循环移位，记为 s'。对于 `s = "bcacaba"`，最小循环移位是 `s' = "ababcac"`。s' 的后缀与它的循环移位比较方式完全相同：`"c" < "cac"` 并且 `"cababca" < "cacabab"`。这是因为 s' 在截断点处起到了隐式 `$` 终止符的作用。由于循环 BWT 对输入旋转不敏感，计算 s' 的后缀 BWT 就能得到正确结果。我是从 Christoph Diegelmann 的 StackOverflow 答案（https://stackoverflow.com/a/31491630/5417677）中学到这个的。最小循环移位可以通过在 s‖s 上运行 Duval 算法（https://cp-algorithms.com/string/lyndon_factorization.html#finding-the-smallest-cyclic-shift）找到。它是线性时间的，但常数比 SA-IS 小，因此速度提升很大。唯一的问题是主索引会被旋转 s 破坏。有两种方法可以恢复它。第一，我们可以构建 s' 的后缀数组而不是直接跳转到 BWT，然后找到 s 在该列表中的位置。这可行，但成本较高。另一种方法是直接按定义计算主索引，即字典序小于 s 的 s 的循环移位数。Z-函数（https://cp-algorithms.com/string/z-function.html）可以线性时间地为每个后缀 (s‖s)[i...] 找到它与 s‖s 的公共前缀长度，这正是确定循环移位 i 是否小于 s 的位置。或者，根据数据的不同，使用二次算法可能更快。 这是我的实现（https://github.com/purplesyringa/computercraft-programs/blob/df3551a486e239f7761312ba8989560229ddd241/initrd-ng/initrd-ng/src/bz.rs#L4-L58）。 ### 结论 那么你应该使用哪一种？在大多数情况下，后缀 BWT 更好：高效的后缀 BWT 编码器实现更简单且拥有更高的顶级性能，而解码仅稍慢，且仅在字母表满时明显。我认为循环移位 BWT 在需要小巧解码器且不关心编码器大小时尤其出色。比如 Demo 场景（https://en.wikipedia.org/wiki/Demoscene）和一般的尺寸编码。在脚本语言中它也更重要，因为解码运行时不主要被内存访问主导，调整会明显减慢过程。 希望我没有在这篇文章中犯任何错误！这个话题我不太熟悉，所以如果我犯了错误，请随时给我发电子邮件。