Uzi Vishkin. Deterministic Sampling–A New Technique for Fast Pattern Matching (1991)

Abstract. 本文是 https://doi.org/10.1137/0220002 的阅读笔记。文中提出了一种确定性的采样方法用于加速字符串匹配，由此得到了操作次数 $O(n)$ 时间 $O(\log^* n)$ 的并行算法。

文中用到了大量的分治、分块、复杂度平衡、递归的技巧，虽然这些无非是思考并行算法的时候的重点考虑方向，但第一次见的时候感觉还是比较奇妙，整个论文读完像做了一道超大的构造题。

• 简介
• 前置知识
  • 字符串的周期性
  • CRCW PRAM 的选择
  • 前缀和问题
• 模式分析
• 基本文本分析
  • 基本常数时间算法
  • 基本最优加速算法
  • 模式串有周期性的情况
• 复杂度平衡后的 $O(\log^* n)$ 文本分析
• 进一步讨论
• 参考文献

简介

字符串匹配问题是一类常见的问题。我们重新叙述其严格定义和相关名词。

字符串匹配问题. 给定一个长度为 $n$ 的字符串 $T[1…n]$，称为文本串，和一个更短的长度为 $m$ 的字符串 $P[1…m]$，称为模式串。求所有的起始位置 $1\leq i\leq n-m+1$，满足 $P = T[i…i+m-1]$。

本文提出的是求解字符串匹配问题的并行算法，采用的计算模型是 CRCW PRAM，鉴于以前没有接触过，给出其定义：

定义 1.1（CRCW PRAM）. CRCW PRAM 系指 Concurrent-Reading-Concurrent-Writing Parallel Random Access Machine，即多台处理器公用一个无限大容量的共享存储器，公用一个同步时钟，允许并发的读写操作。

必须详细考察写操作在并行时如何被采纳。CRCW PRAM 进一步可以分为

• Common CRCW PRAM 只允许多台处理器同时写相同内容到内存单元中，否则机器行为未定义。
• Arbitary CRCW PRAM 任选一台作为优胜者写入。
• Priority CRCW PRAM 选取编号最小的作为优胜者写入。

在涉及到细节时我们再深入讨论计算模型相关的问题。

此问题的朴素算法是 $O(nm)$ 的，容易发现在 有 $O(nm)$ 个处理器的 CRCW PRAM 上可以做到 $O(1)$ 时间。非平凡的算法如 KMP，Boyer Moore 等通常包含了两个阶段：模式分析和文本分析。KMP 中的求 next 就属于模式分析，把文本串放在上面跑就属于文本分析。

本文的主要结果是

1. 一种新的字符串匹配问题的线性算法（串行）。
2. 一种新的 $O(\log^* n)$ 文本分析的并行算法，处理器数最优。
3. 一种需要先做模式分析的 $O(log^2 m / \log \log m)$ 文本分析算法，处理器数最优。
4. 一种随机高概率成功的 $O(\log m)$ 的模式分析算法，处理器数最优。
5. 一种非标准 PRAM 上的最优时间的文本分析算法，处理器数最优。

我们只试图看懂 1，2，3。

这里的 $\log^*$ 是一个没有见过的记号。这是一个增长极慢的函数：

定义 1.2. 我们定义 $\log^{(i)}$ 表示 $\log$ 函数的 $i$ 次复合。$\log^*(r)$ 定义为最小的 $i$ 使得 $\log^{(i)} \leq 2$。注意底数是 $2$。

这篇文章提出的核心科技 deterministic sampling 启发自我们如何“观察”一个模式串是否在文本串中出现：

• 采样. 首先找到一个模式串的局部。
• 搜索. 在主串中找所有的此局部出现的位置。
• 验证. 验证这个局部前后是否是模式串整体。

我们期望采到一个模式串的子集 $[s(1), …, s(l)]$，理想情况下其应当满足

1. 其长度尽可能短，这样我们的搜索量才可能不会太大。
2. 如果这个子集成功匹配，其前后一段距离都不可能产生合法的匹配，这样我们的验证量才可能不会太大。

