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 2，同理由k a y( 3 , 5) = 4得知M35 由M3 4×M55 算出。依此类推，M34 由M3 3×M44 得出。因而此最优乘法算法的步骤为：
M11×M2 2 = M1 2
M3 3×M4 4 = M3 4
M3 4×M5 5 = M3 5
M1 2×M3 5 = M1 5
计算c(i, j) 和k a y (i, j) 的递归代码见程序1 5 - 6。在函数C中，r 为全局一维数组变量， k a y是全局二维数组变量，函数C返回c(i j) 之值且置k a y [a] [b] =k ay (a , b) (对于任何a , b)，其中c(a , b)在计算c(i,j) 时皆已算出。函数Traceback 利用函数C中已算出的k a y值来推导出最优乘法算法的步骤。
设t (q)为函数C的复杂性，其中q=j-i+ 1（即Mij 是q个矩阵运算的结果）。当q为1或2时，t(q) =d，其中d 为一常数；而q> 2时，t (q)=2q-1?k = 1t (k ) +e q，其中e 是一个常量。因此当q＞2时，t(q)＞2t (q- 1 ) +e，所以t (q)= W ( 2q)。函数Traceback 的复杂性为(q)。
程序15-6 递归计算c (i, j) 和kay (i, j)
int C(int i, int j)
{ / /返回c(i,j) 且计算k(i,j) = kay[i][j]
if (i==j) return 0; //一个矩阵的情形
if (i == j-1) { //两个矩阵的情形
kay[i][i+1] = i;
return r[i]*r[i+1]*r[r+2];}
/ /多于两个矩阵的情形
/ /设u为k = i 时的最小值
int u = C(i,i) + C(i+1,j) + r[i]*r[i+1]*r[j+1];
kay[i][j] = i;
/ /计算其余的最小值并更新u
for (int k = i+1; k < t =" C(i,k)" u =" t;" i ="="">） return c[i][j];
/ /若未计算,则进行计算
if(i==j) return 0; //一个矩阵的情形
i f ( i = = j - 1 ) { / /两个矩阵的情形
kay[i][i+1]=i;
c [ i ] [ j ] = r [ i ] * r [ i + 1 ] * r [ i + 2 ] ;
return c[i][j];}
/ /多于两个矩阵的情形
/ /设u为k = i 时的最小值
int u=C(i,i)+C(i+1,j)+r[i]*r[i+1]*r[j+1];
k a y [ i ] [ j ] = i ;
/ /计算其余的最小值并更新u
for (int k==i+1; k
int t=C(i,k)+C(k+1,j)+r[i]*r[k+1]*r[j+1];
if (t
u = t ;
k a y [ i ] [ j ] = k ; }
}
c [ i ] [ j ] = u ;
return u;
}
为分析改进后函数C 的复杂性，再次使用分期计算方法。注意到调用C(1, q) 时每个c (i, j)（1≤i≤j≤q）仅被计算一次。要计算尚未计算过的c(a,b)，需附加的工作量s =j-i＞1。将s 计入第一次计算c (a, b) 时的工作量中。在依次计算c(a, b) 时，这个s 会转计到每个c (a, b) 的第一次计算时间c 中，因此每个c (i, i) 均被计入s。对于每个s，有q-s+ 1个c(i, j) 需要计算，因此总的工作消耗为q-1 ?s=1(q-s+ 1) = (q3 )。
3. 迭代方法
c 的动态规划递归式可用迭代的方法来求解。若按s = 2，3，.，q-1 的顺序计算c (i, i+s)，每个c 和kay 仅需计算一次。
例3-14 考察例3 - 1 3中五个矩阵的情况。先初始化c (i, i) (0≤i≤5) 为0，然后对于i=1, ., 4分别计算c (i, i+ 1 )。c (1, 2)= r1 r2 r3 = 5 0，c (2, 3)= 5 0，c ( 3,4)=20 和c (4, 5) = 2 0 0。相应的k ay 值分别为1，2，3和4。
当s= 2时，可得：
c( 1 , 3 ) = m i n {c( 1 , 1 ) +c(2,3)+ r1 r2 r4 , c( 1 , 2 ) +c( 3 ,3 )+r1r3r4 }=min{0+50+500,50+0+100}
=150
且k a y( 1 , 3 ) = 2。用相同方法可求得c( 2 , 4 )和c( 3 , 5 )分别为3 0和4 0，相应k a y值分别为2和3。
当s= 3时，需计算c(1,4) 和c( 2 , 5 )。计算c(2,5) 所需要的所有中间值均已知(见( 1 5 - 5 )式)，代入计算公式后可得c( 2 , 5 ) = 9 0，k a y( 2 , 5 ) = 2。c( 1 , 4 )可用同样的公式计算。最后，当s= 4时，可直接用（1 5 - 4）式来计算c( 1 , 5 )，因为该式右边所有项都已知。
计算c 和kay 的迭代程序见函数M a t r i x C h a i n（见程序1 5 - 8），该函数的复杂性为(q3 )。计算出kay 后同样可用程序1 5 - 6中的Traceback 函数推算出相应的最优乘法计算过程。
程序15-8 c 和kay 的迭代计算
void MatrixChain(int r[], int q, int **c, int **kay)
{// 为所有的Mij 计算耗费和k a y
// 初始化c[i][i], c[i][i+1]和k a y [ i ] [ i + 1 ]
for (int i = 1; i < q; i++) {
c[i][i] = 0;
c[i][i+1] = r[i]*r[i+1]*r[i+2];
kay[i][i+1] = i;
}
c[q][q] = 0;
/ /计算余下的c和k a y
for (int s = 2; s < q; s++)
for (int i = 1; i <= q - s; i++) {
// k = i时的最小项
c[i][i+s] = c[i][i] + c[i+1][i+s] + r[i]*r[i+1]*r[i+s+1];
kay[i][i+s] = i;
// 余下的最小项
for (int k = i+1; k < i + s; k++) {
int t = c[i][k] + c[k+1][i+s] + r[i]*r[k+1]*r[i+s+1];
if (t < c[i][i+s]) {// 更小的最小项
c[i][i+s] = t;
kay[i][i+s] = k;}
}
}
}
没有评论:
发表评论