程式員編程藝術第十一章：最長公用子序列(LCS)問題

    程式員編程藝術第十一章：最長公用子序列(LCS)問題

0、前言

    程式員編程藝術系列重新開始創作了（前十章，請參考程式員編程藝術第一~十章集錦與總結）。回顧之前的前十章，有些代碼是值得商榷的，因當時的代碼只顧闡述演算法的原理或思想，所以，很多的與代碼規範相關的問題都未能做到完美。日後，會著力修善之。

    搜遍網上，講解這個LCS問題的文章不計其數，但大多給讀者一種並不友好的感覺，稍感晦澀，且代碼也不夠清晰。本文力圖避免此些情況。力保通俗，闡述詳盡。同時，經典演算法研究系列的第三章（三、dynamic programming）寫的極其糟糕，所以，也算是對那文的一種彌補。有任何問題，歡迎不吝賜教。

第一節、問題描述

    什麼是最長公用子序列呢?好比一個數列 S，如果分別是兩個或多個已知數列的子序列，且是所有符合此條件序列中最長的，則S 稱為已知序列的最長公用子序列。

    舉個例子，如：有兩條隨機序列，如 1 3 4 5 5 ，and 2 4 5 5 7 6，則它們的最長公用子序列便是：4 5 5。

第二節、LCS問題的解決思路

• 窮舉法   

    解最長公用子序列問題時最容易想到的演算法是窮舉搜尋法，即對X的每一個子序列，檢查它是否也是Y的子序列，從而確定它是否為X和Y的公用子序列，並且在檢查過程中選出最長的公用子序列。X和Y的所有子序列都檢查過後即可求出X和Y的最長公用子序列。X的一個子序列相應於下標序列{1, 2, …, m}的一個子序列，因此，X共有2m個不同子序列（Y亦如此，如為2^n），從而窮舉搜尋法需要指數時間（2^m * 2^n）。

• 動態規划算法

    事實上，最長公用子序列問題也有最優子結構性質。

記:

Xi=﹤x1，⋯，xi﹥即X序列的前i個字元 (1≤i≤m)（首碼）

Yj=﹤y1，⋯，yj﹥即Y序列的前j個字元 (1≤j≤n)（首碼）

假定Z=﹤z1，⋯，zk﹥∈LCS(X , Y)。

• 若xm=yn（最後一個字元相同），則不難用反證法證明：該字元必是X與Y的任一最長公用子序列Z（設長度為k）的最後一個字元，即有zk = xm = yn 且顯然有Zk-1∈LCS(Xm-1 , Yn-1)即Z的首碼Zk-1是Xm-1與Yn-1的最長公用子序列。此時，問題化歸成求Xm-1與Yn-1的LCS（LCS(X , Y)的長度等於LCS(Xm-1 , Yn-1)的長度加1）。

• 若xm≠yn，則亦不難用反證法證明：要麼Z∈LCS(Xm-1, Y)，要麼Z∈LCS(X , Yn-1)。由於zk≠xm與zk≠yn其中至少有一個必成立，若zk≠xm則有Z∈LCS(Xm-1 , Y)，類似的，若zk≠yn 則有Z∈LCS(X , Yn-1)。此時，問題化歸成求Xm-1與Y的LCS及X與Yn-1的LCS。LCS(X , Y)的長度為：max{LCS(Xm-1 , Y)的長度, LCS(X , Yn-1)的長度}。

    由於上述當xm≠yn的情況中，求LCS(Xm-1 , Y)的長度與LCS(X , Yn-1)的長度，這兩個問題不是相互獨立的：兩者都需要求LCS(Xm-1，Yn-1)的長度。另外兩個序列的LCS中包含了兩個序列的首碼的LCS，故問題具有最優子結構性質考慮用動態規劃法。

    也就是說，解決這個LCS問題，你要求三個方面的東西：1、LCS（Xm-1，Yn-1）+1；2、LCS（Xm-1，Y），LCS（X，Yn-1）；3、max{LCS（Xm-1，Y），LCS（X，Yn-1）}。

    行文至此，其實對這個LCS的動態規劃解法已敘述殆盡，不過，為了成書的某種必要性，下面，我試著再多加詳細闡述這個問題。