这个采样并不能简单地用均匀随机构造，读者可以自行举出反例。但是却可以通过一个确定性的过程精细地构造，这就是所谓的 deterministic sample。

前置知识

字符串的周期性

定义 2.1.1. 令 $u, w$ 是两个字符串，称 $u$ 是 $w$ 的一个周期若 $w$ 是 $u^k$ 的一个前缀。

字符串 $w$ 的最短周期定义为最短的 $u$，$u$ 是 $w$ 的一个周期。若 $w$ 的最短周期长度不超过其本身长度的一半，称 $w$ 是周期性的。

收集一些我们熟知的 border 论。

命题 2.1.2（equivalence fact）. 上述定义等价于 $w$ 是 $uw$ 的前缀。

命题 2.1.3（conflicting occurrence fact）. 对于一个串 $w$，其最短周期是 $u$，若它出现在位置 $i$，则其不可能出现在位置 $j$，其中 $i < j < i + |u|$。

命题 2.1.4（gcd lemma / weak periodicity lemma）. 若 $w$ 有长度为 $p, q$ 的周期，且 $p + q \leq |w|$，那么 $w$ 也有长度为 $\gcd(p, q)$ 的周期。

命题 2.1.5 若 $v$ 出现在位置 $j$ 和 $j + P$，那么

1. $v$ 是周期性的，且其周期 $T$ 不超过 $P$。
2. $v$ 出现在 $j + T$。

下面是本文使用的一个关键定义。

定义 2.1.5（witness 数组）. 考虑一个非周期性的字符串 $P[1…m]$。对于任意的 $1 < i \leq m / 2$，都有 $m - i + 1$ 并非原串的 border。这表明，至少存在一个位置 $1\leq k\leq m - i + 1$ 使得 $P(k)\ne P(i + k - 1)$。$\mathrm{WITNESS}(i)$ 存储这样的一个 $k$。

此数组可以在 $O(\log m)$ [VI85] 甚至 $O(\log\log m)$ [BG88] 时间内求出（共 $O(m)$ 次操作），我们可能会找时间看比较简单的那一个。

CRCW PRAM 的选择

通过增加 $O(n)$ 个处理器，可以在 $O(1)$ 时间内用 Common CRCW PRAM 实现 Priority CRCW PRAM，增加后各种开销变化都是常数级别的。我们只需要考虑如下问题：

给定一个 $n$ 位长的 vector，求最左边的非零位。

原论文 reference 一篇八十年代打字机写的论文，不想看所以自己想一个。

解法.

首先考虑如何 $O(1)$ 时间 $O(m^2)$ 个处理器实现求长度为 $m$ 的 vector 的最小值。

第一步定义一个辅助数组 $g_{i, 0/1}$，全初始化为 $0$。给每个位分配 $m$ 个处理器，第 $i$ 位第 $j$ 个处理器用来比较 $v_i$ 是否小于等于 $v_j$。如果大于等于在 $g_{i, 1}$ 写 $1$，否则在 $g_{i, 0}$ 写 $1$。

第二步使用 $m$ 个处理器，如果仅 $g_{i, 1}$ 为一向返回值写 $v_i$。

然后我们只需要将数组分成 $\sqrt n$ 块，首先统计每块是否有数，然后用上述算法求出最小的有数的块，然后用上述算法求出有数的块中最左边的非零位。

本文中的所有操作可以在 Common CRCW PRAM 上完成。少数操作中最直观的算法需要使用 Priority CRCW PRAM，上述论证表明无额外开销地用 Common CRCW PRAM 实现。

其次，关于并行算法，我们有如下定理：

定理（Brent）. 任何共 $x$ 次基本操作、总时间 $t$ 的同步并行算法（理解为可以同时执行任意数量的计算，但是必须按照依赖关系定时）可以用 $p$ 个处理器在 $\lceil x / p\rceil + t$ 时间内实现。

证明. 假设 $i$ 时刻同时执行的指令有 $x_i$ 条，这 $i$ 条可以被乱序执行，所以用 $p$ 个处理器执行时间为

