文本比较算法Ⅶ——线性空间求最长公共子序列的Nakatsu算法
在参阅《A Longest Common Subsequence Algorithm Suitable for Similar Text Strings》（Narao Nakatsu，Yahiko Kambayashi，Shuzo Yajima著）后。发现该算法可以利用线性空间求出最长公共子序列。该算法的时间占用O（n（m-p+1）），p为最长公共子序列的长度。
字符串A和字符串B，计算LCS(A,B)
定义一：设M=Len(A)，N=Len(B)，不妨设M≤N。
定义二：A=a1a2……aM，表示A是由a1a2……aM这M个字符组成
B=b1b2……bN，表示B是由b1b2……bN这N个字符组成
LCS(i,j)=LCS(a1a2……ai,b1b2……bj)，其中1≤i≤M，1≤j≤N
定义三：L(k,i)表示，所有与字符串a1a2……ai有长度为k的LCS的字符串b1b2……bj中j的最小值。
用公式表示就是：L(k,i)=Min｛j｝ Where LCS(i,j)=k
用一个例子来说明：A="CD"，B="CEFDRT"。
很明显的是LCS(2,1)=1，LCS(2,2)=1，LCS(2,3)=1。
满足LCS(2,j)=1这个条件的j有三个，分别是j=1、j=2、j=3。其中j最小值是1。故L(1,2)=1
为了推导L的计算，有下面几个定理。
定理一：任意的i，1≤i≤M。有L(1,i)＜L(2,i)＜L(3,i)……
定理二：任意的i，1≤i≤M-1。任意的k，1≤k≤M。有L(k,i+1)≤L(k,i)
定理三：任意的i，1≤i≤M-1。任意的k，1≤k≤M－1。有L(k,i)＜L(k+1,i+1)
定理四：如果L(k,i+1)存在，则L(k,i+1)的计算公式为
L(k,i+1)=Min｛Min｛j｝，L(k,i)｝ Where ｛ai+1=bj And j>L(k-1,i)｝
上面四个定理证明从略。可以从上面四个定理推导出L的计算。
故，L的计算公式为
①L(1,1)=Min｛j｝ Where ｛a1=bj｝
②L(1,i)=Min｛Min｛j｝ Where ｛ai=bj｝，L(1,i-1)｝ 此时，1<i≤M
③L(k,i)=Min｛Min｛j｝ Where ｛ai=bj And j＞L(k-1,i-1)｝，L(k,i-1)｝ 此时，1<i≤M，1<k≤M
注：以上公式中，若找不到满足Where后面条件的j，则j=MaxValue
当i＜k时，则L(k,i)=MaxValue
MaxValue是一个常量，表示“不存在”
在实际的算法实现中，为了简化运算，再次提出几个定义。
定义： L(0,i)=0 1≤i≤M
L(k,0)=MaxValue 1<k≤M
MaxValue=N+1 （只要定义一个j不可能取到的值就可以了）
则，在①中的公式可以写成
L(1,1)=Min｛j｝ Where ｛a1=bj｝=Min｛j｝ Where ｛a1=bj And j>0 }
=Min｛Min｛j｝ Where ｛a1=bj And j>0 ｝，MaxValue｝
=Min｛Min｛j｝ Where ｛a1=bj And j>L(0,0) ｝，L(1,0)｝
在②中的公式可以写成
L(1,i)=Min｛Min｛j｝ Where ｛ai=bj｝，L(1,i-1)｝
=Min｛Min｛j｝ Where ｛ai=bj And j>0｝，L(1,i-1)｝
=Min｛Min｛j｝ Where ｛ai=bj And j>L(0,i)｝，L(1,i-1)｝ 此时，1<i≤M
于是，三个公式统一了
④L(k,i)=Min｛Min｛j｝ Where ｛ai=bj And j＞L(k-1,i-1)｝，L(k,i-1)｝ 此时，1≤i≤M，1≤k≤M
且当i＜k时，则L(k,i)=MaxValue
仔细观察④，公式还可以写成如下
⑤ L(k,i)=Min｛j｝ Where ｛ai=bj And L(k-1,i-1)<j<L(k,i-1)｝ 1≤i≤M，1≤k≤M，且j存在
或 L(k,i)=L(k,i-1) 1≤i≤M，1≤k≤M，当j不存在时
写成⑤的目的有两个：一个是简化计算，不计算不必要的值；一个是为了标记，为后面计算最长公共子序列做准备。
接下来，将会用例子来说明：
A：481234781；B：4411327431
第一步：初始化L矩阵
| 4 | 8 | 1 | 2 | 3 | 4 | 7 | 8 | 1 | ||
|---|---|---|---|---|---|---|---|---|---|---|
| i=0 | i=1 | i=2 | i=3 | i=4 | i=5 | i=6 | i=7 | i=8 | i=9 | |
| k=0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| k=1 | V | |||||||||
| k=2 | V | V | ||||||||
| k=3 | V | V | V | |||||||
| k=4 | V | V | V | V | ||||||
| k=5 | V | V | V | V | V | |||||
| k=6 | V | V | V | V | V | V | ||||
| k=7 | V | V | V | V | V | V | V | |||
| k=8 | V | V | V | V | V | V | V | V | ||
| k=9 | V | V | V | V | V | V | V | V | V |
第二步：如表格所示，计算第一条对角线
| 4 | 8 | 1 | 2 | 3 | 4 | 7 | 8 | 1 | ||
|---|---|---|---|---|---|---|---|---|---|---|
| i=0 | i=1 | i=2 | i=3 | i=4 | i=5 | i=6 | i=7 | i=8 | i=9 | |
| k=0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| k=1 | V | 1 | ||||||||
| k=2 | V | V | V | |||||||
| k=3 | V | V | V | V | ||||||
| k=4 | V | V | V | V | V | |||||
| k=5 | V | V | V | V | V | V | ||||
| k=6 | V | V | V | V | V | V | V | |||
| k=7 | V | V | V | V | V | V | V | V | ||
| k=8 | V | V | V | V | V | V | V | V | V | |
| k=9 | V | V | V | V | V | V | V | V | V | V |
运气很差，只有第一个单元格有值，其余的都是V（MaxValue），很显然这不是问题的解。这条对角线满足L(k-1,i-1)<j<L(k,i-1)的只有一个单元格。先把它标记出来。
第三步：如表格所示，计算相邻的第二条对角线
| 4 | 8 | 1 | 2 | 3 | 4 | 7 | 8 | 1 | ||
|---|---|---|---|---|---|---|---|---|---|---|
| i=0 | i=1 | i=2 | i=3 | i=4 | i=5 | i=6 | i=7 | i=8 | i=9 | |
| k=0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| k=1 | V | 1 | 1 | |||||||
| k=2 | V | V | V | 3 | ||||||
| k=3 | V | V | V | V | 6 | |||||
| k=4 | V | V | V | V | V | 9 | ||||
| k=5 | V | V | V | V | V | V | V | |||
| k=6 | V | V | V | V | V | V | V | V | ||
| k=7 | V | V | V | V | V | V | V | V | V | |
| k=8 | V | V | V | V | V | V | V | V | V | V |
| k=9 | V | V | V | V | V | V | V | V | V | V |
运气比较好，有四个值。说明LCS(A,B)至少是4。但由于对角线上有V（MaxValue），所以本条对角线还不是解。同时，满足L(k-1,i-1)<j<L(k,i-1)的条件的单元格有3个。把它们标记出来。
第四步：如表格所示，计算相邻的第三条对角线
| 4 | 8 | 1 | 2 | 3 | 4 | 7 | 8 | 1 | ||
|---|---|---|---|---|---|---|---|---|---|---|
| i=0 | i=1 | i=2 | i=3 | i=4 | i=5 | i=6 | i=7 | i=8 | i=9 | |
| k=0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| k=1 | V | 1 | 1 | 1 | ||||||
| k=2 | V | V | V | 3 | 3 | |||||
| k=3 | V | V | V | V | 6 | 5 | ||||
| k=4 | V | V | V | V | V | 9 | 8 | |||
| k=5 | V | V | V | V | V | V | V | V | ||
| k=6 | V | V | V | V | V | V | V | V | V | |
| k=7 | V | V | V | V | V | V | V | V | V | V |
| k=8 | V | V | V | V | V | V | V | V | V | V |
| k=9 | V | V | V | V | V | V | V | V | V | V |
同理，这条对角线也不是解。把满足L(k-1,i-1)<j<L(k,i-1)的单元格标记出来。一共是两个。
第五步：如表格所示，计算相邻的第四条对角线
| 4 | 8 | 1 | 2 | 3 | 4 | 7 | 8 | 1 | ||
|---|---|---|---|---|---|---|---|---|---|---|
| i=0 | i=1 | i=2 | i=3 | i=4 | i=5 | i=6 | i=7 | i=8 | i=9 | |
| k=0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| k=1 | V | 1 | 1 | 1 | 1 | |||||
| k=2 | V | V | V | 3 | 3 | 3 | ||||
| k=3 | V | V | V | V | 6 | 5 | 5 | |||
| k=4 | V | V | V | V | V | 9 | 8 | 7 | ||
| k=5 | V | V | V | V | V | V | V | V | V | |
| k=6 | V | V | V | V | V | V | V | V | V | V |
| k=7 | V | V | V | V | V | V | V | V | V | V |
| k=8 | V | V | V | V | V | V | V | V | V | V |
| k=9 | V | V | V | V | V | V | V | V | V | V |
很遗憾，这个还不是解。满足L(k-1,i-1)<j<L(k,i-1)的解就只有1个。标记出来。
第六步：如表格所示，计算相邻的第五条对角线
| 4 | 8 | 1 | 2 | 3 | 4 | 7 | 8 | 1 | ||
|---|---|---|---|---|---|---|---|---|---|---|
| i=0 | i=1 | i=2 | i=3 | i=4 | i=5 | i=6 | i=7 | i=8 | i=9 | |
| k=0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| k=1 | V | 1 | 1 | 1 | 1 | 1 | ||||
| k=2 | V | V | V | 3 | 3 | 3 | 2 | |||
| k=3 | V | V | V | V | 6 | 5 | 5 | 5 | ||
| k=4 | V | V | V | V | V | 9 | 8 | 7 | 7 | |
| k=5 | V | V | V | V | V | V | V | V | V | 10 |
| k=6 | V | V | V | V | V | V | V | V | V | V |
| k=7 | V | V | V | V | V | V | V | V | V | V |
| k=8 | V | V | V | V | V | V | V | V | V | V |
| k=9 | V | V | V | V | V | V | V | V | V | V |
这条对角线上的值都不是V（MaxValue）。因此，问题的解出来了。注意到这条对角线对应的k=5，因此LCS(A,B)=5。
在表格中，满足L(k-1,i-1)<j<L(k,i-1)的单元格一共有9个（包括本条对角线的2个）。最大公共子序列的下标就在这9个单元格中。获取的依据是：每行取一个，每行的单元格都在上一行单元格的右边
按照这个依据，最长公共子序列的下标分别是
1；3；5；7；10 C=b1b3b5b7b10=41371
此时最佳匹配为： A：48123478_1
B：4411327431
1；3；6；7；10 C=b1b3b6b7b10=41271
此时最佳匹配为： A：481__23478_1
B：441132__7431
1；3；6；9；10 C=b1b3b6b9b10=41231
此时最佳匹配为： A：481__2__34781
B：441132743___1
1；3；5；8；10 C=b1b3b5b8b10=41341
此时最佳匹配为： A：48123__4781
B：441132743_1
1；3；6；8；10 C=b1b3b6b8b10=41241
此时最佳匹配为： A：481__234781
B：441132743_1
可以看出，Nakatsu算法很好的找到了所有的最长的公共子序列。但是，要找到所有的最长公共子序列，存储的空间是不会少于O（MN）的。因此，退而求其次，只求一个最长公共子序列。这个在空间的占用上是可以优化到O（M）的。
Nakatsu算法的时间分析：Nakatsu算法并没有把整个L矩阵计算出来，只是不停的尝试计算对角线，当遇到满足条件的对角线，就退出计算。假设LCS(A,B)=P。则可以判断出一共计算了M-P+1条对角线。而在每一条对角线的计算中，计算的次数满足下列的条件，假设计算第t条对角线
S=(L(1+t-1,1)-L(t-1,0))+(L(2+t-1,2)-L(1+t-1,1))……+(L(M+1-t,M)-L(M-t,M-1))
=L(M+1-t,M)-L(t-1,0)
=L(M+1-t,M)≤V=N+1
故Nakatsu算法的时间占用为O（（N+1）（M-P+1））。当P越是接近M，Nakatsu算法的速度越快。
没有优化过的Nakatsu算法占用空间O（MN）。因此，接下来的步骤就是对算法进行优化。前面提到，如果要计算出所有的最长公共子序列，空间是不可能压缩到线性空间的。因此，退而求其次，只求出一个最长公共子序列。通过观察，发现每一行第一个不是MaxValue的值组成的下标就是一个最长公共子序列。在本篇文章中就是1、3、6、9、10，对应的最长公共子序列是b1b3b6b9b10=41231。
接下来就是优化过程。注意到计算的过程是沿着对角线进行的。且每一个值的计算只和本条对角线和上一条对角线相关。因此，优化的思路就是把对角线转化为一维数组，也就是线性空间
第一步：初始化数组LL()和P()；
LL(0)=0
LL(i)=V 1≤i≤M
P(i)=V 1≤i≤M
此时，LL(0)表示L(0,0)；LL(1)表示L(1,0)；LL(2)表示L(2,1)；……
第二步：依次计算第一条对角线上的元素
用临时变量T计算L(1,1)；T=F(L(0,0),L(1,0))=F(LL(0),LL(1))。注：F表示某种运算关系
将T的值赋给LL(1)。此时LL(1)表示LL(1,1)，LL(2)表示L(2,1)；
重复上面的计算，直到计算完本条对角线
如果是第k行的第一个不为V的值，将该值赋给P(k)
第三步：第一条对角线计算完之后，此时，LL(0)表示L(0,1)；LL(1)表示L(1,1)；LL(2)表示L(2,2)；……
如果，这条对角线不是解，重复第二步，计算下一条对角线。直到遇到解为止。
不过要注意的是。第i条对角线只有m-i+1个元素。所以只计算到LL(m-i+1)。
某条对角线，某个元素是V的话，则这条对角线之后的元素都是V，就不需要计算了。
贴上VB2008的代码。 用这个代码计算上面的示例时，得到的最长公共子序列是41231。在优化到线性空间的情况下，只能得出一个最长公共子序列。这也是Nakatsu算法的本质——计算最长公共子序列。至于要得到最佳匹配，还得继续研究。
代码如下：
Public Class clsNakatsu
Private _A() As Char, _B() As Char
Public Sub New(ByVal A As String, ByVal B As String)
_A = A.ToCharArray
_B = B.ToCharArray
End Sub
Public Function GetLCS() As String
Dim L(_A.Length) As Integer, P(_A.Length) As Integer
Dim i As Integer, j As Integer, T As Integer
L(0) = -1
For i = 1 To _A.Length
L(i) = Integer.MaxValue
P(i) = Integer.MaxValue
Next
For i = 0 To _A.Length - 1
For j = 1 To _A.Length - i
T = GetP(i + j - 1, L(j - 1), L(j))
If T <> Integer.MaxValue Then L(j) = T
If P(j) = Integer.MaxValue AndAlso L(j) <> Integer.MaxValue Then P(j) = L(j)
If L(j) = Integer.MaxValue Then Exit For
Next
If L(_A.Length - i) <> Integer.MaxValue Then
Dim tS As New System.Text.StringBuilder
For j = 1 To _A.Length - i
tS.Append(_B(P(j)))
Next
Return tS.ToString
End If
Next
Return ""
End Function
Private Function GetP(ByVal Index As Integer, ByVal StartP As Integer, ByVal EndP As Integer)
Dim i As Integer
If EndP = Integer.MaxValue Then EndP = _B.Length
For i = StartP + 1 To EndP - 1
If _A(Index) = _B(i) Then Return i
Next
Return Integer.MaxValue
End Function
End Class
