动态规划近似串匹配

1、计算机算法设计与分析导论刘璟1Chapter 7. 动态规划（Dynamic Programming） 7.1 动态规划的基本原理 7.2 最优二分搜索树(Optimal Binary Search Tree) 7.3 近似串匹配(Approximate String Matching)问题27.1 动态规划的基本原理7.1.1 Fibonacci数的计算 Fibonacci数又称为Fibonacci数列，定义为：F0=0, F1=1, Fn=Fn-1 + Fn-2 (n 2)计算Fibonacci数列可由递归函数fibo完成。 递归函数fibo 由此可知，函数fibo的计算量随n的增加而急剧

2、增加，n=6时需25次调用，n=10时需177次调用，n=15时需1974次调用。进一步的研究表明，调用次数 An = 2Fn+1 - 1，其中 ， 。可以估计， ，其计算量是n的指数函数。3从Fig.7.1中可以看出，大量的调用过程是重复的，此算法可以改进。 函数fibo的改进函数fib这个程序的时间代价为O(n)阶。 47.1.2 矩阵连乘的顺序问题 1. 一个实例四个矩阵A1、A2、A3、A4相乘，设其阶数分别为A1：301，A2：140，A3：4010，A4：1025。 因为矩阵相乘满足结合律，所以可有下面五种（实际为六种）不同的运算次序，而且不同的运算次序所需的元素间乘法的次数不同，

3、如下面所列： ( ( A1 A2 ) A3 ) A4 30140+304010+301025=20700( A1 A2 )( A3 A4 ) 30140+401025+304025=41200( A1( A2 A3 ) ) A4 14010+30110+301025=8200A1( ( A2 A3 ) A4 ) 14010+11025+30125=1400A1( A2( A3 A4 ) ) 401025+14025+30125=117505五种运算次序的计算结果相同，但所花费的时间代价相差极大。如果某一应用问题需要经常进行矩阵连乘运算，应首先确定最优的矩阵连乘的次序。2. 最优矩阵连乘次序（O

4、ptimal Matrix Multiplication Order）问题给定n个矩阵A1,A2,.,An的维数为d0,d1,d2,.,dn，即：Ai的维数为di-1di(1in)。求一种连乘次序，使得计算时所需的元素乘法的总数最少。 3. 算法的思路如同上面实例中所做的那样，把所有可能的运算次序所需的计算量全计算出来，选择最小者，这种方法称为穷举法。不过，当n值较大时，这种方法的计算量过高。n个矩阵相乘应有（n-1）！种不同的运算次序，计算每种运算次序所需的时间代价需要2(n-1)次乘法和n-2次加法运算。当n=11时需要7257.6万次乘法。6由于不能事先决定最初在何处进行划分，所以分治策

5、略难以解决这个问题。但这个问题满足最优子结构性质：即，当选择了进行最外层相乘的位置之后，其左右两边的矩阵相乘序列都必须是时间代价最小的，可以考虑采用贪心策略。 一种使用贪心策略的解决方法是：每次优先选择其相乘代价最小的两个矩阵。例如在本实例中，A1A2A3A4有三个可能的相邻矩阵乘积，其中A2A3的代价400为最小。那么首先完成A2A3，在余下的两个可能的相邻矩阵乘积中：30110和11025相比较，后者较小。于是得到的解为A1( ( A2 A3 ) A4 )，与穷举法的最优结果一致。不过该贪心策略不能保证得到最优解。例如反例： 矩阵A、B、C，维数分别为：401，120，2050， (AB)

6、C 40120+402050=40800 A(BC) 12050+40150=30007另外一种采用贪心策略的方法是，对于n个矩阵A1. An的维数序列d0.dn，每次从d1.dn-1中取最大值di，首先进行Ai与Ai+1的乘法使最大的维数仅参加一次乘法运算，这样做有可能有利于减少矩阵连乘的计算量。使用这一策略，对上面的两个简单实例进行计算，其结果都是正确的。但是这个贪心算法仍然不能总是找出最优解。 可以从Fibonacci数的计算得到启发，采用动态规划的方法设计算法。最简单的递归算法描述如下： 递归算法MinCost 8此算法的递归过程中，存在与fibo函数相似的地方，即会有大量的重复调用，