第三節、動態規划算法解LCS問題

3.1、最長公用子序列的結構

    最長公用子序列的結構有如下表示：

    設序列X=<x1, x2, …, xm>和Y=<y1, y2, …, yn>的一個最長公用子序列Z=<z1, z2, …, zk>，則：

1. 若xm=yn，則zk=xm=yn且Zk-1是Xm-1和Yn-1的最長公用子序列；
2. 若xm≠yn且zk≠xm ，則Z是Xm-1和Y的最長公用子序列；
3. 若xm≠yn且zk≠yn ，則Z是X和Yn-1的最長公用子序列。

    其中Xm-1=<x1, x2, …, xm-1>，Yn-1=<y1, y2, …, yn-1>，Zk-1=<z1, z2, …, zk-1>。

3、2.子問題的遞迴結構

    由最長公用子序列問題的最優子結構性質可知，要找出X=<x1, x2, …, xm>和Y=<y1, y2, …, yn>的最長公用子序列，可按以下方式遞迴地進行：當xm=yn時，找出Xm-1和Yn-1的最長公用子序列，然後在其尾部加上xm(=yn)即可得X和Y的一個最長公用子序列。當xm≠yn時，必須解兩個子問題，即找出Xm-1和Y的一個最長公用子序列及X和Yn-1的一個最長公用子序列。這兩個公用子序列中較長者即為X和Y的一個最長公用子序列。

    由此遞迴結構容易看到最長公用子序列問題具有子問題重疊性質。例如，在計算X和Y的最長公用子序列時，可能要計算出X和Yn-1及Xm-1和Y的最長公用子序列。而這兩個子問題都包含一個公用子問題，即計算Xm-1和Yn-1的最長公用子序列。

    與矩陣連乘積最優計算次序問題類似，我們來建立子問題的最優值的遞迴關係。用c[i,j]記錄序列Xi和Yj的最長公用子序列的長度。其中Xi=<x1, x2, …, xi>，Yj=<y1, y2, …, yj>。當i=0或j=0時，空序列是Xi和Yj的最長公用子序列，故c[i,j]=0。其他情況下，由定理可建立遞迴關係如下：

3、3.計算最優值

    直接利用上節節末的遞迴式，我們將很容易就能寫出一個計算c[i,j]的遞迴演算法，但其計算時間是隨輸入長度指數增長的。由於在所考慮的子問題空間中，總共只有θ(m*n)個不同的子問題，因此，用動態規划算法自底向上地計算最優值能提高演算法的效率。

    計算最長公用子序列長度的動態規划算法LCS_LENGTH(X,Y)以序列X=<x1, x2, …, xm>和Y=<y1, y2, …, yn>作為輸入。輸出兩個數組c[0..m ,0..n]和b[1..m ,1..n]。其中c[i,j]儲存Xi與Yj的最長公用子序列的長度，b[i,j]記錄指示c[i,j]的值是由哪一個子問題的解達到的，這在構造最長公用子序列時要用到。最後，X和Y的最長公用子序列的長度記錄於c[m,n]中。

Procedure LCS_LENGTH(X,Y);<br />begin<br /> m:=length[X];<br /> n:=length[Y];<br /> for i:=1 to m do c[i,0]:=0;<br /> for j:=1 to n do c[0,j]:=0;<br /> for i:=1 to m do<br /> for j:=1 to n do<br /> if x[i]=y[j] then<br /> begin<br /> c[i,j]:=c[i-1,j-1]+1;<br /> b[i,j]:="";<br /> end<br /> else if c[i-1,j]≥c[i,j-1] then<br /> begin<br /> c[i,j]:=c[i-1,j];<br /> b[i,j]:="↑";<br /> end<br /> else<br /> begin<br /> c[i,j]:=c[i,j-1];<br /> b[i,j]:="←"<br /> end;<br /> return(c,b);<br />end; 

    由於每個數組單元的計算耗費Ο(1)時間，演算法LCS_LENGTH耗時Ο(mn)。

