全文检索看起来是"老问题"：给一段长文本建索引，回答"模式 P 出现了几次、在哪些位置、上下文是什么"。倒排索引（inverted index）在词粒度上无敌，但它要求先分词，对 DNA 序列、二进制字节流、列式存储的字典列这类"没有词边界"的数据无能为力。而把每个后缀都排序存下来（后缀数组）虽然能精确回答子串查询，却要 O(n) 空间且查询要二分 O(m·log n)。2000 年 Paolo Ferragina 与 Giovanni Manzini 提出的 FM-index，用一个巧妙的观察把后缀数组"压缩"进了它自己变换出的 BWT 字符串里——只存 BWT 加几个小辅助结构，就能在接近 n·Hk 的极低空间下，完成 count / locate / extract 三类查询，且完全随机访问。本文从空间困境讲起，逐层剖析 BWT、LF 映射、反向搜索与位级压缩，并给出一份可运行的 Python 教学版实现。
一、全文索引的空间困境：为什么不能把每个后缀都存下来
子串查询最朴素的做法是后缀数组：把所有 n 个后缀排字典序，存下每个后缀在原文的起始位置。查询模式 P（长度 m）时，在排序后的后缀数组上二分，找到"以 P 为前缀"的后缀区间 [sp, ep]，区间长度就是出现次数，区间里每个值就是出现位置。这没错，但代价是：
- 空间 O(n log n) 比特：每个后缀位置要用 log n 位，n 个就是 n·log n 位。对 3Gb 的人类基因组，这要几百 MB 仅存下标，还不算原文。
- 查询 O(m log n) 时间：每次二分都比较 m 个字符。
- 随机访问昂贵：要取回某次出现的上下文（extract），还得回原文切片，意味着原文得另外存着。
FM-index 的灵魂问题是：后缀数组里的"顺序信息"能不能用更小的结构表达，使得我们既不必存整个数组、又能还原出任意区间与位置？答案是把后缀数组喂进 Burrows-Wheeler Transform，再看变换后字符串的局部性。
二、后缀数组（Suffix Array）：O(n) 空间定位子串起点
先把后缀数组说清楚，因为它就是 FM-index 的"压缩对象"。给定文本 S（长度 n），其后缀数组 SA 是一个整数数组，SA[i] 表示字典序第 i 小的后缀在 S 中的起始下标。例如 S = "banana$"（'$' 是小于所有字母的哨兵，标记串尾）：
S = b a n a n a $ 后缀 起始 字典序 $ 6 最小 a$ 5 ana$ 3 anana$ 1 banana$ 0 na$ 4 nana$ 2 最大 SA = [6, 5, 3, 1, 0, 4, 2]
任意子串查询都可以转成对 SA 的区间查询：模式 P 出现的位置，恰好就是"以 P 为前缀"的那些后缀的起始下标，也就是 SA 里一段连续区间。难点只在于：这段区间怎么在不二分、不逐字符比较的情况下快速求得？FM-index 用 BWT 解决了它。
三、BWT 变换：把重复字符搬进相邻列
Burrows-Wheeler Transform 的构造非常直白：把文本的所有循环移位排成矩阵，按字典序排序，取每一行的最后一个字符拼起来，就是 BWT。对后缀数组视角，这等价于"取每个后缀的前一个字符"——因为 SA 已经把后缀排好序了。我们称排好序后的首列为 F、末列（BWT）为 L：
| 行 | 排序后的后缀 | F（首列） | L = BWT（上一行末字符） |
|---|---|---|---|
| 0 | $ | $ | a |
| 1 | a$ | a | n |
| 2 | ana$ | a | n |
| 3 | anana$ | a | b |
| 4 | banana$ | b | $ |
| 5 | na$ | n | a |
| 6 | nana$ | n | a |
BWT 的魔力在于局部性：原文中连续出现的相同字符，其"前一个字符"往往也相近，于是排好序后 L 列会出现长串同质字符（这里的 "nn"、两个 "a"）。这种同质块是后续压缩（游程、位图）能大显身手的根。注意 BWT 是可逆的：仅由 L 就能还原出原文（因为 F 就是 L 排序后的结果），这正是"压缩而不丢信息"的前提。
四、从 BWT 矩阵到 LF 映射：Last-to-First 重排
排序矩阵有个关键不变量：第 i 行的首字符，在第 i-1 行里是作为末字符出现的，且它在该字符的"第几个"保持不变。换句话说，把 L 中第 i 个位置的字符，映射回 F 中"同字符的第几个"，得到的新行号就叫 LF(i)（Last-to-First）。对 "banana$"：
L = "annb$aa" F = "$aaabnn"（F 就是 L 排序）
C['$']=0 C['a']=1 C['b']=4 C['n']=5 # C[c]：严格小于 c 的字符总数
rank('a', i) = L[0..i] 中 'a' 的出现次数
LF(i) = C[ L[i] ] + rank( L[i], i-1 )
LF[0]=C['a']+0=1
LF[1]=C['n']+0=5
LF[2]=C['n']+1=6
LF[3]=C['b']+0=4
LF[4]=C['$']+0=0
LF[5]=C['a']+1=2
LF[6]=C['a']+2=3
LF = [1, 5, 6, 4, 0, 2, 3]LF 映射的深层含义是 SA[LF[i]] = (SA[i] - 1) mod n：它把"后缀起点减 1"这个操作编码进了 BWT 的局部性里。也就是说，只要能算 LF，就能沿着后缀数组往"前一个字符"走，而完全不需要真的存 SA。这正是 FM-index 把 O(n log n) 位压缩到接近 n·Hk 的核心。
五、LF 映射与字符计数：用 rank 回答"前面有几个 c"
实现 LF 只需要两样东西：C[c]（字符频率前缀和，常数大小）和 rank(c, i)（BWT 前 i 个位置里 c 出现了几次）。rank 看似要扫描，实则可以用累计出现数组 O(1) 查到——这是 FM-index 唯一需要"额外存"的结构。工程上有两级优化：
- 原始 occ 数组：对每个字符维护一个长度为 n 的前缀计数数组，O(1) 查询、O(n·Σ) 空间。
- Wavelet Tree / 位图 + RRR：把每个字符的出现位图压成稀疏位图，用 select/rank 在 O(log Σ) 时间内回答，空间降到接近 n·H0（零阶熵）。这就是"压缩"二字落地的位置。
只要 rank 够快够小，LF 就是 O(1) 操作，整个索引的查询骨架就立住了。
六、反向搜索（Backward Search）：从模式末尾倒着走
怎么不求区间地直接算出"以 P 为前缀的后缀区间"？技巧是从模式的最后一个字符往回处理。每加入一个字符 c（从后往前），区间 [sp, ep] 被收窄到"在 BWT 的 [sp, ep] 这一段里，字符 c 所占据的那段新区间"。公式只用 C 和 rank：
sp, ep = 0, n-1
for c in reversed(P): # 从模式末尾往前
lo = rank(c, sp-1) # BWT[0..sp-1] 中 c 的个数
hi = rank(c, ep) # BWT[0..ep] 中 c 的个数
if hi <= lo: return 空 # 区间内没有 c，模式不存在
sp = C[c] + lo
ep = C[c] + hi - 1
# 结束后 [sp, ep] 即"以 P 为前缀"的后缀在 SA 中的区间用 "ana" 在 "banana$" 上走一遍：起始 sp=0,ep=6。先处理 'a'：lo=0, hi=3 → sp=1,ep=3（所有以 'a' 开头的后缀）；再处理 'n'：在 [1,3] 内 'n' 的 lo=0,hi=2 → sp=5,ep=6（"na$"、"nana$"）；最后处理 'a'：在 [5,6] 内 'a' 的 lo=1,hi=3 → sp=2,ep=3。落点区间 [2,3] 正好对应后缀 "ana$" 与 "anana$"，出现次数 2。整个过程没有比较任何原文字符，只做了 3 次 O(1) 的 rank/LF 计算。
七、三种核心查询：count / locate / extract
反向搜索结束得到的区间 [sp, ep] 直接给出三大查询：
- count：区间长度 ep - sp + 1，即出现次数。O(m) 时间，且不需要回原文。
- locate：区间内每个行号 i 对应的后缀起点 SA[i]。由不变量
SA[LF[i]] = (SA[i] - 1) mod n，从行 i 沿 LF 一直走到 '$' 所在的首行（行 0），步数 k 满足SA[i] = (SA[0] + k) mod n。SA[0] 是哨兵位置（恒为 n-1），所以每定位一个出现只需 O(n) 步——生产系统会对 SA 做稀疏采样（每隔 √n 或 log n 行存一个真值）把这一步压到 O(log n)。 - extract：拿到某个出现的起点后，取回它的上下文子串。教学版直接回原文切片；压缩版则沿 LF 向前走 length 步逐字符拼回（注意 LF 给出的是"前一个字符"，方向要反转）。
def count(self, pattern):
sp, ep = self.backward_search(pattern)
return max(0, ep - sp + 1)
def _sa_of_row(self, i):
steps, pos = 0, i
while pos != 0: # 走回 '$' 所在的首行
pos = self.lf(pos)
steps += 1
return (self.sa[0] + steps) % self.n # SA[0] 即哨兵位置
def locate(self, pattern):
sp, ep = self.backward_search(pattern)
return sorted(self._sa_of_row(i) for i in range(sp, ep + 1))
def extract(self, sa_pos, length):
return self.text[sa_pos: sa_pos + length] # 教学版回原文切片八、压缩骨架：Wavelet Tree 与 occ 的位级压缩
FM-index 之所以能"压缩"，关键在于 rank 用的 occ 数组本身可以被压。最经典的落法是 Wavelet Tree：把字母表递归二分，每一层用一张位图标记"落到左半还是右半"，rank 退化成在位图上数 1 的个数。位图再经 RRR（Raman–Raman–Rao）或 SDArray 编码压到接近零阶熵。这样整体空间：
# Wavelet Tree 上的 rank 示意（两层示例）
# 字母表 {a,b,n,$} 二分：左={$,a} 右={b,n}
bitmap = "0100110" # 1 表示落在右半
# rank(c, i)：先按位图定位在右半的第几个，再递归到对应子层
def wt_rank(c, i):
if 叶子: return occ[c][i]
bit = 1 if c in 右半 else 0
k = rank1(bitmap, i) if bit else (i - rank1(bitmap, i))
return wt_rank(子层字符(c), k)配合对 SA 的稀疏采样（存 √n 个真值）与对原文可选的 Golomb/Rice 编码，FM-index 在实测中能把英文文本压到 约 1.2–2 字节/字符，同时保留全部子串查询能力——这比存完整后缀数组（约 4 字节/字符仅下标）加上原文要省得多。
九、空间账：为何能压到接近 1.2–2 字节/字符
把几类全文索引方案摊开对比，FM-index 在"子串随机访问 + 高压缩"这一格是唯一的：
| 方案 | 空间（字节/字符） | 子串查询 | 随机 extract | 典型用途 |
|---|---|---|---|---|
| 后缀数组（裸） | ~4（仅下标） | O(m log n) | 需回原文 | 教学、小规模 |
| 后缀树 | ~8–20（指针密集） | O(m) | 需回原文 | 多模式、LLM |
| 倒排索引 | 依赖分词 | 词级 O(1) | 不支持子串 | 搜索引擎 |
| FM-index（压缩） | ~1.2–2 | O(m) | 支持 | 基因组、字节流 |
注意空间随文本熵走而非随 n 走：文本越有规律（DNA、日志、列式字典列），压缩比越高。这是它天然适合"高重复、无词界"数据的根本原因。
十、工程落地：BWA、Bowtie2 与列式存储的字典压缩
FM-index 不是论文玩具，它是现代生命科学与数据系统的地基：
- 基因组比对：BWA、Bowtie2 用2 阶 FM-index（2FM-index）在数十 Gbps 的人类中做短读长比对，靠的就是 count/locate 在压缩空间上的 O(m) 查询。人类基因组 3Gb 的索引约 3–4GB，且能随机访问任意子串。
- 字符串去重与模糊匹配：rsync 风格的 rolling hash、模糊搜索（近似 count）都建立在反向搜索的区间收窄上。
- 列式存储字典压缩：ClickHouse、Parquet 的字典编码列在做"前缀/子串谓词"时，可直接对字典列建 FM-index，把
LIKE '%abc%'变成一次 count 而非全列扫描。 - 后缀自动机/后缀树退化场景：当内存紧张又必须保留子串查询时，FM-index 是后缀树的高压缩替代。
十一、与后缀数组 / 后缀树的取舍
FM-index 不是"更好"，而是"在不同约束下更合适"：
- 要最小空间 + 能忍 O(m) 查询 → FM-index。它把后缀数组"隐含"在 BWT 里，用计算换存储。
- 要 O(m) 查询且频繁做多种模式/子树操作（如后缀链接） → 后缀树/后缀自动机更直观，但指针密集、空间大。
- 要最快的 locate → 裸后缀数组 + 原文，二分 O(m log n) 但实现最简单，适合内存宽裕的小数据。
- 要词级检索 → 倒排索引仍是王者，FM-index 不擅长"词"语义。
实战里常见组合：用 FM-index 做底层压缩子串索引，上面再铺一层倒排或采样，兼顾空间与查询语义。
十二、结语
FM-index 的优雅在于它把"排好序的后缀数组"这一看似必须显式存储的结构，重新编码进文本自己的 BWT 变换之中：BWT 带来局部性、LF 映射把"前缀缩减"变成 O(1) 的 rank 查询、反向搜索让 count 完全不碰原文、稀疏采样的 SA 与 Wavelet Tree 把 locate/extract 也压进接近熵的空间。它证明了一件事——很多"必须存下来"的顺序信息，其实可以靠可逆变换 + 局部性 + 代数映射隐式表达。这也是为什么从人类基因组到列式存储的字典列，FM-index 始终是"高压缩、可随机访问子串索引"这一需求的第一选择。下面的教学实现把上述全部机制浓缩进约 40 行 Python，可直接跑通 count / locate / extract。
class FMIndex:
def __init__(self, text):
self.text = text
self.n = len(text)
self.sa = sorted(range(self.n), key=lambda i: text[i:])
self.bwt = ''.join(text[(self.sa[i] - 1) % self.n] for i in range(self.n))
from collections import Counter
cnt = Counter(text)
self.chars = sorted(cnt)
self.C = {}
tot = 0
for c in self.chars:
self.C[c] = tot
tot += cnt[c]
self.occ = {c: [0] * (self.n + 1) for c in self.chars}
for i, ch in enumerate(self.bwt):
for c in self.chars:
self.occ[c][i + 1] = self.occ[c][i]
self.occ[ch][i + 1] += 1
def rank(self, ch, i):
return self.occ[ch][i + 1] if i >= 0 else 0
def lf(self, i):
return self.C[self.bwt[i]] + self.rank(self.bwt[i], i - 1)
def backward_search(self, pattern):
sp, ep = 0, self.n - 1
for ch in reversed(pattern):
if ch not in self.C: return (0, -1)
lo = self.rank(ch, sp - 1)
hi = self.rank(ch, ep)
if hi <= lo: return (0, -1)
sp = self.C[ch] + lo
ep = self.C[ch] + hi - 1
return (sp, ep)
def count(self, pattern):
sp, ep = self.backward_search(pattern)
return max(0, ep - sp + 1)
def _sa_of_row(self, i):
steps, pos = 0, i
while pos != 0:
pos = self.lf(pos)
steps += 1
return (self.sa[0] + steps) % self.n
def locate(self, pattern):
sp, ep = self.backward_search(pattern)
return sorted(self._sa_of_row(i) for i in range(sp, ep + 1))
def extract(self, sa_pos, length):
return self.text[sa_pos: sa_pos + length]
# 端到端示例
idx = FMIndex("banana$")
print(idx.count("ana")) # 2
print(idx.locate("ana")) # [1, 3]
print(idx.extract(1, 3)) # "ana"
print(idx.count("na")) # 2
print(idx.locate("na")) # [2, 4]
发表评论 取消回复