$$
\sum_{i=1}^t \lceil x_i / p\rceil \leq \sum_{i=1}^t x_i / p + 1 \leq \lceil x / p\rceil + t
$$

此外可以想象如果有 $n$ 次计算，要求时间为 $t$，那么需要的处理器数至少是 $n / t$ 个。若算法达到了此下界（允许乘常数），则称处理器数最优。我们对上面的问题的解法是处理器数最优的。

前缀和问题

前缀和问题. 给定数组 $a[1…n]$，求数组 $b[1…n]$ 使得 $b[i] = \sum_{j=1}^i a[j]$。

这个问题好像其实没有那么平凡。[CV89] 指出存在一个 $O(\log n / \log\log n)$ 的算法，我们有空再看。

容易发现下面的问题弱于上面的问题：

给定一个 $0/1$ 数组 $a[1…n]$，提取其中标记为 $1$ 的编号到一个较短的数组中。

本文提出的算法经常使用前缀和算法来处理这个问题。

模式分析

模式分析的首要目标是构造 deterministic sample（DS）。这个采样的构造需要辅助问题 column sample。步骤 2 计算 column sample，步骤 3 用 column sample 构造 deterministic sample。

不失一般性地假设 $m$ 是偶数。然后考虑模式串的 $m / 2$ 个复制，按照下面形式排列：

$$
\begin{matrix}
\text{Copy $c_{m / 2}$} & & & & - & \cdots & - & - & \cdots & - \\
\vdots & & & & \vdots & & \vdots & \vdots & & \\
\text{Copy $c_2$} & & - & \cdots & - & \cdots & - & - & & \\
\text{Copy $c_1$} & - & - & \cdots & - & \cdots & - & & & \\
\text{Column id} & 1 & 2 & \cdots & m / 2 & \cdots & m & m + 1 & \cdots & m + m / 2 - 1 \\
\end{matrix}
$$

Column sample. 选择至多 $\log m - 1$ 个列 $\overline{ds}(1), …, \overline{ds}(l)$（其中 $l < \log m$），附上一个字符 $car(\overline{ds}(i))$，满足

1. 恰好存在一个复制 $c_j$，满足
   1. $c_j$ 和所有选出的列都相交。即 $j\leq \overline{ds}(i) < j + m$。
   2. 相交处的字符恰好是附上的字符。即 $car(\overline{ds}(i)) = P(\overline{ds}(i) - j + 1)$。
2. 对于其他的赋值，至少有一列与其相交，并且交点处的字符与列所附的字符不同。即对于任意的 $c_j\ne c_k$ 都存在 $\overline{ds}(i)$ 满足 $k\leq \overline{ds}(i) < k + m$，且 $P(\overline{ds}(i) - k + 1)\ne car(\overline{ds}(i))$。

其解的存在性由步骤 2 的算法给出。假设满足上述条件 $1$ 的复制为 $c_x$。

Deterministic sample. 无非是将 $\overline{ds}$ 用 $c_x$ 中的位置来表达。具体地，$DS$ 是一个数组 $[ds(1), ds(2), …, ds(l)]$，其中 $l\leq \log m - 1$，且 $ds(j) = \overline{ds}(j) - x + 1$。

---

步骤 1. 判断模式串是否是周期性的。同时求 witness 数组。如果是周期性的，求出周期。

注意若原串是周期性的（周期为 $p$），那么其长度为 $2p - 1$ 的前缀必然是非周期性的，在第 4.3 节我们讲如何做周期性的情况。

此过程可以被 [VI85] 完成，没有必要卷到 [BG88] 的水平，因为大头不在这里。

---

步骤 2. 此步骤归纳地构造集合 $A_1, …, A_l$，其中 $A_0 = \{c_1, …, c_{m / 2}\}$，且 $|A_{i + 1}|\leq |A_i| / 2$，$|A_l| = 1$。

下文中，对于复制 $c_i, c_j$，若 $i < j$，称 $c_i$ 在 $c_j$ 左边，$c_j$ 在 $c_i$ 右边。

若 $A_i$ 中只有一个元素，那么此步骤结束，我们求得了 $l = i$，可以转向步骤 3。