7、其调用关系如Fig7.2所示，其中n=4。 9递归函数MinCost的调用过程在n=4时共27次，n=5时为81次，n=6时为249次，当n增大时，计算量急剧增加。如果采用自底向上的计算方式，函数MinCost(1,n)的计算代价将大大减少。其计算过程为：MinCost(1,1)=MinCost(2,2)=MinCost(3,3)=MinCost(4,4)=0MinCost(1,2)=MinCost(1,1)+MinCost(2,2)+d0d1d2=30140=1200类似地，Mincost(2,3)=d1d2d3=14010=400；Mincost(3,4)=d2d3d4=401025=10

8、000；在第二步计算MinCost(1,3)和MinCost(2,4)，最后计算MinCost(1,4)。由此可以看出，n=4时函数MinCost(1,4)的计算量是4+3+2+1=10次，n=5时为15次，n=6时为21次。 104. 矩阵连乘的DP算法设矩阵个数为n，维数为dimn+1（n+1个值），同时用数组MultiOrdern保存程序运行时得到的最优乘法顺序，可得： 算法7.1 最优矩阵连乘算法 MatrixOrder 函数MatrixOrder返回最优连乘次序所需的最小时间代价，同时将最优次序保存在数组MultiOrder中，从此数组中得到最优乘法次序的算法为：ExtractOrd

9、er 用上述算法计算本节给出的实例，其结果为：最小代价为cost0n=1400，最优乘法次序MultiOrder1.3=2,3,1，即A1(A2A3)A4)。 117.1.3 动态规划（DP）算法的基本条件 1. 最优子结构性质最优化原理。其特征是：当要求一个问题的最优解时，构成整体解的子问题的解也必须是最优的。例如，为了使计算n个矩阵连乘A1A2.An的代价最小，无论最后一次乘法的位置k在何处（1kn），其左右两部分A1.Ak，Ak+1.An的乘积也必须是代价最小的。这也可用最短路径问题来说明：若V1V2.Vn是一条从V1到Vn的最短路径，那么这条路径的任何一段，比如从Vi到Vj（1ijn）

10、也必须是一条最短路径。2. 子结构重迭性质 简单的递归程序解法都是一种自顶向下进行递归分解的过程，其中包含了大量的重复调用，在这种情况下采用动态规划方法特别有效。因此，问题中这种子结构重迭性质是采用动态规划方法的另一个条件。动态规划算法的一个特征是自底向上，它可以大幅度减少计算代价。12解决这类具有子结构重迭性质的问题还有另外一种被称为备忘录（memorization）方法的自顶向下的递归方式。其思想是：由程序设置“备忘录”，每计算出一个新的子结构的解时，都保存起来。当遇到一次递归时，判断是否已经计算，如果已经计算，只需取出先前保存的结果即可。这种方式与动态规划算法相比，要增加保存与检索中间结

11、果的代价。算法7.2 求Fibonacci数的备忘录算法 MemoFib 137.2 最优二分搜索树（Optimal Binary Search Tree）7.2.1 OBST问题 二分字典树是构造数据存储与检索系统的一种方便形式，例如一个翻译系统需要一个字典数据库。可以按照单词的字典顺序构造一个二分搜索树，能获得较高的搜索效率。 假定字典中只包含15个单词：and，cabbage，come，has，king，of，pig，ring，said，talk，the，thing，time，walrus，wing。按照字典顺序插入显然形成一个退化的二分树即线性链表，其平均搜索代价是(1+15)/2=8

12、次单词比较。14采用Fig7.3中的二分搜索树，最多仅需4次比较，平均需要(1+22+34+48)/15=3.17次比较，因此完全二分搜索树具有最小的平均搜索代价。 15在实用系统中，还要考虑到不同的单词在实际使用过程中被搜索的概率不同，例如，and、the等经常被搜索，而pig、cabbage等则很少被搜索。因此若把and、the放在树的最底层必然增大词典的平均搜索代价，这时的平均搜索代价应为：即，二分搜索树T的平均搜索代价应为从根到各个单词节点在T中的路长Ci与其查找概率Pi乘积之和。当各个单词的查找概率不同时，完全二分搜索树就不一定最优了。例如，假定词典仅由4个单词cat、come、of