3、4.構造最長公用子序列

    由演算法LCS_LENGTH計算得到的數組b可用於快速構造序列X=<x1, x2, …, xm>和Y=<y1, y2, …, yn>的最長公用子序列。首先從b[m,n]開始，沿著其中的箭頭所指的方向在數組b中搜尋。當b[i,j]中遇到""時（意味著xi=yi是LCS的一個元素），表示Xi與Yj的最長公用子序列是由Xi-1與Yj-1的最長公用子序列在尾部加上xi得到的子序列；當b[i,j]中遇到"↑"時，表示Xi與Yj的最長公用子序列和Xi-1與Yj的最長公用子序列相同；當b[i,j]中遇到"←"時，表示Xi與Yj的最長公用子序列和Xi與Yj-1的最長公用子序列相同。這種方法是按照反序來找LCS的每一個元素的。

    下面的演算法LCS(b,X,i,j)實現根據b的內容列印出Xi與Yj的最長公用子序列。通過演算法的調用LCS(b,X,length[X],length[Y])，便可列印出序列X和Y的最長公用子序列。

Procedure LCS(b,X,i,j);<br />begin<br /> if i=0 or j=0 then return;<br /> if b[i,j]="" then<br /> begin<br /> LCS(b,X,i-1,j-1);<br /> print(x[i]); {列印x[i]}<br /> end<br /> else if b[i,j]="↑" then LCS(b,X,i-1,j)<br /> else LCS(b,X,i,j-1);<br />end;

在演算法LCS中，每一次的遞迴調用使i或j減1，因此演算法的計算時間為O(m+n)。

例如，設所給的兩個序列為X=<A，B，C，B，D，A，B>和Y=<B，D，C，A，B，A>。由演算法LCS_LENGTH和LCS計算出的結果如所示：

    我來說明下此圖（參考演算法導論）。在序列X={A，B，C，B，D，A，B}和 Y={B，D，C，A，B，A}上，由LCS_LENGTH計算出的表c和b。第i行和第j列中的方塊包含了c[i，j]的值以及指向b[i，j]的箭頭。在c[7,6]的項4，表的右下角為X和Y的一個LCS<B，C，B，A>的長度。對於i，j>0，項c[i，j]僅依賴於是否有xi=yi，及項c[i-1，j]和c[i，j-1]的值，這幾個項都在c[i，j]之前計算。為了重構一個LCS的元素，從右下角開始跟蹤b[i，j]的箭頭即可，這條路徑標示為陰影，這條路徑上的每一個“”對應於一個使xi=yi為一個LCS的成員的項（高亮標示）。

    所以根據上述圖所示的結果，程式將最終輸出：“B C B A”。

3、5.演算法的改進

    對於一個具體問題，按照一般的演算法設計策略設計出的演算法，往往在演算法的時間和空間需求上還可以改進。這種改進，通常是利用具體問題的一些特殊性。

    例如，在演算法LCS_LENGTH和LCS中，可進一步將數組b省去。事實上，數組元素c[i,j]的值僅由c[i-1,j-1]，c[i-1,j]和c[i,j-1]三個值之一確定，而數組元素b[i,j]也只是用來指示c[i,j]究竟由哪個值確定。因此，在演算法LCS中，我們可以不藉助於數組b而藉助於數組c本身臨時判斷c[i,j]的值是由c[i-1,j-1]，c[i-1,j]和c[i,j-1]中哪一個數值元素所確定，代價是Ο(1)時間。既然b對於演算法LCS不是必要的，那麼演算法LCS_LENGTH便不必儲存它。這一來，可節省θ(mn)的空間，而LCS_LENGTH和LCS所需要的時間分別仍然是Ο(mn)和Ο(m+n)。不過，由於數組c仍需要Ο(mn)的空間，因此這裡所作的改進，只是在空間複雜性的常數因子上的改進。

    另外，如果只需要計算最長公用子序列的長度，則演算法的空間需求還可大大減少。事實上，在計算c[i,j]時，只用到數組c的第i行和第i-1行。因此，只要用2行的數組空間就可以計算出最長公用子序列的長度。更進一步的分析還可將空間需求減至min(m, n)。

第四節、編碼實現LCS問題

    動態規劃的一個計算最長公用子序列的方法如下，以兩個序列 X、Y 為例子：