否则考虑 $A_i$ 中包含的复制中最左的和最右的，由于其差距不超过 $m / 2$，witness 数组给出了一个列，这个列上包含了两个不同的字符 $c, c’$。那么和此列处字符等于 $c$ 的复制和等于 $c’$ 的复制中至少有一类是不超过一半的，将不超过一半的一类对应的字符作为 $\overline{ds}(i)$，将该字符置为 $car(\overline{ds}(i))$，将该类作为 $A_{i + 1}$。

例子. 考虑字符串 $\texttt{100111}$，其排列为

$$
\begin{matrix}
c_3 & & & 1 & 0 & 0 & 1 & 1 & 1 & & \\
c_2 & & 1 & 0 & 0 & 1 & 1 & 1 & & & \\
c_1 & 1 & 0 & 0 & 1 & 1 & 1 & & & & \\
\end{matrix}
$$

考察 $\text{WITNESS}(2)$，这个值可能是 $1$，这给出了一个列，$c_1$ 和 $c_3$ 在此列上不同：

$$
\begin{matrix}
c_3 & & & \color{blue}{1} & 0 & 0 & 1 & 1 & 1 & & \\
c_2 & & 1 & \color{blue}{0} & 0 & 1 & 1 & 1 & & & \\
c_1 & 1 & 0 & \color{blue}{0} & 1 & 1 & 1 & & & & \\
\end{matrix}
$$

这列上为 $\texttt{0}$ 的有 $c_1, c_2$，为 $\texttt{1}$ 的有 $c_3$，后者大小不超过一半，将 $\overline{ds}(1)$ 置为 $3$，$car(\overline{ds}(1))$ 置为 $\texttt{0}$，$A_2$ 置为 $\{c_3\}$，继续递归。

$A_2$ 的大小为 $1$，本步骤完成。

可以发现这一步骤的核心就是选出某一位等于 $c$ 或者 $c’$ 的元素，因此每一轮迭代时间复杂度是 $O(\log |A_i| / \log \log |A_i|)$，总操作量是 $O(|A_i|)$，处理器数最优。集合的大小是指数级下降的，于是总复杂度是 $O(\log^2 m / \log \log m)$，总操作量是 $O(m)$。

---

步骤 3. 求出 DS。

本步骤是平凡的。

---

整个模式分析流程中的复杂度大头是步骤 2，总复杂度为 $O(\log^2 m / \log\log m)$。

基本文本分析

本节给出两种文本分析的算法，一种是基本常数时间算法，需要 $O(1)$ 时间 $O(n\log m)$ 次操作；另一种是基本最优加速算法，需要 $O(\log m)$ 时间 $O(n)$ 次操作。这两种算法展示了两种极端的情况，文章最后会将其平衡一个 $O(\log^* n)$ 时间 $O(n)$ 操作次数的算法。

首先假设模式串是非周期性的，稍后再谈如何处理周期性串的情形。

我们对串分块，块长为 $m/2$。如果把每个块和其后面地 $m/2$ 个元素视作一组资料，那么不同的组之间完全没有相互影响。下面的两个算法都是将所有块并行处理的。

基本常数时间算法

步骤 1. 对于位置 $1\leq i \leq n - m + 1$，逐一检查是否有 $T(i - 1 + ds(j)) = P(ds(j))$ 均成立。若是，称 $i$ 为候选。

显然此步骤可以在 $O(1)$ 时间，$O(n\log m)$ 次操作内完成。

步骤 2. 只保留每一块的第一个和最后一个候选。

此步骤可以在 $O(1)$ 时间，$O(n)$ 次操作内完成。

步骤 3. 对于保留的候选，逐字节验证。

此步骤可以在 $O(1)$ 时间，$O(n)$ 次操作内完成（$n / (m / 2)\times m = O(n)$）

可能对于步骤 2，我们需要着重分析一下正确性。

命题（the ricochet property）. 假设模式分析中得到的 $A_l$ 中的唯一元素是 $c_x$，那么若 $i$ 是候选，则其前面 $x - 1$ 个元素不可能是答案，其后面 $m / 2 - x$ 个元素不可能是答案。

