最大の全体的なサブシーケンスを見つける

この記事では、最大共通サブシーケンスまたはLCS（最長共通シーケンス）を見つける問題を解決するための一般的なアルゴリズムを確認します。 実際の使用ではなく学習に重点が置かれているため、実装の言語としてPythonが選択されました。これにより、コード量が削減され、主要なアイデアに集中できます。

定義

シーケンスは、要素の順序付けられたセットです。 文字列はシーケンスの特殊なケースであり、単純化のためにさらに例を検討しますが、変更することなく、任意のテキストまたはその他の一貫性のあるものにも使用できます。

要素x ₁ x ₂ ... x _mで構成されるシーケンスxと、要素y ₁ y ₂ ... y _nで構成されるシーケンスyが _あるとします。 zが取得される元の要素xのインデックスの厳密に増加するセットが存在する場合、 zはxのサブシーケンスです。

xとyの 共通のサブシーケンスは 、 xの サブシーケンスとyのサブシーケンスの両方であるシーケンスzです。

最大合計サブシーケンスは、最大長の合計サブシーケンスです。 さらに本文では、略語LCSを使用します。

例として、 x = HA B R AHA BR 、 y = HARB OU R 、この場合LCS（x、y）= HARBRとします。 すでにLCS計算アルゴリズムに直接アクセスできますが、なぜこれが必要なのかを理解しておくといいでしょう。

実用化

最も一般的な用途は、GNU diffなどのファイル比較プログラムでの使用です。 2つのテキストのLCSを見つけたら、基本的な変更のリストをコンパイルしてxをyに、またはその逆に変換するのは簡単な作業です。 ボーナスとして、サブシーケンス全体の長さに基づいて、2つのシーケンスの類似性を判断するためのメトリックを指定できます。 それだけです、今では間違いなくビジネスに取りかかることができます。

最初のアプローチまたは民芸

いくつかの観察を開始するには：

1. シーケンスxおよびyについてLCS（x、y）= zをすでに計算している場合、同じ要素を追加してxおよびyから取得したシーケンスのLCSは、zとこの追加された要素で構成されます。
2. シーケンスxとyに1つの異なる要素を追加する場合、取得されたxaとybの場合、LCSはLCS（x、yb）またはLCS（xa、y）の2つのうち大きい方になります。

これらの観察はすでに再帰を実装するのに十分です：

def LCS_RECURSIVE(x, y): if len(x) == 0 or len(y) == 0: return [] if x[-1] == y[-1]: return LCS_RECURSIVE(x[:-1], y[:-1]) + [x[-1]] else: left = LCS_RECURSIVE(x[:-1], y) right = LCS_RECURSIVE(x, y[:-1]) return left if len(left) > len(right) else right
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

これは、この実装では当てはまらないと思うかもしれません。 最悪の場合、xとyの間に同じLCS_RECURSIVE要素がない場合、再帰レベルlen（x）+ len（y）で2 ^{len（x）+ len（y）}回が呼び出されます。 パフォーマンスに目を閉じても、コードは次のようになります。

 from uuid import uuid4 x = [uuid4().hex for _ in xrange(500)] y = [uuid4().hex for _ in xrange(500)] print LCS_RECURSIVE(x,y)
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

追加の呼び出しがなければ、sys.setrecursionlimitは失敗します。 最も成功した実装ではありません。

動的プログラミングまたは64 kbで十分です

問題のアルゴリズムは、Needleman-Wunschアルゴリズムとも呼ばれます。

全体のアプローチは、行がx要素で、列がy要素である、マトリックスの段階に要約されます。 記入時には、すでに行われた観察から生じる2つの規則が適用されます。

1.要素x _iがy _jに等しい場合、セル（i、j）にセル（i-1、j-1）の値が1を加えて書き込まれます

2.要素x _iが y _jと等しくない場合、値（i-1、j）および（i、j-1）の最大値がセル（i、j）に記録されます。

