最长公共子序列与最小编辑距离-你有更快的算法么？

前言

工作以后很少用到算法，时间长了脑袋也不灵活。但是最近，我有幸在工作中又开始用到算法，这使我欣喜若狂，这两个算法主要用在计算文本相似度的场景中。小伙伴看标题中的“最长公共子序列”和“最小编辑距离”算法想当然都认为是烂大街的算法了，有什么好讲的。当然那种朴素的$ O(N^{2}) $的算法本文不去讲，而是讲更快速的算法，而且还有变种问题。

最长公共子序列

对于字符串P和T，求其最长公共子序列。

动态规划算法

(dp_{ij} = egin{cases} & dp_{i-1,j-1}+1 ext{ if p[i]=t[j] } \\ & max(dp_{i-1,j}, dp_{i, j-1}) ext{ if p[i]!=t[j] } end{cases})

这个动规方程不讲了，大伙都知道原理。

转换成最长单调递增子序列

P="abdba"
T="dbaaba"

• 步骤1：对于P中的每个字符，找到其在T中的位置。

pos(a) = {3,4,6}
pos(b) = {2,5}
pos(d) = {1}

• 步骤2：把P用字母在T中的位置替换，得到新串C，注意位置要由大到小排列。

C = [<6,4,3>,<5,2>,<1>,<5,2>,<6,4,3>]

• 步骤3：求C的严格的单调递增子序列，子序列长度即为最长公共子序列长度。

本例严格最长单调递增序列为:

4,5,6  => a,b,a   
2,5,6  => b,b,a  
3,5,6  => a,b,a  

每个严格最长单调递增序列都对应一个最长公共子序列，这个算法的平均时间复杂度为$ O(nlogn) $。

变种最长公共子序列

我在工作中需要解决的问题不是这么简单的最长公共子序列问题，而是有些变化。重新看这个问题，可以把最长公共子序列看成是两个字符串的匹配，每个字符如果匹配上则得1分，我们本质是求这个最大匹配。那么如果每个字符匹配得分不是1分而是任意分将如何做呢？这里描述一下这个变种问题：

P: 模式串
a b c d c d  <- 这行是字符
1 1 1 2 2 3 <- 这行是权重
T: 文本串
a c b d  <- 这行是字符
1 2 1 3  <- 这行是权重

对于模式串和文本串的每个字符，在其下方都标有权重，匹配规则如下：

• 字符不匹配得0分
• 字符匹配但是权重不匹配得1分
• 字符与权重都匹配得权重分数

可以利用动态规划算法：
(dp_{ij} = max(dp_{i-1,j},dp_{i,j-1},dp_{i-1,j-1}+score(p_{i},t_{j})))，时间复杂度(O(N^2))，求得结果为6。

如何利用最长公共子序列算法呢？我们发现如上的转换，只适合匹配权重为1的情况，如果想适应到各种灵活的权重，需要做一个拆分和一个转换。
拆分规则：

• 对于P中的每个字符，把其重复W次，W为其权重。

P: 模式串
a b c d c d => a b c d d c c d d d
1 1 1 2 2 3    1 1 1 1 1 1 1 1 1 1 <= 这行是权重，拆分后都为1
0 1 2 3 4 5    0 1 2 3 4 5 6 7 8 9 <= 这行是下标

序列转换规则：

• 如果T中字符与P中字符相同但是权重不同则下标重复1次。
• 如果T中字符与P中字符相同权重也相同则下标重复W次，W为权重。
• 转换时，T从前往后，P从后向前匹配，重复下标为拆分后下标。

此时得到的序列是：
[<0> <6,5> <6,5> <2> <1> <9,8,7> <9,8,7> <9,8,7> <4,3>] 

最长严格上升子序列为:
0 5 6 7 8 9 => a c d
               1 2 3 <= 这行为权重

在权重较小时（我遇到的实际问题中，w<=3），算法的平均时间复杂度还是仍然为(O(nlogn))。

最小编辑距离算法

动态规划算法

egin{cases}
& ext{ if p[i] = t[j]: } dp[i][j] = dp[i-1][j-1] \
& ext{ if p[i] !=t[j]: } dp[i][j] = min(dp[i][j-1],dp[i-1][j],dp[i-1][j-1]) + 1
end{cases}

最普通烂大街算法，时间复杂度(O(mn))。

列划分算法

参考文献：《Theoretical and empirical comparisons of approximate string matching algorithms》。In Proceedings of the 3rd Annual Symposium on Combinatorial Pattern Matching, number 664 in Lecture Notes in Computer Science, pages 175~184. Springer-Verlag, 1992。Author:WI Chang ，J Lampe，时间复杂度为(O(mn/sqrt[2]{b}))，其中b为字符集大小。这篇论文足足看了一个星期，终于看明白了其中的原理，而且这个算法基本上是解决最小编辑距离问题的最快速算法。观察编辑矩阵发现，矩阵每一列都包含多个子序列，子序列内编辑距离连续，子序列之间编辑距离相差1。