证明. 反证，发现与 column sample 的定义矛盾。

发现第一个候选和最后一个候选联合起来可以排除中间所有的候选，因此只需要考虑这两个候选。

---

容易发现本算法可以在总共 $O(n\log m)$ 次操作，用 $O(n\log m)$ 个处理器可以在常数时间内完成。

基本最优加速算法

现在我们被允许总时间 $O(\log n)$，但是总操作次数限制为 $O(n)$。那么思路其实很自然，我们将上述全并行的算法改为增量地执行即可。需要点明的是，为了使得总操作次数为 $O(n)$，每一轮的候选位置数量必须至少折半。

依旧对于每一个块分开考虑。初始时块中所有位置都是候选。我们顺序地考虑 $ds(1), ds(2), …, ds(n)$，每次能够排除一些元素，但是这远远不够。因此我们继续用上 ricochet property，我们发现核对完 $ds(1)$ 之后，和最左边的侯选位置的相对位置不包含在 $A_1$ 中的都不可能成为答案，而 $A_1$ 的大小自然能够将候选位置的数量给 bound 住。

理顺思路之后我们考虑并行实现。

步骤 1. 计算所有满足如下条件的位置 $i$：$T(i - 1 + ds(1)) = P(ds(1))$。

步骤 2-1. 找到所有位置中最左边和最右边的一个，记为 $T(a)$，$T(b)$。定义左辅助变量 $lg = a-1 + ds(1)$ 和右辅助变量 $rg = a - 1 + ds(1)$。

现在用 $A_1$ 来构造集合 $T_{lg}$ 和 $T_{rg}$，使得 $\overline{ds}(1)$ 分别被平移至与 $lg$ 和 $rg$ 重合，其他元素可以根据此偏移量算出。

命题（key correctness observation）. 不在 $T_{lg}\cup T_{rg}$ 中的元素不可能成为答案。

证明. 和 ricochet property 类似。

$T_{lg}\cup T_{rg}$ 中可能包含一些非候选的位置。将其中的候选位置称作真候选位置。

步骤 2-2. 用 $A_1$ 构造 $T_{lg}$ 和 $T_{rg}$，方法如上。

对于 $ds(\alpha), 2\leq \alpha \leq l$，容易将上述过程修改为：

步骤 2-1’. 计算 $T_{lg} \cup T_{rg}$ 中所有真候选位置中满足 $T(i - 1 + ds(\alpha))$ 最左边的和最右边的 $T(a), T(b)$，并重算 $lg, rg$。

步骤 2-2’. 用 $A_\alpha$ 更新 $T_{lg}, T_{rg}$。

这个过程需要仔细考虑处理器分配。为了做到最优处理器数，对于每一个块我们只有 $O(m/\log m)$ 个处理器。我们每一轮都只给 $T_{lg}$ 和 $T_{rg}$ 分配处理器，处理器将读取此位置是否为真候选，若不是则其将空转。然后发现前 $\log\log m$ 轮我们的处理器数是不够用的，但是这部分的总复杂度仍然是 $O(\log m)$ 的（因为每轮的处理时间指数级下降），后面 $\log m - \log \log m$ 轮处理器足够，每轮 $O(1)$，自然总复杂度是 $O(\log m)$，计算次数为 $O(m)$，处理器数最优。

步骤 3. 逐字符检查候选是否形成匹配。

综上所述，本算法复杂度为 $O(\log m)$，总操作数 $O(n)$，处理器数最优。

模式串有周期性的情况

不妨设模式串 $P = u^kv$，其中 $|u| = p$，$v$ 是 $u$ 的真前缀。

鉴于 $|v| < |u|$ 以及命题 2.1.3，算出 $u^k$ 的出现位置再验证后面 $|v|$ 位的复杂度正确，我们只需要先算出 $u^k$ 的出现位置。

首先用上述算法计算 $p’ = P[1, 2p - 1]$ 在 $T$ 中出现的位置。那么 $u^k$ 在某位置 $i$ 出现当且仅当 $p’$ 在 $i, i + p, …, i + (k - 2)p$ 均出现。上述条件 $\bmod p$ 同余的等价类上有一个连续的长度为 $k - 1$ 的区间。一切归结为如下问题：