13、、the组成，它们的查找概率分别为：cat:0.12，come:0.08，of:0.35，the:0.4516由这四个单词可以组成许多种不同的二分字典树，Fig7.4中给出了其中的三个。 按照上面给出的计算方法，三棵树（a,b,c）的平均搜索代价分别为：(a) 0.08+2(0.12+0.35)+30.45=2.37次比较；(a) 0.35+2(0.08+0.45)+30.12=1.77次比较；(a) 0.35+2(0.12+0.45)+30.08=1.73次比较。17平均搜索代价的差别说明，不同的二分搜索树在一定的查找概率条件下，性能是不同的。因此，在单词集合及其查找概率给定的条件下，构造一

14、个平均搜索代价最小的二分搜索树的意义是很明显的。 在实际问题中，还要考虑不成功的搜索，即查找给定单词集合之外的单词。这时可以把二分搜索树加以扩充，例如Fig7.4中的(c)，可以为这棵4个节点的树增加5个外部节点，形成一棵新树，如Fig7.5所示。其中4个内部节点表示单词集的4个单词，5个外部节点则表示检索其它单词时不成功的搜索路径。例如：搜索单词as，将到达最左边的外部节点，其代价为两次比较。 18最优二分搜索树（OBST）问题可以描述为： 已知：n个单词的查找概率p1,p2,.,pn和对应于n+1个外部节点的不成功查找概率q0,q1,q2,.,qn， 。 求：构造一种二分搜索树T，使平均搜

15、索代价A(T)最小。 解这个问题最简单的方法是把由n个单词（节点）组成的所有二分搜索树的平均搜索代价全算出来，取其最小者。不过，4个单词的二分搜索树有14种，7个单词的二分搜索树有429种就可知道，n较大时，这种方法是根本行不通的。 7.2.2 动态规划算法的思路设n个单词为a1,a2,.,an，则由它们构成二分搜索树T1,n可以表示为Fig7.6的形式。其中ak(1kn)为T1,n的根，其左子树T1,k-1由a1.ak组成，右子树Tk+1,n由ak+1.an组成。 19用代价函数Cost(low,high)表示由alow.ahigh及相应的外部节点blow-1,blow,.,bhigh组成的

16、二分搜索树的最小代价。特别地，当high=low-1时表示空树，Cost(i,i-1)=0。 Cost(low,high,r)表示指定以ar为根时，组成的二分搜索树的最小搜索代价。如果T1,n是最优的，那么T1,k-1，Tk+1,n两个子树也必然是最优的，因此，该问题具有最优子结构性质；另一方面，单词序列中任意几个相邻单词组成的子树将同时出现在多个包含它们的更大的树中，因此，子结构重迭性质也是很明显的。T1,k-1是由a1至ak-1元素构成子树中最优者20权函数W(low,high)表示相应二分搜索树的所有节点的搜索概率之和，称为该树的权：特别地：W(i,i)=pi+qi-1+qi；W(i,i

17、-1)=qi-1。（1in）。 以上各式有如下的递推关系：Cost(low,high,r) =pr+W(low,r-1)+Cost(low,r-1)+W(r+1,high)+Cost(r+1,high) =W(low,high)+Cost(low,r-1)+Cost(r+1，high)Cost(low,high) =MinCost(low,high,r)|lowrhigh依照自底向上的顺序，很容易计算出Cost(1,n)以及最优树的各级子树的根。(注意路径的长度如何计算的,如,根结点r,在计算中最终结果为pr*1,而其他结点,如在根的下层,则计算两次为pj+pj=pj*2) 第i个单词为树的权