看图发现，编辑矩阵可以分割成黄蓝间隔的子序列，我们只要能够计算每个序列的起始位置以及长度，就可以计算最小编辑距离，而序列内编辑距离连续可以直接用加法跳过。
在这里直接给出作者的结论。

• 定义1：对于编辑矩阵第j列上的某个连续子序列，当i-dp[i][j]=h，我们就说dp[i][j]属于子序列h。
• 定义2：子序列h截止位置为i，当且仅当dp[i+1][j]属于另一个子序列t，t>h。
• 定义3：子序列h截止位置为i，长度为0，当且仅当dp[i+1][j] < dp[i][j]。
• 推论1：如果编辑矩阵列i上，子序列h长度等于0，且截止位置为j，则列i+1上的子序列h截止位置在j+1。
• 推论2：如果编辑距离矩阵列i上，子序列h长度L大于0，且截止位置为j，若在q=[j-L+2,j+1]上，存在最小的坐标r，p[r] = t[i+1]，则列i+1上的子序列h截止位置为r-1；如果没有这样的r并且列i上的子序列h+1长度大于0，则列i+1上的子序列h的截止位置为j+1；除此之外列i+1上的子序列h截止位置为j。
• 推论3：没有连续的两个子序列长度都是0。

这样我们就可以只追踪列划分点而不用追踪编辑矩阵所有元素来计算编辑距离。下面给出一个示例代码，比较丑，能看即可：

inline int loc(int find[][200], int *len, int ch, int pos) {
  for(int i = 0; i < len[ch]; ++i) {
    if(find[ch][i] >= pos)  return find[ch][i];
  }
  return -1;
}

int new_column_partition(char *p, char *t) {
  int len_p = strlen(p);
  int len_t = strlen(t);
  int find[26][200];
  int len[26] = {0};
  int part[200];  
  //生成loc表，用来快速查询
  for(int i = 0; i < len_p; ++i) {
    find[p[i] - 'a'][len[p[i] - 'a']++] = i + 1;
  }
  
  int pre_cn = 0, next_cn = 1, min_v = len_p;
  part[0] = len_p;
  
  for(int i = 0; i < len_t; ++i) {
    //前一列partition数
    pre_cn = next_cn;
    next_cn = 0;
    int l = part[0] + 1;
    int b = 1;
    int e = l;
    int tmp;
    int tmp_value = 0;
    int pre_v = part[0];
    //前一列第0个partition长度肯定>=1
    if(len[t[i] - 'a'] >0 && (tmp = loc(find, len, t[i] - 'a', b)) != -1 && tmp <= e) {
      part[next_cn++] = tmp - 1;
    } else if(pre_cn >= 2 && part[1] - part[0] != 0){
      part[next_cn++] = part[0] + 1;
    } else {
      part[next_cn++] = part[0];
    }
    //每列第一个partition尾值
    tmp_value = part[0];

    //遍历前一列剩下的partition
    for(int j = 1; j < pre_cn && part[next_cn - 1] < len_p; ++j) {
      int x = part[j], y = pre_v;
      pre_v = part[j];
      l = x - y;
      if(l == 0) {
        part[next_cn++] = x + 1;
      } else {
        b = x - l + 2;
        e = x + 1;
        if(b <= len_p && len[t[i] - 'a'] > 0 && (tmp = loc(find, len, t[i] - 'a', b)) != -1 && tmp <= e) {
          part[next_cn++] = tmp - 1;
        } else if(j + 1 < pre_cn && part[j + 1] - x != 0) {
          part[next_cn++] = x + 1;
        } else {
          part[next_cn++] = x;
        }
      }
      l = part[j] - part[j - 1];
      if(l == 0) {
        //新得到的partition长度为0，那么下一个partition的起始值比上一个partition尾值少1
        tmp_value -= 1;
      } else {
        tmp_value += l - 1;
      }
    }
    
    if(part[next_cn - 1] != len_p) {
      part[next_cn++] = len_p;  
      tmp_value += len_p - part[next_cn - 2] - 1;
      if(tmp_value < min_v) {
        min_v = tmp_value;
      }
    } else {
      min_v = min_v < tmp_value ? min_v : tmp_value;
    }
  }
  return min_v;
}

总结

这篇博客中的两个算法算是我今年重点研究的算法了，带来的收益很不错，直接把检索的某个模块性能提高50%，再加上我的这篇博客《我是怎么用跳表优化搜索引擎的？》 中写到的优化，使得检索系统的整个成本降低50%。嗯，不说了，我得去领年终奖去啦。。。。。