哪些位置 $i$ 满足：从 $i$ 起，连续 $k - 1$ 位均为 1？

我们将串按长度 $k - 1$ 分块，则可以计算每块的第一个 $0$ 和最后一个 $0$（由此得到前缀 1 和后缀 1）的长度。现在对于位置为 $i$ 的 $1$，若其是后缀 1，且 $i + k - 2$ 打中了下一块的某个前缀 $1$，则其满足条件。于是上述过程可以在 $O(1)$ 时间、$O(n)$ 总操作内用最优处理器数完成。

复杂度平衡后的 $O(\log^* n)$ 文本分析

我们重新审视一下两个基本文本分析算法：

• 基本常数时间在一轮内并行验证所有的采样。
• 基本最优加速算法在 $O(\log n)$ 轮内逐步验证采样，每轮验证一个。

那么平衡的思路无非是，分成 $o(\log n)$ 轮逐步验证采样，每轮并行验证多个采样。需要注意的是，基本常数时间算法是带 $\log$ 的，所以需要预先像基本最优加速算法一样验证前几个元素。

我们的算法分为三个阶段：

阶段 1. 执行 $\log \log^* n + 1$ 轮基本最优加速算法。令 $\delta = \log\log^* n + 1$（现在已经处理了 $DS$ 的前几个元素）。

执行完上述操作后 $T_{lg}$ 和 $T_{rg}$ 中的元素个数上界为 $n / \log^* n$。

本阶段有 $O(n)$ 次操作，时间为 $O(\log \log^* n) < O(\log^* n)$。

阶段 2. 令 $count$ 从 $\log^*$ 递减到 $1$，每一轮用 $ds(\delta + 1), …, ds(\delta + \log^{(count)} n)$ 验证 $T_{lg}$ 和 $T_{rg}$ 中的真候选位置是否还是真候选。选出当前最左和最右的真候选，更新 $T_{lg}$ 和 $T_{rg}$。

每一轮我们会对每个候选位置执行 $O(\log^{(count)} n)$ 次验证操作，总时间为 $O(1)$，而候选位置的数量为 $n / 2^{\delta + \log^{(count)}n} = n / (\log^* n\log^{(count - 1)}n) \leq n / (\log^* n \log^{(count)} n)$，操作次数总共不超过 $O(n / \log^* n)$。所以阶段 2 的总时间为 $O(\log^* n)$，总操作次数为 $O(n)$。

阶段 3. 将候选位置前后文逐字符比较。

本阶段总操作次数 $O(n)$，时间为 $O(1)$。

综上，我们得到了总时间 $O(\log^* n)$，总操作次数为 $O(n)$，处理器数最优的文本分析算法。

进一步讨论

是否有亚线性匹配算法？

理论上很难。

是否能够向高维拓展？

虽然高维上没有类似于周期的概念，但是如果输入满足某种 conflict occurrence assumption 即任意将模式串摆在两个足够近的位置（每一维差距不超过 $m/2$）都可以发现重叠部分有一位不同（类似于 witness），做法就可以拓展到高维。总的来说难度很大。

是否能够处理近似匹配？

还是需要某种 conflict occurrence assumption，因此难度也很大。

参考文献

[VI85] Uzi Vishkin, Optimal parallel pattern matching in strings, Inform. and Control, 67 (1985), 91–113

[BG88] Dany Breslauer, Zvi Galil, An optimal $O(\log\log n)$ time parallel string matching algorithm, SIAM J. Comput., 19 (1990), 1051–1058

[CV89] Richard Cole, Uzi Vishkin, Faster optimal parallel prefix sums and list ranking, Inform. and Comput., 81 (1989), 334–352

[CV86] Cole, U. Vishkin, Deterministic coin tossing and accelerating cascades: micro and macro techniques for designing parallel algorithms, Proc. 18th Annual ACM Symposium on Theory of Computing, Association for Computing Machinery, New York, 1986, 206–219