18、第i个单词前一个未找到的为树的权217.2.3 算法OBST已知：n个单词a1.an的成功与不成功查找概率分别为：p1.pn，q0.qn。即：P1.n，Q0.n。求：由a1.an组成的最优二分搜索树。为此设置二维数组Cnn，Wnn，Rnn。算法7.3 最优二分搜索树构造算法 OBST 算法OBST的结果是数组C和数组R，其中C1n为求得的最优二分搜索树的搜索代价，R中包含了构造该树的信息。假定单词序列A1.n类型为word，树节点的类型为：struct node word key; node * left,* right;可得到：由数组R构造搜索树的算法 BuildTree227.2.4 算法

19、OBST的复杂度分析算法7.3中，T(n)=O(n3)，再有T(n)=(n3)。算法BuildTree表面上看是一个双递归程序，但由于算法每调用一次即生成一个（内部）节点，因此至多调用n次程序就结束，其复杂度是O(n)阶的。算法的空间代价主要体现在计算过程中所必需的二维数组，故有S(n)=(n2)。 7.2.5 讨论1. DP算法应用于OBST构造问题的效果非常明显。n个单词的二分字典树的个数是n的一个增长很快的指数函数。例如，10个节点的二分字典树就有16796种，因此穷举法是任何高速计算机都不能接受的。但DP算法把计算量降到了(n3)阶，而且可以估计出n3项的系数大约为1/6，当n较大时其

20、加速的倍数已经是天文数字了。由此可以看出DP算法的强大威力。 232. 构造二分搜索树的问题与矩阵连乘的顺序问题十分相似，关键是代价函数递推关系的确定，这也是许多应用动态规划方法设计有效算法的共同之处。同时，由于问题的不同，处理细节时又是有差别的。 3. 该算法仍然可以改进，D.Knuth提出了一个重要的改进方法，把最核心的关系式：Cost(low,high)=MinCost(low,high,r)|lowrhigh改进为：Cost(low,high)=MinCost(low,high,r)|R(low,high-1)rR(low+1,high)。也就是在算法7.3中的内层循环 for(r=l

21、ow;r=high;r+).改变为 for(r=Rlowhigh-1;r=Rlow+1high;r+).。24这一改动虽然使程序变化不大，但却明显地改变了计算量。这可以通过一个简化实例来说明：设有100个单词a1.a100，如果每个单词的查找概率相同，则应有R199=50，R2100=51，即原算法变量r从r=1到r=100循环100次，改进之后变量r从r=50到r=51循环2次，计算量减少几十倍。 改进后算法的时间复杂度可由下式估计，内层循环的计算次数为：因此算法7.3的时间复杂度为O(n2)。Knuth的改进，使计算复杂度从(n3)降为(n2)。25这一成果更为重要的部分是证明其结果的正确