充填はiとjの二重サイクルで行われ、値が1ずつ増加します。したがって、各反復でこのステップで必要なセル値はすでに計算されています。

 def fill_dyn_matrix(x, y): L = [[0]*(len(y)+1) for _ in xrange(len(x)+1)] for x_i,x_elem in enumerate(x): for y_i,y_elem in enumerate(y): if x_elem == y_elem: L[x_i][y_i] = L[x_i-1][y_i-1] + 1 else: L[x_i][y_i] = max((L[x_i][y_i-1],L[x_i-1][y_i])) return L
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

何が起こっているのかの実例：



直接増加したセルは強調表示されます。 これらのセルを接続してマトリックスに入力すると、目的のLCSが得られます。 この場合、セルが強調表示されている場合は最大インデックスから最小インデックスに移動しながら接続する必要があります。そうでない場合は、対応する要素をLCSに追加し、マトリックス内の大きい値に応じて上または左に移動します：

 def LCS_DYN(x, y): L = fill_dyn_matrix(x, y) LCS = [] x_i,y_i = len(x)-1,len(y)-1 while x_i >= 0 and y_i >= 0: if x[x_i] == y[y_i]: LCS.append(x[x_i]) x_i, y_i = x_i-1, y_i-1 elif L[x_i-1][y_i] > L[x_i][y_i-1]: x_i -= 1 else: y_i -= 1 LCS.reverse() return LCS
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

アルゴリズムの複雑さはO（len（x）* len（y））、同じメモリスコアです。 したがって、100,000行の2行のファイルを1行ずつ比較する場合は、10 ^10個の要素のマトリックスをメモリに保存する必要があります。 つまり 実際に使用すると、ほとんど突然のMemoryErrorを受け取る恐れがあります。 次のアルゴリズムに進む前に、要素xの各反復で行列Lを埋めるとき、前の移動で取得した行のみが必要であることに注意してください。 つまり LCS自体を計算することなくLCSの長さを見つけることだけにタスクが制限されている場合、メモリ使用量をO（len（y））に減らして、行列Lの2行のみを保存することができます。

場所またはHirschbergアルゴリズムのために戦う

このアルゴリズムの背後にある考え方は単純です。入力シーケンスx = x ₁ x ₂ ... x _mを、境界インデックスiの任意の2つの部分にxb = x ₁ x ₂ ... x _iおよびxe = x _{i + 1で}分割する場合x _{i + 2} ... x _m 、次に同様にyをパーティション分割する方法があります（yをyb = y ₁ y ₂ ... y _jとye = y _{j + 1} y _{j + 2} ...に分割するインデックスjがありますy _n ）LCS（x、y）= LCS（xb、yb）+ LCS（xe、ye）

このパーティションyを見つけるために提案されます：

1. xsとyのfill_dyn_matrixで動的行列Lを塗りつぶします。 L1計算された行列Lの最後の行と同等
2. * xeおよび* yの行列Lに入力します（Pythonの観点から、xeおよびyの逆シーケンス：リスト（反転（x））、リスト（反転（y）））。 * L2は計算された行列Lの最後の行と等しく、L2は* L2からの川です
3. jを、合計L1 [j] + L2 [j]が最大になるインデックスと等しくします

これは、2つの反対側の端から行列Lを満たすことで表すことができます。



マトリックスLの最後の行にのみ必要があるため、計算中にマトリックス全体を保存する必要がないことに注意してください。 fill_dyn_matrixの実装をわずかに変更すると、次のようになります。

 def lcs_length(x, y): curr = [0]*(1 + len(y)) for x_elem in x: prev = curr[:] for y_i, y_elem in enumerate(y): if x_elem == y_elem: curr[y_i + 1] = prev[y_i] + 1 else: curr[y_i + 1] = max(curr[y_i], prev[y_i + 1]) return curr
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

