最長公共子序列與最小編輯距離-你有更快的演算法麼?

haolujun發表於2017-12-20

前言

工作以後很少用到演算法,時間長了腦袋也不靈活。但是最近,我有幸在工作中又開始用到演算法,這使我欣喜若狂,這兩個演算法主要用在計算文字相似度的場景中。小夥伴看標題中的“最長公共子序列”和“最小編輯距離”演算法想當然都認為是爛大街的演算法了,有什麼好講的。當然那種樸素的$ O(N^{2}) $的演算法本文不去講,而是講更快速的演算法,而且還有變種問題。

最長公共子序列

對於字串P和T,求其最長公共子序列。

動態規劃演算法

\(dp_{ij} = \begin{cases} & dp_{i-1,j-1}+1 \text{ if p[i]=t[j] } \\\\ & max(dp_{i-1,j}, dp_{i, j-1}) \text{ 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)\)

最小編輯距離演算法

動態規劃演算法

\begin{cases}
& \text{ if p[i] = t[j]: } dp[i][j] = dp[i-1][j-1] \\
& \text{ 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%。嗯,不說了,我得去領年終獎去啦。。。。。

相關文章