22、性，即证明：Rij-1RijRi+1j,直观上看是显然的，但严格证明却相当困难。4. 关于二分字典树的研究还有许多成果，例如可以用贪心算法来构造几乎最优的二分字典树，还有最优二分分裂树（Optimal Binary Split Tree）、偏斜二分搜索树（Biased Search Tree）、自调节二分搜索树(Self-adjusting BST）等等。 267.3 近似串匹配问题（Approximate String Matching）7.3.1 近似串匹配问题 设样本为P：p1p2.pm，文本为T：t1t2.tn。K-近似匹配（K-approximate match）：对于非负整数K，样

23、本P在文本T中的K-近似匹配，是指P在T中的包含至多K个差别的匹配。这里所指的差别（difference）是指下列三种情况之一： 修改（revise）：P与T中的对应字符不同； 删去（delete）：T中含有一个未出现在P中的字符； 插入（insert）：T中不含有出现在P中的一个字符。27例如，下面是一个包含上述三种差别的3-近似匹配，一般的情形是样本P是正确的，文本T中包含上述三类错误，一般称之为编辑错误。事实上，能够指出上例中的两个字符串有三个差别，并不是一件容易的事，因为不同的对应方法可以得到不同的K值。例如，把两者从字符a开始，顺序一一对应，可以计算出有6个修改（revise）错误。

24、改变对应方法，就会产生不同的结论。因此应该指出，P与T为K-近似匹配，包含下面两层含义： 1 二者的差别数至多为K； 2 差别数是指二者在所有不同匹配对应方式下的最小编辑错误总数。P:a p p r o x i m a t l yT:a p r o x i o m a l l y28因此，一般的K-近似匹配问题可以描述为：已知：字符串P： p1p2.pm（称为Pattern），字符串T：t1t2.tn（称为Text），和一个正整数K。求：样式P在文本T上的K-近似匹配的第一次出现或所有出现。在完全串匹配中，P在T上的一次匹配检查十分简单，用至多m次字符比较即可完成，问题的关键是确定一次匹配检查

25、之后样本P的移动方法。最简单的方法是移动一位，更巧妙的方法则可以在不丢解的条件下移动较大的距离。K-近似匹配的关键是P在T上的K-近似匹配检查，即计算样本P与当前对应的文本T之子串的最小差别数。如果该值小于等于K，则找到结果，否则样本P右移。过程中的计算集中在计算最小差别数上，这实际上是一个优化问题，比完全匹配问题复杂得多。29近似串匹配的一个简化问题是比较两个字符串的差别。例如，OCR（光学字符识别）系统是通过在计算机上运行识别算法，将扫描仪扫描得到的字符文本的图像信息转化为作为识别结果的字符串。为了确定OCR系统的识别率，需要对原始字符串和识别结果字符串进行比较。这种比较，实际上就是简化的

26、近似串匹配问题，它应计算出结果字符串与原始字符串间的最小编辑错误数。7.3.2 DP算法的思路 近似串匹配问题同样具有最优子结构性质和子结构重迭性质。 如果样本p1p2.pm在文本T的某一位置上有一最优（差别数最小）的对应关系，那么样本P的任意一个子串pi.pj（1ijm）与T的对应关系也必然是最优（差别数最小）的。(最优子结构)计算样本P对应文本T某一位置的最小差别数时，对P的任一子串pi.pj与T相应位置上的最小差别数的计算，必然包含计算任何包含pi.pj的P及其子串的最小差别数的过程中。(子结构覆盖)30因此，可以考虑用动态规划方法设计出快速的近似串匹配算法。问题的关键是如何找出代价函数

27、的自底向上的递推关系。为此可定义一个代价函数Dij，0im，0jn。Dij表示样本子串p1.pi与文本子串t1.tj之间的最小差别数。由此可知：Dmj表示样本P在文本T的位置j处的最小差别数，如果DmjK，说明P在tj处找到了K-近似匹配。代价函数的初始值很容易确定： D0j=0，这是因为样本P为空串，与文本t1.tj有0个字符不同；Di0=i，这是因为样本串p1.pi与空文本串相比，有i个字符不同。31可以通过下面的简单分析得出代价函数的递推关系。当样本子串p1.pi与文本子串t1.tj对应时，有四种可能的情况： (1) 字符pi与tj相对应且pi=tj，总差别数为Di-1j-1； (2)

28、字符pi与tj相对应且pitj，总差别数为Di-1j-1+1； (3) 字符tj为多余，即tj对应于空格，总差别数为Dij-1+1； (4) 字符pi对应于tj后的空格，总差别数为Di-1j+1。根据代价函数Dij的定义，它应该取所有可能值之中的最小者，即： 32代价函数Dij实际上是一个二维数组，其初始值和Dij的计算方法由Fig7.7表示，每个Dij的值总是由其左上方的三个已知值来确定。 337.3.3 DP算法算法7.4 近似串匹配算法 ASM 7.3.4 算法的复杂度分析与实例在算法7.4的内层循环语句中，进行了一次字符比较和取三个整数中最小值的操作，因此，算法的时间复杂度为T(n)=

29、O(nm)。算法的空间代价显然由代价函数矩阵Dmn决定，因此也是O(nm)。与近似串匹配问题本身的复杂性相对照，动态规划的解决算法十分简明和有效，这一点可以通过下面的实例看出。已知：样本P=happy，K=1，T=Have a hsppy day.是一个可能有编辑错误的文本。求：P在T中的K-近似出现。34按照算法7.4解这一问题的计算过程，实际上就是逐列计算二维数组D1.m1.n的过程，当在第m=5行出现小于或等于K的值时，计算即可终止。这个过程可用Fig7.8中给出的二维数组表示。 35在计算过程中，首先D0j全部置为0，Di0置为i。由于认为大写的“H”与小写“h”是不同的，因此D11在

30、D00+1=1、D01+1=1、D10+1=2三者中取最小值，故D11=1。在逐列的计算过程中，总是根据Dij对应的两个字符pi、tj是否相同，以及位于Dij左上方、正上方、正左方的三个已计算值决定Dij的取值。例如：计算D26时，因为样本与文本的对应字符都为“a”，所以应在1、1+1=2、2+1=3三个值中取最小者，故D26=1。计算D412时，样本字符为“p”，文本字符为“y”，这时三个相关的已知值为2、3、1，所以应在3、4、2中取最小值，故D412=2。36由于D512=1=K，且m=5，所以在T12处找到了差别数为1的K-近似匹配。这时文本T和样本P的对应关系为：T： Have a

31、hsppyP： happy事实上，当K=2时，在T11、T12、T13三处都可以找到2-近似匹配。7.3.5 讨论1. 近似串匹配问题由于允许有三种类型的差别，两个字符串的对应关系可以有很多种选择，不同的选择方法总数为样本长度m的指数函数。如果采用穷举算法，计算各种不同对应关系的差别数，再求最小值，计算量是很大的。然而，通过使用动态规划方法，使得计算复杂度仅为O(nm)，十分简捷。372. 算法7.4只计算出样本P在文本T各个位置的最小差别数，以及相应的位置Ti，没有给出得到最小差别的字符对应关系。做到这一点并不困难，可以再设置一个mn阶的二维数组B1.m1.n，Bij的值可以在计算Dij的过

32、程中求得。在计算最小值时，Bij的值依据Di-1j-1、Dij-1、Di-1j而分别取为r（revise）、i（insert）、d（delete）。当Dmj=K时，即可由Bmj回溯，求出这一匹配的字符对应关系。3. 近似串匹配问题可以定义为求样本P在文本T中所有的K-近似匹配，实现方法是把二维数组Dnm的所有值全计算出来即可。这只需对算法7.4稍作修改：将原程序中的“if(Dmj=K) return j;”改为“if(Dmj=K) coutK=Dmj:jendl;”在Fig7.8中，虚线右侧的值就是由修改后的程序计算出来的。384. 另一个很有实用价值的近似串匹配问题，是比较两个文本文件f1和

33、f2的差别。例如，检验电脑录入人员的文字输入的准确率，总是给出一个标准文本f，请录入员在指定时间内输入得到一个结果文件f，由计算机判定其出错数。实际上就是计算字符串f与标准串f之间的差别数。这个问题并不简单，它是近似串匹配问题的一个变种，其中两个字符串的地位是平等的，其算法应与算法7.4有所不同。5. 字符串的近似匹配问题可以看作最简单的模式识别问题。当把字符串扩展为数据串时，可以发现，语音识别问题、在线（on-line）手写字符的识别问题与近似串匹配问题在本质上有很大的相似性，因此，动态规划算法在语音识别和手写字符识别领域可以发挥很大的作用。笔者曾经设计一种二维的动态规划算法用以解决平面手写

34、字符的脱机弹性匹配问题中，也取得了成功。第七章完3940递归函数fibo int fibo(int n) int f; if(n2) f=n; else f=fibo(n-1)+fibo(n-2); return f;返回41函数fibo的改进函数fib int fib(int n) int i,f,f0=0,f1=1; for(i=2;i=0;i-) for(j=i+1;j=n;j+) if(j-i=1) mcost=0; mlast=-1; else mcost=maxint;返回下页44算法7.1 最优矩阵连乘算法MatrixOrder（续） for(k=i+1;kj;k+) a=costik; b=costkj; c=dimi*dimk*dimj; if( (a+b+c)1) k=lastlowhigh; ExtractOrder(low,k,l