次に、アルゴリズム自体について直接説明します。 古典的な分割統治です：シーケンスxには要素がありますが、xを半分に分割し、yの適切なパーティションを見つけて、シーケンスのペア（x​​b、yb）と（xe、ye）の再帰呼び出しの合計を返します。 些細なケースでは、xが1つの要素で構成され、yにある場合、この単一の要素xから単純にシーケンスを返すことに注意してください。

 def LCS_HIRSHBERG(x, y): x_len = len(x) if x_len == 0: return [] elif x_len == 1: if x[0] in y: return [x[0]] else: return [] else: i = x_len // 2 xb, xe = x[:i], x[i:] L1 = lcs_length(xb, y) L2 = reversed(lcs_length(xe[::-1], y[::-1])) SUM = (l1 + l2 for l1,l2 in itertools.izip(L1, L2)) _, j = max((sum_val, sum_i) for sum_i, sum_val in enumerate(SUM)) yb, ye = y[:j], y[j:] return LCS_HIRSHBERG(xb, yb) + LCS_HIRSHBERG(xe, ye)
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

O（len（y））のメモリ要件があり、計算に必要な時間は2倍になりましたが、O（len（x）len（y））は漸近的です。

ハント-シマンスキーアルゴリズム

最初に必要なことは、マッチリストテーブルのハッシュを作成することです。これには、2番目のシーケンスの要素のインデックスが含まれます。 これに必要な時間はO（len（y））です。 次のPythonコードがこれを実行します。

 y_matchlist = {} for index, elem in enumerate(y): indexes = y_matchlist.setdefault(elem, []) indexes.append(index) y_matchlist[elem] = indexes
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

HARBORシーケンスの場合、ハッシュは次のようになります{'A'：[1]、 'B'：[3]、 'H'：[0]、 'O'：[4]、 'R'：[2、6] 、「U」：[5]}。



さらに、THRESH [k-1] <y_indexおよびy_index <THRESH [k]の場合、THRESHのk番目の要素の値がインデックスy_indexになるように、シーケンスxの要素をソートして、用意されたマッチリストの対応するインデックスでTHRESH配列を埋めます。 したがって、THRESH配列はいつでもソートされ、バイナリ検索を使用して適切なkを見つけることができます。 THRESH要素を更新するとき、y_indexインデックスに対応するシーケンス要素もLCSに追加します。 次のコードで明確にできます。

 x_length, y_length = len(x), len(y) min_length = min(x_length, y_length) THRESH = list(itertools.repeat(y_length, min_length+1)) LINK_s1 = list(itertools.repeat(None, min_length+1)) LINK_s2 = list(itertools.repeat(None, min_length+1)) THRESH[0], t = -1, 0 for x_index,x_elem in enumerate(x): y_indexes = y_matchlist.get(x_elem, []) for y_index in reversed(y_indexes): k_start = bisect.bisect_left(THRESH, y_index, 1, t) k_end = bisect.bisect_right(THRESH, y_index, k_start, t) for k in xrange(k_start, k_end+2): if THRESH[k-1] < y_index and y_index < THRESH[k]: THRESH[k] = y_index LINK_x[k] = (x_index, LINK_x[k-1]) if k > t: t = k
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

その後、LINK_xからのみシーケンス自体を収集できます。

このアルゴリズムの実行時間はO（（r + n）log n）です。ここで、nは大きいシーケンスの長さで、rは一致の数です。最悪の場合、r = n * nの場合、前のソリューションよりも実行時間が長くなります。 メモリ要件は、O（r + n）のオーダーのままです。

合計

この問題を解決するアルゴリズムは多数ありますが、漸近的には、最も効率的なアルゴリズム（MasekとPaterson）のO（n * n / log n）時間推定があります。 LCS計算の全体的な効率が低いことを考えると、実際には、アルゴリズムの実行前に最も単純な準備が実行されます。たとえば、シーケンスの最初と最後で同一の要素を破棄し、シーケンス間の些細な違いを検索します。 また、動作時間の漸近的な動作に影響を与えないビット演算を使用した最適化がいくつかあります。

サンプルを含むすべてのコードをダウンロードする