設有二維數組 f[i][j] 表示
X 的 i 位和
Y 的 j 位之前的最長公用子序列的長度，則有：

f[1][1] = same(1,1)

f[i][j] = max{f[i − 1][j − 1] +same(i,j), f[i − 1][j] ,f[i][j
− 1]}

其中，same(a,b)當
X 的第 a 位與
Y 的第 b 位完全相同時為“1”，否則為“0”。

此時，f[i][j]中最大的數便是
X 和 Y 的最長公用子序列的長度，依據該數組回溯，便可找出最長公用子序列。

該演算法的空間、時間複雜度均為O(n2)，經過最佳化後，空間複雜度可為O(n)，時間複雜度為O(nlogn)。

以下是此演算法的java代碼：

<br />import java.util.Random;</p><p>public class LCS{<br /> public static void main(String[] args){</p><p> //設定字串長度<br /> int substringLength1 = 20;<br /> int substringLength2 = 20; //具體大小可自行設定</p><p> // 隨機產生字串<br /> String x = GetRandomStrings(substringLength1);<br /> String y = GetRandomStrings(substringLength2);</p><p> Long startTime = System.nanoTime();<br /> // 構造二維數組記錄子問題x[i]和y[i]的LCS的長度<br /> int[][] opt = new int[substringLength1 + 1][substringLength2 + 1];</p><p> // 動態規劃計算所有子問題<br /> for (int i = substringLength1 - 1; i >= 0; i--){<br /> for (int j = substringLength2 - 1; j >= 0; j--){<br /> if (x.charAt(i) == y.charAt(j))<br /> opt[i][j] = opt[i + 1][j + 1] + 1; //參考上文我給的公式。<br /> else<br /> opt[i][j] = Math.max(opt[i + 1][j], opt[i][j + 1]); //參考上文我給的公式。<br /> }<br /> }</p><p> -------------------------------------------------------------------------------------</p><p> 理解上段，參考上文我給的公式：</p><p> 根據上述結論，可得到以下公式，</p><p> 如果我們記字串Xi和Yj的LCS的長度為c[i,j]，我們可以遞迴地求c[i,j]：</p><p> / 0 if i<0 or j<0<br /> c[i,j]= c[i-1,j-1]+1 if i,j>=0 and xi=xj<br /> / max(c[i,j-1],c[i-1,j] if i,j>=0 and xi≠xj</p><p> -------------------------------------------------------------------------------------</p><p> System.out.println("substring1:"+x);<br /> System.out.println("substring2:"+y);<br /> System.out.print("LCS:");</p><p> int i = 0, j = 0;<br /> while (i < substringLength1 && j < substringLength2){<br /> if (x.charAt(i) == y.charAt(j)){<br /> System.out.print(x.charAt(i));<br /> i++;<br /> j++;<br /> } else if (opt[i + 1][j] >= opt[i][j + 1])<br /> i++;<br /> else<br /> j++;<br /> }<br /> Long endTime = System.nanoTime();<br /> System.out.println(" Totle time is " + (endTime - startTime) + " ns");<br /> }</p><p> //取得定長隨機字串<br /> public static String GetRandomStrings(int length){<br /> StringBuffer buffer = new StringBuffer("abcdefghijklmnopqrstuvwxyz");<br /> StringBuffer sb = new StringBuffer();<br /> Random r = new Random();<br /> int range = buffer.length();<br /> for (int i = 0; i < length; i++){<br /> sb.append(buffer.charAt(r.nextInt(range)));<br /> }<br /> return sb.toString();<br /> }<br />}
第五節、面試考察的熱點

    在很多的面試過程中，考察的都是細節問題，如本文所闡述的LCS問題就經常出現在各大公司的面試中。請看下題：

56.最長公用字串
題目：如果字串一的所有字元按其在字串中的順序出現在另外一個字串二中，則字串一稱之為字串二的子串。注意，並不要求子串（字串一）的字元必須連續出現在字串二中。請編寫一個函數，輸入兩個字串，求它們的最長公用子串，並列印出最長公用子串。
例如：輸入兩個字串 BDCABA 和 ABCBDAB，字串 BCBA 和 BDAB 都是是它們的最長公用子串，則輸出它們的長度 4，並列印任意一個子串。
分析：求最長公用子串（Longest Common Subsequence, LCS）是一道非常經典的動態規劃題，因此一些重視演算法的公司像 MicroStrategy 都把它當作面試題。

    ok，此問題的分析上文已經敘述殆盡，代碼也已經給出，不再贅述。僅僅補上一個LCS問題的c/c++版本，如下（有任何問題或錯誤，懇請指正）：

#include<iostream><br />#include<stdio.h><br />#include<stdlib.h><br />#define LEN 501<br />using namespace std;<br />int b[LEN][LEN];<br />void lcs(int i,int j,char x[],int b[][LEN]);<br />int lcslength(char x[],char y[],int b[][LEN]){<br />int lenx=strlen(x),leny=strlen(y);<br />int c[LEN][LEN]={0};<br />for(int i=1;i<=lenx;i++){<br />for(int j=1;j<=leny;j++){<br />if(x[i-1]==y[j-1]){<br />c[i][j]=c[i-1][j-1]+1;<br />b[i][j]=1;<br />}<br />else{<br />if(c[i-1][j]>=c[i][j-1]){<br />c[i][j]=c[i-1][j];<br />b[i][j]=2;<br />}<br />else{<br />c[i][j]=c[i][j-1];<br />b[i][j]=3;<br />}<br />}<br />}<br />}<br />return c[lenx][leny];<br />} </p><p>int main(){<br /> char x[LEN],y[LEN];<br /> while(scanf("%s%s",x,y)!=EOF){<br /> printf("%d\n",lcslength(x,y,b));<br /> lcs(strlen(x),strlen(y),x,b);<br /> }<br /> return 0;<br />}<br />void lcs(int i,int j,char x[],int b[][LEN]){<br /> if(i==0 || j==0)return;<br /> if(b[i][j]==1){<br /> lcs(i-1,j-1,x,b);<br /> printf("%c",x[i-1]);<br /> }<br /> else if(b[i][j]==2)lcs(i-1,j,x,b);<br /> else lcs(i,j-1,x,b);<br />}
第六節、改進的演算法

    下面咱們來瞭解一種不同於動態規劃法的一種新的求解最長公用子序列問題的方法,該演算法主要是把求解公用字串問題轉化為求解矩陣L(p,m)的問題，在利用定理求解矩陣的元素過程中（1）while(i<k),L(k,i)=null，
                  （2）while(L(k,i)=k),L(k,i+1)=L(k,i+2)=…L(k,m)=k；

    求出每列元素，一直到發現第p+1 行都為null 時退出迴圈，得出矩陣L(k,m)後，B[L(1,m-p+1)]B[L(2,m-p+2)]…B[L(p,m)]即為A 和B 的LCS，其中p 為LCS 的長度。

6.1 主要定義及定理

• 定義 1 子序列(Subsequence)：給定字串A=A[1]A[2]…A[m]，(A[i]是A 的第i 個字母，A[i]∈字元集Σ，l<= i<m = A ， A 表示字串A 的長度)，字串B 是A 的子序列是指B=A[ 1 i ]A[ 2 i ]…A[ k i ],其中1 i < 2 i <…< k i 且k<=m.
• 定義2 公用子序列(Common Subsequence)：給定字串A、B、C，C 稱為A 和B 的公用子序列是指C 既是A 的子序列，又是B 的子序列。
• 定義3 最長公用子序列(Longest Common Subsequence 簡稱LCS)：給定字串A、B、C，C 稱為A 和B 的最長公用子序列是指C 是A 和B 的公用子序列，且對於A 和B 的任意公用子序列D，都有D <= C 。給定字串A 和B，A =m，B =n，不妨設m<=n，LCS 問題就是要求出A 和B 的LCS。
• 定義4 給定字串A=A[1]A[2]…A[m]和字串B=B[1]B[2]…[n]，A( 1:i)表示A 的連續子序列A[1]A[2]…A[i]，同樣B(1:j)表示B 的連續子序列B[1]B[2]…[j]。Li(k)表示所有與字串A(1:i) 有長度為k 的LCS 的字串B(l:j) 中j 的最小值。用公式表示就是Li(k)=Minj(LCS(A(1:i)，B(l:j))=k) [3]。

定理1 ∀ i∈[1，m]，有Li(l)<Li(2)<Li(3)<…<Li(m) .
定理2 ∀i∈[l，m-1]，∀k∈[l，m]，有i 1 L + (k)<= i L (k).
定理3 ∀ i∈[l，m-1]， ∀ k∈[l，m-l]，有i L (k)< i 1 L + (k+l).
以上三個定理都不考慮Li(k)無定義的情況。
定理4[3] i 1 L + (k)如果存在，那麼它的取值必為: i 1 L + (k)=Min(j, i L (k))。這裡j 是滿足以下條件的最小整數:A[i+l]=B[j]且j> i L (k-1)。

    矩陣中元素L(k，i)=Li(k)，這裡(1<i<=m，1<k<=m)，null 表示L(k,i)不存在。當i<k 時，顯然L(k，i)不存在。
    設p=Maxk(L(k ， m) ≠ null) ， 可以證明L 矩陣中L(p,m) 所在的對角線,L(1,m-p+1),L(2,m-p+2)…L(p-1,m-1),L(p,m) 所對應的子序列B[L(1,m-p+1)]B[L(2,m-p+2)]…B[L(p,m)]即為A 和B 的LCS，p 為該LCS 的長度。這樣，LCS 問題的求解就轉化為對m m L × 矩陣的求解。

6.2 演算法思想
    根據定理,第一步求出第一行元素,L(1,1),L(1,2),…L(1,m),第二步求第二行,一直到發現第p+1 行都為null 為止。在計算過程中遇到i<k 時,L(k,i)=null, 及L(k,i)=k時,L(k,i+1)=L(k,i+2)=…L(k,m)=k。這樣,計算每行的時間複雜度為O(n),則整個時間複雜度為O(pn)。在求L 矩陣的過程中不用儲存整個矩陣,只需儲存當前行和上一行即可。空間複雜度為O(m+n)。

    下面給出一個例子來說明:給定字串A 和B，A=acdabbc，B=cddbacaba，(m= A =7，n= B =9)。按照定理給出的遞推公式，求出A 和B 的L 矩陣2，其中的$表示NULL。

    則A 和B 的LCS 為B[1]B[2]B[4]B[6]=cdbc,LCS 的長度為4。

6.3 演算法虛擬碼
演算法 L(A,B,L)
輸入 長度分別為m,n 的字串A,B
輸出 A,B 的最長公用子序列LCS

L(A,B,L){//字串A，B，所求矩陣L<br /> for(k=1;k<=m;k++){ //m 為A 的長度<br /> for(i=1;i<=m;i++){<br /> if(i<k) L[k][i]=N;//i<k 時,L(k,i)=null，N 代表無窮大<br /> if(L[k][i]==k)//L(k,i)=k 時,L(k,i+1)=L(k,i+2)=…L(k,m)=k<br /> for(l=i+1;l<=m;l++)<br /> { L[k][l]=k;<br /> Break;}<br /> for(j=1;j<=n;j++){//定理4 的實現<br /> if(A[i+1]==B[j]&&j>L[k-1][i]){<br /> L[k][i+1]=(j<L[k][i]?j:L[k][i]);<br /> break;<br /> }<br /> if(L[k][i+1]==0)<br /> L[k][i]=N;<br /> }<br /> if(L[k][m]==N)<br /> {p=k-1;break;}<br /> }<br /> p=k-1;<br />}<br />

6.4 結語
    本節主要描述區別於動態規劃法的一種新的求解最長公用子序列問題的方法，在不影響精確度的前提下，提高序列匹配的速度，根據定理i 1 L + (k)=Min(j, i L (k))得出矩陣，在求解矩陣的過程中對最耗時的L(p,m)進行條件約束最佳化。我們在Intel(R) Core(TM)2 Quad 雙核處理器、1G 記憶體，軟體環境：windows xp 下實驗結果證明，本文演算法與其他經典的比對演算法相比,不但能夠取得準確的結果,而且速度有了較大的提高（本節參考了劉佳梅女士的論文）。

    若有任何問題，懇請不吝指正。謝謝各位。完。