3.1 算法思想
和贪婪算法一样,在动态规划中,可将一个问题的解决方案视为一系列决策的结果。不同的是,在贪婪算法中,每采用一次贪婪准则便做出一个不可撤回的决策,而在动态规划中,还要考察每个最优决策序列中是否包含一个最优子序列。
例3-1 [最短路经] 考察图1 2 - 2中的有向图。假设要寻找一条从源节点s= 1到目的节点d= 5的最短路径,即选择此路径所经过的各个节点。第一步可选择节点2,3或4。假设选择了节点3,则此时所要求解的问题变成:选择一条从3到5的最短路径。如果3到5的路径不是最短的,则从1开始经过3和5的路径也不会是最短的。例如,若选择的子路径(非最短路径)是3,2,5 (耗费为9 ),则1到5的路径为1,3,2,5 (耗费为11 ),这比选择最短子路径3,4,5而得到的1到5的路径1,3,4,5 (耗费为9) 耗费更大。
所以在最短路径问题中,假如在的第一次决策时到达了某个节点v,那么不管v 是怎样确定的,此后选择从v 到d 的路径时,都必须采用最优策略。
例3-2 [0/1背包问题] 考察1 3 . 4节的0 / 1背包问题。如前所述,在该问题中需要决定x1 .. xn的值。假设按i = 1,2,.,n 的次序来确定xi 的值。如果置x1 = 0,则问题转变为相对于其余物品(即物品2,3,.,n),背包容量仍为c 的背包问题。若置x1 = 1,问题就变为关于最大背包容量为c-w1 的问题。现设r?{c,c-w1 } 为剩余的背包容量。
在第一次决策之后,剩下的问题便是考虑背包容量为r 时的决策。不管x1 是0或是1,[x2 ,.,xn ] 必须是第一次决策之后的一个最优方案,如果不是,则会有一个更好的方案[y2,.,yn ],因而[x1,y2,.,yn ]是一个更好的方案。
假设n=3, w=[100,14,10], p=[20,18,15], c= 11 6。若设x1 = 1,则在本次决策之后,可用的背包容量为r= 116-100=16 。[x2,x3 ]=[0,1] 符合容量限制的条件,所得值为1 5,但因为[x2,x3 ]= [1,0] 同样符合容量条件且所得值为1 8,因此[x2,x3 ] = [ 0,1] 并非最优策略。即x= [ 1,0,1] 可改进为x= [ 1,1,0 ]。若设x1 = 0,则对于剩下的两种物品而言,容量限制条件为11 6。总之,如果子问题的结果[x2,x3 ]不是剩余情况下的一个最优解,则[x1,x2,x3 ]也不会是总体的最优解。
例3-3 [航费] 某航线价格表为:从亚特兰大到纽约或芝加哥,或从洛杉矶到亚特兰大的费用为$ 1 0 0;从芝加哥到纽约票价$ 2 0;而对于路经亚特兰大的旅客,从亚特兰大到芝加哥的费用仅为$ 2 0。从洛杉矶到纽约的航线涉及到对中转机场的选择。如果问题状态的形式为(起点,终点),那么在选择从洛杉矶到亚特兰大后,问题的状态变为(亚特兰大,纽约)。从亚特兰大到纽约的最便宜航线是从亚特兰大直飞纽约,票价$ 1 0 0。而使用直飞方式时,从洛杉矶到纽约的花费为$ 2 0 0。不过,从洛杉矶到纽约的最便宜航线为洛杉矶-亚特兰大-芝加哥-纽约,其总花费为$ 1 4 0(在处理局部最优路径亚特兰大到纽约过程中选择了最低花费的路径:亚特兰大-芝加哥-纽约)。
如果用三维数组(t a g,起点,终点)表示问题状态,其中t a g为0表示转飞, t a g为1表示其他情形,那么在到达亚特兰大后,状态的三维数组将变为( 0,亚特兰大,纽约),它对应的最优路径是经由芝加哥的那条路径。
当最优决策序列中包含最优决策子序列时,可建立动态规划递归方程( d y n a m i c -programming recurrence equation),它可以帮助我们高效地解决问题。
例3-4 [0/1背包] 在例3 - 2的0 / 1背包问题中,最优决策序列由最优决策子序列组成。假设f (i,y) 表示例1 5 - 2中剩余容量为y,剩余物品为i,i + 1,.,n 时的最优解的值,即:和利用最优序列由最优子序列构成的结论,可得到f 的递归式。f ( 1 ,c) 是初始时背包问题的最优解。可使用( 1 5 - 2)式通过递归或迭代来求解f ( 1 ,c)。从f (n, * )开始迭式, f (n, * )由(1 5 - 1)式得出,然后由( 1 5 - 2)式递归计算f (i,*) ( i=n- 1,n- 2,., 2 ),最后由( 1 5 - 2)式得出f ( 1 ,c)。
对于例1 5 - 2,若0≤y<1 0,则f ( 3 ,y) = 0;若y≥1 0,f ( 3 ,y) = 1 5。利用递归式(1 5 - 2),可得f (2, y) = 0 ( 0≤y<10 );f(2,y)= 1 5(1 0≤y<1 4);f(2,y)= 1 8(1 4≤y<2 4)和f(2,y)= 3 3(y≥2 4)。因此最优解f ( 1 , 11 6 ) = m a x {f(2,11 6),f(2,11 6 - w1)+ p1} = m a x {f(2,11 6),f(2,1 6)+ 2 0 } = m a x { 3 3,3 8 } = 3 8。
现在计算xi 值,步骤如下:若f ( 1 ,c) =f ( 2 ,c),则x1 = 0,否则x1 = 1。接下来需从剩余容量c-w1中寻求最优解,用f (2, c-w1) 表示最优解。依此类推,可得到所有的xi (i= 1.n) 值。
在该例中,可得出f ( 2 , 11 6 ) = 3 3≠f ( 1 , 11 6 ),所以x1 = 1。接着利用返回值3 8 -p1=18 计算x2 及x3,此时r = 11 6 -w1 = 1 6,又由f ( 2 , 1 6 ) = 1 8,得f ( 3 , 1 6 ) = 1 4≠f ( 2 , 1 6 ),因此x2 = 1,此时r= 1 6 -w2 = 2,所以f (3,2) =0,即得x3 = 0。
动态规划方法采用最优原则( principle of optimality)来建立用于计算最优解的递归式。所谓最优原则即不管前面的策略如何,此后的决策必须是基于当前状态(由上一次决策产生)的最优决策。由于对于有些问题的某些递归式来说并不一定能保证最优原则,因此在求解问题时有必要对它进行验证。若不能保持最优原则,则不可应用动态规划方法。在得到最优解的递归式之后,需要执行回溯(t r a c e b a c k)以构造最优解。
编写一个简单的递归程序来求解动态规划递归方程是一件很诱人的事。然而,正如我们将在下文看到的,如果不努力地去避免重复计算,递归程序的复杂性将非常可观。如果在递归程序设计中解决了重复计算问题时,复杂性将急剧下降。动态规划递归方程也可用迭代方式来求解,这时很自然地避免了重复计算。尽管迭代程序与避免重复计算的递归程序有相同的复杂性,但迭代程序不需要附加的递归栈空间,因此将比避免重复计算的递归程序更快。
3.2 应用
3.2.1 0/1背包问题
1. 递归策略
在例3 - 4中已建立了背包问题的动态规划递归方程,求解递归式( 1 5 - 2)的一个很自然的方法便是使用程序1 5 - 1中的递归算法。该模块假设p、w 和n 为输入,且p 为整型,F(1,c) 返回f ( 1 ,c) 值。
程序15-1 背包问题的递归函数
int F(int i, int y)
{// 返回f ( i , y ) .
if (i == n) return (y < w[n]) ? 0 : p[n];
if (y < w[i]) return F(i+1,y);
return max(F(i+1,y), F(i+1,y-w[i]) + p[i]);
}
程序1 5 - 1的时间复杂性t (n)满足:t ( 1 ) =a;t(n)≤2t(n- 1)+b(n>1),其中a、b 为常数。通过求解可得t (n) =O( 2n)。
例3-5 设n= 5,p= [ 6 , 3 , 5 , 4 , 6 ],w=[2,2,6,5,4] 且c= 1 0 ,求f ( 1 , 1 0 )。为了确定f ( 1 , 1 0 ),调用函数F ( 1 , 1 0 )。递归调用的关系如图1 5 - 1的树型结构所示。每个节点用y值来标记。对于第j层的节点有i=j,因此根节点表示F ( 1 , 1 0 ),而它有左孩子和右孩子,分别对应F ( 2 , 1 0 )和F ( 2 , 8 )。总共执行了2 8次递归调用。但我们注意到,其中可能含有重复前面工作的节点,如f ( 3 , 8 )计算过两次,相同情况的还有f ( 4 , 8 )、f ( 4 , 6 )、f ( 4 , 2 )、f ( 5 , 8 )、f ( 5 , 6 )、f ( 5 , 3 )、f (5,2) 和f ( 5 , 1 )。如果保留以前的计算结果,则可将节点数减至1 9,因为可以丢弃图中的阴影节点。
正如在例3 - 5中所看到的,程序1 5 - 1做了一些不必要的工作。为了避免f (i,y)的重复计算,必须定义一个用于保留已被计算出的f (i,y)值的表格L,该表格的元素是三元组(i,y,f (i,y) )。在计算每一个f (i,y)之前,应检查表L中是否已包含一个三元组(i,y, * ),其中*表示任意值。如果已包含,则从该表中取出f (i,y)的值,否则,对f (i,y)进行计算并将计算所得的三元组(i,y,f (i,y) )加入表L。L既可以用散列(见7 . 4节)的形式存储,也可用二叉搜索树(见11章)的形式存储。
2. 权为整数的迭代方法
当权为整数时,可设计一个相当简单的算法(见程序1 5 - 2)来求解f ( 1 ,c)。该算法基于例3 - 4所给出的策略,因此每个f (i,y) 只计算一次。程序1 5 - 2用二维数组f [ ][ ]来保存各f 的值。而回溯函数Tr a c e b a c k用于确定由程序1 5 - 2所产生的xi 值。函数K n a p s a c k的复杂性为( n c),而Tr a c e b a c k的复杂性为( n )。
程序15-2 f 和x 的迭代计算
template<class T>
void Knapsack(T p[], int w[], int c, int n, T** f)
{// 对于所有i和y计算f [ i ] [ y ]
// 初始化f [ n ] [ ]
for (int y = 0; y <= yMax; y++)
f[n][y] = 0;
for (int y = w[n]; y <= c; y++)
f[n][y] = p[n];
// 计算剩下的f
for (int i = n - 1; i > 1; i--) {
for (int y = 0; y <= yMax; y++)
f[i][y] = f[i+1][y];
for (int y = w[i]; y <= c; y++)
f[i][y] = max(f[i+1][y], f[i+1][y-w[i]] + p[i]);
}
f[1][c] = f[2][c];
if (c >= w[1])
f[1][c] = max(f[1][c], f[2][c-w[1]] + p[1]);
}
template<class T>
void Traceback(T **f, int w[], int c, int n, int x[])
{// 计算x
for (int i = 1; i < n; i++)
if (f[i][c] == f[i+1][c]) x[i] = 0;
else {x[i] = 1;
c -= w[i];}
x[n] = (f[n][c]) ? 1 : 0;
}
3. 元组方法( 选读)
程序1 5 - 2有两个缺点:1) 要求权为整数;2) 当背包容量c 很大时,程序1 5 - 2的速度慢于程序1 5 - 1。一般情况下,若c>2n,程序1 5 - 2的复杂性为W (n2n )。可利用元组的方法来克服上述两个缺点。在元组方法中,对于每个i,f (i, y) 都以数对(y, f (i, y)) 的形式按y的递增次序存储于表P(i)中。同时,由于f (i, y) 是y 的非递减函数,因此P(i) 中各数对(y, f (i, y)) 也是按f (i, y) 的递增次序排列的。
例3-6 条件同例3 - 5。对f 的计算如图1 5 - 2所示。当i= 5时,f 由数对集合P( 5 ) = [ ( 0 , 0 ) , ( 4 , 6 ) ]表示。而P( 4 )、P( 3 )和P( 2 )分别为[ ( 0 , 0 ) , ( 4 , 6 ) , ( 9 , 1 0 ) ]、[ ( 0 , 0 ) ( 4 , 6 ) , ( 9 , 1 0 ) , ( 1 0 , 11)] 和[ ( 0 , 0 ) ( 2 , 3 ) ( 4 , 6 ) ( 6 , 9 ) ( 9 , 1 0 ) ( 1 0 , 11 ) ]。
为求f ( 1 , 1 0 ),利用式(1 5 - 2)得f ( 1 , 1 0 ) = m a x{f ( 2 , 1 0 ),f ( 2 , 8 ) + p 1}。由P( 2 )得f ( 2 , 1 0 ) = 11、f (2,8)=9 (f ( 2 , 8 ) = 9来自数对( 6,9 ) ),因此f ( 1 , 1 0 ) = m a x{11 , 1 5}= 1 5。现在来求xi 的值,因为f ( 1 , 1 0 ) =f ( 2 , 6 ) +p1,所以x1 = 1;由f ( 2 , 6 ) =f ( 3 , 6 - w 2 ) +p2 =f ( 3 , 4 ) +p2,得x2 = 1;由f ( 3 , 4 ) =f ( 4 , 4 ) =f ( 5 , 4 )得x3=x4 = 0;最后,因f ( 5 , 4 )≠0得x5= 1。
检查每个P(i) 中的数对,可以发现每对(y,f (i,y)) 对应于变量xi , ., xn 的0/1 赋值的不同组合。设(a,b)和(c,d)是对应于两组不同xi , ., xn 的0 / 1赋值,若a≥c且b<d,则(a, b) 受(b, c) 支配。被支配者不必加入P(i)中。若在相同的数对中有两个或更多的赋值,则只有一个放入P(i)。假设wn≤C,P(n)=[(0,0), (wn , pn ) ],P(n)中对应于xn 的两个数对分别等于0和1。对于每个i,P(i)可由P(i+ 1 )得出。首先,要计算数对的有序集合Q,使得当且仅当wi≤s≤c且(s-wi , t-pi )为P(i+1) 中的一个数对时,(s,t)为Q中的一个数对。现在Q中包含xi = 1时的数对集,而P(i+ 1 )对应于xi = 0的数对集。接下来,合并Q和P(i+ 1 )并删除受支配者和重复值即可得到P(i)。
例3-7 各数据同例1 5 - 6。P(5)=[(0,0),(4,6)], 因此Q= [ ( 5 , 4 ) , ( 9 , 1 0 ) ]。现在要将P( 5 )和Q合并得到P( 4 )。因( 5 , 4 )受( 4 , 6 )支配,可删除( 5 , 4 ),所以P(4)=[(0,0), (4,6), (9,10)]。接着计算P( 3 ),首先由P( 4 )得Q=[(6,5), (10,11 ) ],然后又由合并方法得P(3)=[(0,0), (4,6), (9,10), (10,11 ) ]。最后计算P( 2 ):由P( 3 )得Q= [ ( 2 , 3 ),( 6 , 9 ) ],P( 3 )与Q合并得P(2)=[(0,0), (2,3), (4,6), (6,9), (9,10). (10,11 ) ]。因为每个P(i) 中的数对对应于xi , ., xn 的不同0 / 1赋值,因此P(i) 中的数对不会超过2n-i+ 1个。计算P(i) 时,计算Q需消耗( |P(i+ 1 ) |)的时间,合并P(i+1) 和Q同样需要( |P(i+ 1 ) | )的时间。计算所有P(i) 时所需要的总时间为: (n ?i=2|P(i + 1)|= O ( 2n )。当权为整数时,|P(i) |≤c+1, 此时复杂性为O ( m i n {n c, 2n } )。
如6 . 4 . 3节定义的,数字化图像是m×m的像素阵列。假定每个像素有一个0 ~ 2 5 5的灰度值。因此存储一个像素至多需8位。若每个像素存储都用最大位8位,则总的存储空间为8m2 位。为了减少存储空间,我们将采用变长模式( variable bit scheme),即不同像素用不同位数来存储。像素值为0和1时只需1位存储空间;值2、3各需2位;值4,5,6和7各需3位;以此类推,使用变长模式的步骤如下:
1) 图像线性化根据图15-3a 中的折线将m×m维图像转换为1×m2 维矩阵。
2) 分段将像素组分成若干个段,分段原则是:每段中的像素位数相同。每个段是相邻像素的集合且每段最多含2 5 6个像素,因此,若相同位数的像素超过2 5 6个的话,则用两个以上的段表示。
3) 创建文件创建三个文件:S e g m e n t L e n g t h, BitsPerPixel 和P i x e l s。第一个文件包含在2 )中所建的段的长度(减1 ),文件中各项均为8位长。文件BitsPerPixel 给出了各段中每个像素的存储位数(减1),文件中各项均为3位。文件Pixels 则是以变长格式存储的像素的二进制串。
4) 压缩文件压缩在3) 中所建立的文件,以减少空间需求。
上述压缩方法的效率(用所得压缩率表示)很大程度上取决于长段的出现频率。
例3-8 考察图15-3b 的4×4图像。按照蛇形的行主次序,灰度值依次为1 0,9,1 2,4 0,5 0,3 5,1 5,1 2,8,1 0,9,1 5,11,1 3 0,1 6 0和2 4 0。各像素所需的位数分别为4,4,4,6,6,6,4,4,4,4,4,4,4,8,8和8,按等长的条件将像素分段,可以得到4个段[ 1 0,9,1 2 ]、[ 4 0,5 0,3 5 ]、[15, 12, 8, 10, 9, 15, 11] 和[130, 160, 240]。因此,文件SegmentLength 为2,2,6,2;文件BitsPerSegment 的内容为3,5,3,7;文件P i x e l s包含了按蛇形行主次序排列的1 6个灰度值,其中头三个各用4位存储,接下来三个各用6位,再接下来的七个各用4位,最后三个各用8位存储。因此存储单元中前3 0位存储了前六个像素:
1010 1001 1100 111000 110010 100011
这三个文件需要的存储空间分别为:文件SegmentLength 需3 2位;BitsPerSegment 需1 2位;Pixels 需8 2位,共需1 2 6位。而如果每个像素都用8位存储,则存储空间需8×1 6 = 1 2 8位,因而在本例图像中,节省了2位的空间。
假设在2) 之后,产生了n 个段。段标题(segment header)用于存储段的长度以及该段中每个像素所占用的位数。每个段标题需11位。现假设li 和bi 分别表示第i 段的段长和该段每个像素的长度,则存储第i 段像素所需要的空间为li *bi 。在2) 中所得的三个文件的总存储空间为11 n+n ?i = 1li bi。可通过将某些相邻段合并的方式来减少空间消耗。如当段i 和i+ 1被合并时,合并后的段长应为li +li + 1。此时每个像素的存储位数为m a x {bi,bi +1 } 位。尽管这种技术增加了文件P i x e l s的空间消耗,但同时也减少了一个段标题的空间。
例3-9 如果将例1 5 - 8中的第1段和第2段合并,合并后,文件S e g m e n t L e n g t h变为5,6,2,BitsPerSegment 变为5,3,7。而文件Pixels 的前3 6位存储的是合并后的第一段:001010 001001 001100 111000 110010 100011其余的像素(例1 5 - 8第3段)没有改变。因为减少了1个段标题,文件S e g m e n t L e n g t h和BitsPerPixel 的空间消耗共减少了11位,而文件Pixels 的空间增加6位,因此总共节约的空间为5位,空间总消耗为1 2 1位。
我们希望能设计一种算法,使得在产生n 个段之后,能对相邻段进行合并,以便产生一个具有最小空间需求的新的段集合。在合并相邻段之后,可利用诸如L Z W法(见7 . 5节)和霍夫曼编码(见9 . 5 . 3节)等其他技术来进一步压缩这三个文件。
令sq 为前q 个段的最优合并所需要的空间。定义s0 = 0。考虑第i 段(i>0 ),假如在最优合并C中,第i 段与第i- 1,i- 2,.,i-r+1 段相合并,而不包括第i-r 段。合并C所需要的空间消耗等于:第1段到第i-r 段所需空间+ l s u m (i-r+ 1 ,i) * b m a x (i-r+ 1 ,i) + 11
其中l s u m(a, b)=b ?j =a
lj, bmax (a, b)= m a x {ba , ..., bb }。假如在C中第1段到第i-r 段的合并不是最优合并,那么需要对合并进行修改,以使其具有更小的空间需求。因此还必须对段1到段i-r 进行最优合并,也即保证最优原则得以维持。故C的空间消耗为:
si = si-r +l s u m(i-r+1, i)*b m a x(i-r+1, i)+ 11
r 的值介于1到i 之间,其中要求l s u m不超过2 5 6 (因为段长限制在2 5 6之内)。尽管我们不知道如何选择r,但我们知道,由于C具有最小的空间需求,因此在所有选择中, r 必须产生最小的空间需求。
假定k a yi 表示取得最小值时k 的值,sn 为n 段的最优合并所需要的空间,因而一个最优合并可用kay 的值构造出来。
例3-10 假定在2) 中得到五个段,它们的长度为[ 6,3,1 0,2,3 ],像素位数为[ 1,2,3,2,1 ],要用公式(1 5 - 3)计算sn,必须先求出sn-1,.,s0 的值。s0 为0,现计算s1:s1 =s0 +l1 *b1+ 11 = 1 7k a y1 = 1s2 由下式得出:
s2 = m i n {s1 +l2 b2 , s0 + (l1 +l2 ) * m a x {b1 , b2} } + 11 = m i n { 1 7 + 6 , 0 + 9 * 2 } + 11 = 2 9
k a y2 = 2
以此类推,可得s1.s5 = [ 1 7,2 9,6 7,7 3,82] ,k a y1.k a y5 = [ 1,2,2,3,4 ]。因为s5 = 8 2,所以最优空间合并需8 2位的空间。可由k a y5 导出本合并的方式,过程如下:因为k a y5 = 4,所以s5 是由公式(1 5 - 3)在k=4 时取得的,因而最优合并包括:段1到段( 5 - 4 ) = 1的最优合并以及段2,3,4和5的合并。最后仅剩下两个段:段1以及段2到段5的合并段。
1. 递归方法
用递归式(1 5 - 3)可以递归地算出si 和k a yi。程序1 5 - 3为递归式的计算代码。l,b,和k a y是一维的全局整型数组,L是段长限制( 2 5 6),h e a d e r为段标题所需的空间( 11 )。调用S ( n )返回sn 的值且同时得出k a y值。调用Tr a c e b a c k ( k a y, n )可得到最优合并。
现讨论程序1 5 - 3的复杂性。t( 0 ) =c(c 为一个常数): (n>0),因此利用递归的方法可得t (n) = O ( 2n )。Tr a c e b a c k的复杂性为(n)。
程序15-3 递归计算s , k a y及最优合并
int S(int i)
{ / /返回S ( i )并计算k a y [ i ]
if (i == 0 ) return 0;
//k = 1时, 根据公式( 1 5 - 3)计算最小值
int lsum = l[i],bmax = b[i];
int s = S(i-1) + lsum * bmax;
kay[i] = 1;
/ /对其余的k计算最小值并求取最小值
for (int k = 2; k <= i && lsum+l[i-k+1] <= L; k++) {
lsum += l[i-k+1];
if (bmax < b[i-k+1]) bmax = b[i-k+1];
int t = S(i-k);
if (s > t + lsum * bmax) {
s = t + lsum * bmax;
kay[i] = k;}
}
return s + header;
}
void Traceback(int kay[], int n)
{// 合并段
if (n == 0) return;
Tr a c e b a c k ( k a y, n-kay[n]);
cout << "New segment begins at " << (n - kay[n] + 1) << endl;
}
2. 无重复计算的递归方法
通过避免重复计算si,可将函数S的复杂性减少到(n)。注意这里只有n个不同的si。
例3 - 11 再考察例1 5 - 1 0中五个段的例子。当计算s5 时,先通过递归调用来计算s4,.,s0。计算s4 时,通过递归调用计算s3,.,s0,因此s4 只计算了一次,而s3 计算了两次,每一次计算s3要计算一次s2,因此s2 共计算了四次,而s1 重复计算了1 6次!可利用一个数组s 来保存先前计算过的si 以避免重复计算。改进后的代码见程序1 5 - 4,其中s为初值为0的全局整型数组。
程序15-4 避免重复计算的递归算法
int S(int i)
{ / /计算S ( i )和k a y [ i ]
/ /避免重复计算
if (i == 0) return 0;
if (s[i] > 0) return s[i]; //已计算完
/ /计算s [ i ]
/ /首先根据公式(1 5 - 3)计算k = 1时最小值
int lsum = l[i], bmax = b[i];
s[i] =S(i-1) + lsum * bmax;
kay[i] = 1;
/ /对其余的k计算最小值并更新
for (int k = 2; k <= i && lsum+l[i-k+1] <= L; k++) {
lsum += l[i-k+1];
if (bmax < b[i-k+1]) bmax = b[i-k+1];
int t = S(i-k);
if (s[i] > t + lsum * bmax) {
s[i] = t + lsum * bmax;
kay[i] = k;}
}
s[i] += header;
return s[i];
}
为了确定程序1 5 - 4的时间复杂性,我们将使用分期计算模式( amortization scheme)。在该模式中,总时间被分解为若干个不同项,通过计算各项的时间然后求和来获得总时间。当计算si 时,若sj 还未算出,则把调用S(j) 的消耗计入sj ;若sj 已算出,则把S(j) 的消耗计入si (这里sj依次把计算新sq 的消耗转移至每个sq )。程序1 5 - 4的其他消耗也被计入si。因为L是2 5 6之内的常数且每个li 至少为1,所以程序1 5 - 4的其他消耗为( 1 ),即计入每个si 的量是一个常数,且si 数目为n,因而总工作量为(n)。
3. 迭代方法
倘若用式(1 5 - 3)依序计算s1 , ., sn,便可得到一个复杂性为(n)的迭代方法。在该方法中,在si 计算之前, sj 必须已计算好。该方法的代码见程序1 5 - 5,其中仍利用函数Tr a c e b a c k(见程序1 5 - 3)来获得最优合并。
程序15-5 迭代计算s和k a y
void Vbits (int l[], int b[], int n, int s[], int kay[])
{ / /计算s [ i ]和k a y [ i ]
int L = 256, header = 11 ;
s[0] = 0;
/ /根据式(1 5 - 3)计算s [ i ]
for (int i = 1; i <= n; i++) {
// k = 1时,计算最小值
int lsum = l,
bmax = b[i];
s[i] = s[i-1] + lsum * bmax;
kay[i] = 1;
/ /对其余的k计算最小值并更新
for (int k=2; k<= i && lsum+l[i-k+1]<= L; k++) {
lsum += l[i-k+1];
if (bmax < b[i-k+1]) bmax = b[i-k+1];
if (s[i] > s[i-k] + lsum * bmax){
s[i] = s[i-k] + lsum * bmax;
kay[i] = k; }
}
s[i] += header;
}
}
3.2.3 矩阵乘法链
m× n矩阵A与n×p矩阵B相乘需耗费(m n p)的时间(见第2章练习1 6)。我们把m n p作为两个矩阵相乘所需时间的测量值。现假定要计算三个矩阵A、B和C的乘积,有两种方式计算此乘积。在第一种方式中,先用A乘以B得到矩阵D,然后D乘以C得到最终结果,这种乘法的顺序可写为(A*B) *C。第二种方式写为A* (B*C) ,道理同上。尽管这两种不同的计算顺序所得的结果相同,但时间消耗会有很大的差距。
例3-12 假定A为1 0 0×1矩阵,B为1×1 0 0矩阵,C为1 0 0×1矩阵,则A*B的时间耗费为10 0 0 0,得到的结果D为1 0 0×1 0 0矩阵,再与C相乘所需的时间耗费为1 000 000,因此计算(A*B) *C的总时间为1 010 000。B*C的时间耗费为10 000,得到的中间矩阵为1×1矩阵,再与A相乘的时间消耗为1 0 0,因而计算A*(B*C)的时间耗费竟只有10 100!而且,计算( A*B)*C时,还需10 000个单元来存储A*B,而A*(B*C)计算过程中,只需用1个单元来存储B*C。
下面举一个得益于选择合适秩序计算A*B*C矩阵的实例:考虑两个3维图像的匹配。图像匹配问题的要求是,确定一个图像需旋转、平移和缩放多少次才能逼近另一个图像。实现匹配的方法之一便是执行约1 0 0次迭代计算,每次迭代需计算1 2×1个向量T:
T=?A(x, y, z) *B(x, y, z)*C(x, y, z )
其中A,B和C分别为1 2×3,3×3和3×1矩阵。(x , y, z) 为矩阵中向量的坐标。设t 表示计算A(x , y, z) *B(x , y, z) *C(x , y, z)的计算量。假定此图像含2 5 6×2 5 6×2 5 6个向量,在此条件中,这1 0 0个迭代所需的总计算量近似为1 0 0 * 2 5 63 * t≈1 . 7 * 1 09 t。若三个矩阵是按由左向右的顺序相乘的,则t = 1 2 * 3 * 3 + 1 2 * 3 *1= 1 4 4;但如果从右向左相乘, t = 3 * 3 * 1 + 1 2 * 3 * 1 = 4 5。由左至右计算约需2 . 4 * 1 011个操作,而由右至左计算大概只需7 . 5 * 1 01 0个操作。假如使用一个每秒可执行1亿次操作的计算机,由左至右需4 0分钟,而由右至左只需1 2 . 5分钟。
在计算矩阵运算 A*B*C时,仅有两种乘法顺序(由左至右或由右至左),所以可以很容易算出每种顺序所需要的操作数,并选择操作数比较少的那种乘法顺序。但对于更多矩阵相乘来说,情况要复杂得多。如计算矩阵乘积M1×M2×.×Mq,其中Mi 是一个ri×ri + 1 矩阵( 1≤i≤q)。不妨考虑q=4 的情况,此时矩阵运算A*B*C*D可按以下方式(顺序)计算:
A* ( (B*C) *D) A* (B* (C*D)) (A*B) * (C*D) (A* (B*C) ) *D
不难看出计算的方法数会随q 以指数级增加。因此,对于很大的q 来说,考虑每一种计算顺序并选择最优者已是不切实际的。
现在要介绍一种采用动态规划方法获得矩阵乘法次序的最优策略。这种方法可将算法的时间消耗降为(q3 )。用Mi j 表示链Mi×.×Mj (i≤j)的乘积。设c(i,j) 为用最优法计算Mi j 的消耗,k a y(i, j) 为用最优法计算Mi j 的最后一步Mi k×Mk+1, j 的消耗。因此Mij 的最优算法包括如何用最优算法计算Mik 和Mkj 以及计算Mik×Mkj 。根据最优原理,可得到如下的动态规划递归式:k a y(i,i+s)= 获得上述最小值的k. 以上求c 的递归式可用递归或迭代的方法来求解。c( 1,q) 为用最优法计算矩阵链的消耗,k a y( 1 ,q) 为最后一步的消耗。其余的乘积可由k a y值来确定。
1. 递归方法
与求解0 / 1背包及图像压缩问题一样,本递归方法也须避免重复计算c (i, j) 和k a y(i, j),否则算法的复杂性将会非常高。
例3-13 设q= 5和r =(1 0 , 5 , 1 , 1 0 , 2 , 1 0),式中待求的c 中有四个c的s= 0或1,因此用动态规划方法可立即求得它们的值: c( 1 , 1 ) =c( 5 , 5 ) = 0 ;c(1,2)=50; c( 4 , 5 ) = 2 0 0。现计算C( 2,5 ):c( 2 , 5 ) = m i n {c( 2 , 2 ) +c(3,5)+50, c( 2 , 3 ) +c(4,5)+500, c( 2 , 4 ) +c( 5 , 5 ) + 1 0 0 } (1 5 - 5)其中c( 2 , 2 ) =c( 5 , 5 ) = 0;c( 2 , 3 ) = 5 0;c(4,5)=200 。再用递归式计算c( 3 , 5 )及c( 2 , 4 ) :c( 3 , 5 ) = m i n {c( 3 , 3 ) +c(4,5)+100, c( 3 , 4 ) +c( 5 , 5 ) + 2 0 } = m i n { 0 + 2 0 0 + 1 0 0 , 2 0 + 0 + 2 0 } = 4 0c( 2 , 4 ) = m i n {c( 2 , 2 ) +c( 3 , 4 ) + 1 0 ,c( 2 , 3 ) +c( 4 , 4 ) + 1 0 0 } = m i n { 0 + 2 0 + 1 0 , 5 0 + 1 0 + 2 0 } = 3 0由以上计算还可得k a y( 3 , 5 ) = 4,k ay( 2 , 4 ) = 2。现在,计算c(2,5) 所需的所有中间值都已求得,将它们代入式(1 5 - 5)得:
c(2,5)=min{0+40+50, 50+200+500, 30+0+100}=90且k a y( 2 , 5 ) = 2
再用式(1 5 - 4)计算c( 1 , 5 ),在此之前必须算出c( 3 , 5 )、c(1,3) 和c( 1 , 4 )。同上述过程,亦可计算出它们的值分别为4 0、1 5 0和9 0,相应的k a y 值分别为4、2和2。代入式(1 5 - 4)得:
c(1,5)=min{0+90+500, 50+40+100, 150+200+1000, 90+0+200}=190且k a y( 1 , 5 ) = 2
此最优乘法算法的消耗为1 9 0,由k a y(1,5) 值可推出该算法的最后一步, k a y(1,5) 等于2,因此最后一步为M1 2×M3 5,而M12 和M35 都是用最优法计算而来。由k a y( 1 , 2 ) = 1知M12 等于M11×M2