动态规划新版.docx

第6章 动态规划
动态规划是我们要介绍的几种算法设计方法中难度最大的一种，它建立在最优原则的基础上。采用动态规划方法，可以高效地解决许多用贪婪算法或分而治之算法无法解决的问题。
在介绍动态规划的原理之后，将分别考察动态规划方法在解决背包问题、图象压缩、矩阵乘法链、最短途径、无交叉子集等方面的应用。
1算法思想
和贪婪算法同样，在动态规划中，可将一个问题的解决方案视为一系列决策的结果。不同的是，在贪婪算法中，每采用一次贪婪准则便做出一个不可撤回的决策，而在动态规划中，还要考察每个最优决策序列中是否包含一个最优子序列。

例6-1[最短路经]考察下图中的有向图。假设要寻找一条从源节点s=1到目的节点d=5的最短途径，即选择此途径所通过的各个节点。第一步可选择节点2，3或4。假设选择了节点3，则此时所规定解的问题变成：选择一条从3到5的最短途径。假如3到5的途径不是最短的，则从1开始通过3和5的途径也不会是最短的。

所以在最短途径问题中，假如在的第一次决策时到达了某个节点v，那么不管v是如何拟定的，此后选择从v到d的途径时，都必须采用最优策略。

例6-2[0/1背包问题]考察前面的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=116。若设x1=1，则在本次决策之后，可用的背包容量为r=116-100=16。[x2，x3]=[0,1]符合容量限制的条件，所得值为15，但由于[x2，x3]=[1，0]同样符合容量条件且所得值为18，因此[x2，x3]=[0，1]并非最优策略。即x=[1，0，1]可改善为x=[1，1，0]。若设x1=0，则对于剩下的两种物品而言，容量限制条件为116。总之，假如子问题的结果[x2，x3]不是剩余情况下的一个最优解，则[x1，x2，x3]也不会是总体的最优解。

例6-3[航费]假设某航线价格表为：从亚特兰大到纽约或芝加哥，或从洛杉矶到亚特兰大的费用为$100；从芝加哥到纽约票价$20；而对于路经亚特兰大的旅客，从亚特兰大到芝加哥的费用仅为$20。从洛杉矶到纽约的航线涉及到对中转机场的选择。
洛杉矶
芝加哥
亚特兰大
纽约
100
100
20
100
20

假如问题状态的形式为（起点，终点），那么在选择从洛杉矶到亚特兰大后，问题的状态变为（亚特兰大，纽约）。从亚特兰大到纽约的最便宜航线是从亚特兰大直飞纽约，票价$100。而使用直飞方式时，从洛杉矶到纽约的花费为$200。但是，从洛杉矶到纽约的最便宜航线为洛杉矶-亚特兰大-芝加哥-纽约，其总花费为$140（在解决局部最优途径亚特兰大到纽约过程中选择了最低花费的途径：亚特兰大-芝加哥-纽约）。

假如用三维数组（tag，起点，终点）表达问题状态，其中tag为0表达转飞，tag为1表达其他情形，那么在到达亚特兰大后，状态的三维数组将变为（0，亚特兰大，纽约），它相应的最优途径是经由芝加哥的那条途径。

当最优决策序列中包含最优决策子序列时，可建立动态规划递归方程通常是根据最优准则，得出动态规划计算的“递归”计算式。

通常，这个递归式也是解决问题的关键，写出递归公式，然后用递归程序实现，并消除反复计算

在大多数情况下，求出递归方程是实际工作中的难点
而消除反复计算则成为问题是否可以实际求解的关键
（dynamic-programming recurrence equation），它可以帮助我们高效地解决问题。

例6-4[0/1背包]在例6-2的0/1背包问题中，最优决策序列由最优决策子序列组成。假设f(i,y)先把这个表达的地意义弄清楚?

表达例6-2中剩余容量为y，剩余物品为i，i+1，…，n时的最优解的值，即：
这个递归方程事实上将所有的情形所有考虑进去了，然后再从中找到一个最优解

运用最优序列由最优子序列构成的结论，可得到f的递归式。f(1,c)是初始时背包问题的最优解。可通过递归或迭代来求解f(1,c)。从f(n,*)开始迭式，先求出f(n,*)，然后递归计算f(i,*)(i=n-1，n-2，.，2)，最后得出f(1,c)。
对于例6-2（w=[100,14,10],p=[20,18,15]）可得：
f(3,y)=0 0≤y＜10
f(3,y)=15 y≥10
f(2,y)=0 0≤y＜10
f(2，y)=15 10≤y＜14
f(2，y)=18 14≤y＜24 和
f(2，y)=33 y≥24
剩2、3两个物品的情况

因此最优解：
f(1,116)= max{f（2，116），f（2，116-w1）+p1}
= max{f（2，116），f（2，16）+20}
= max{33，38} = 38。
现在计算xi值，环节如下：若f(1,c)=f(2,c)，则x1=0，否则x1=1。接下来需从剩余容量c-w1中寻求最优解，用f(2,c-w1)表达最优解。依此类推，可得到所有的xi(i=1…n)值。
在该例中，可得出f(2,116)=33≠f(1,116)，所以x1=1。接着运用返回值38-p1=18计算x2及x3，此时r=116-w1=16，又由f(2,16)=18，得f(3,16)=14≠f(2,16)，因此x2=1，此时r=16-w2=2，所以f(3,2)=0，即得x3=0。
动态规划方法采用最优原则（principle of optimality）来建立用于计算最优解的递归式。
所谓最优原则即不管前面的策略如何，此后的决策必须是基于当前状态（由上一次决策产生）的最优决策。由于对于有些问题的某些递归式来说并不一定能保证最优原则，因此在求解问题时有必要对它进行验证。若不能保持最优原则，则不可应用动态规划方法。在得到最优解的递归式之后，需要执行回溯（traceback）以构造最优解。

编写一个简朴的递归程序来求解动态规划递归方程是一件很诱人的事。然而，正如我们将在下文看到的，假如不努力地去避免反复计算，递归程序的复杂性将非常可观。假如在递归程序设计中解决了反复计算问题时，复杂性将急剧下降。动态规划递归方程也可用迭代方式来求解，这时很自然地避免了反复计算。事实上，这也是动态规划求解的两个基本方法：
1、 消除反复的递归
2、 采用迭代来代替递归

尽管迭代程序与避免反复计算的递归程序有相同的复杂性，但迭代程序不需要附加的递归栈空间，因此将比避免反复计算的递归程序更快。
2应用
2.1 0/1背包问题
背包问题的递归方程可以用几种方法具体的软件实现，即前面的递归公式求解，还是一个比较复杂的问题

求得。
1）递归策略
在例6-4中已建立了背包问题的动态规划递归方程，求解递归式（6-2）的一个很自然的方法便是使用程序6-1中的递归算法。该模块假设p、w和n为输入，且p为整型，F(1,c)返回f(1,c)值。
程序6-1背包问题的递归函数
实际程序中，w和p以及n需要“输入”到递归函数中

消除反复计算，假如y是整数，则可以考虑用一个二维数组保存F(i,y)，从而消除反复计算

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);
else
return max(F(i+1,y),F(i+1,y-w[i])+p[i]);
}

程序6-1的时间复杂性为：t(n)=Ο(2n)。
在所有f(i,y)计算完毕后，再用下列程序计算xi
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;
}

例6-5设n=5，p=[6,3,5,4,6]，w=[2,2,6,5,4]且c=10,求f(1,10)。为了拟定f(1,10)，调用函数F(1,10)。递归调用的关系如图6-1的树型结构所示。每个节点用y值来标记。对于第j层的节点有i=j，因此根节点表达F(1,10)，而它有左孩子和右孩子，分别相应F(2,10)和F(2,8)。总共执行了28次递归调用。
但我们注意到，其中也许具有反复前面工作的节点，如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)。假如保存以前的计算结果，则可将节点数减至19，由于可以丢弃图中的阴影节点。

图6-1 递归调用树

正如在例6-5中所看到的，程序6-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既可以用散列的形式存储，也可用二叉搜索树的形式存储。

2）权为整数的迭代方法
当权即 w
为整数时，可设计一个相称简朴的算法（称为Knapsack算法，见程序6-2）来求解f(1,c)。该算法基于例6-4所给出的策略，因此每个f(i,y)只计算一次。为简朴起见，程序6-2用二维数组f[][]来保存各f的值。
而回溯函数Traceback用于拟定由程序6-2所产生的xi值。
函数Knapsack的复杂性为(nc)，而Traceback的复杂性为(n)。
程序6-2 f和x的迭代计算下面程序中，要注意二维数组的真正实现，c语言中二维数组的传递需要特别注意

注意下列三种表达方式的差异：
T f[][N]
T **f
T *f[N]

template<class T>
T max( T x, T y )
{
return( (x>y)?x:y );
}

// 注意规定yMax用宏定义，大于c
template<class T>
void Knapsack( T p[],int w[],int c,int n,T (*f)[yMax] )
{//对于所有i和y计算f[i][y]
int i, y;
//初始化f[][]
for( y=0;y<yMax;y++)
f[n][y]=0;
for( y=w[n];y<=c;y++)
f[n][y]=p[n];
//计算剩下的f
for( i=n-1;i>1;i--)
{
for( y=0;y<yMax;y++)
f[i][y]=f[i+1][y];
for( 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)[yMax],int w[],int c,int n,int x[])
{//计算x
int i;
for( 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;
}

程序的具体实现如下（为简化起见，采用整型变量，同时，由于二维数组的传递比较麻烦，将数组f作为整形全局数组）：

#include <stdio.h>
const yMax=20; // 规定比c大
const n=5; // 物体数量
int f[n+1][yMax];

int max( int x, int y )
{
return( (x>y)?x:y );
}
void Knapsack(int p[],int w[],int c)
{//对于所有i和y计算f[i][y]
//初始化f[n][]
int i, y;
for( y=0;y<yMax;y++)
f[n][y]=0;
for( y=w[n];y<=c;y++)
f[n][y]=p[n];
//计算剩下的f
for(i=n-1;i>1;i--)
{
for( y=0;y<yMax;y++)
f[i][y]=f[i+1][y];
for( 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]);
}

void Traceback(int w[],int c,int x[])
{//计算x
int i;
for( 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;
}

void main()
{
int w[]={0,2,2,6,4,5}, p[]={0,6,3,5,4,6}, x[n+1];
int i,c=10;
Knapsack( p, w, c);
Traceback( w, c, x);
for( i=1; i<=n; i++ )
printf( "%d ", x[i] );
}

程序有两个缺陷：
l 规定权为整数；
l 当背包容量c很大时，程序6-2的速度慢于程序6-1。由于实际的计算中，f中的很多元素见前面实例计算 F（1，116）
按6-2 需要计算3*116个元素，而实际只用到其中大约7个元素

是多余计算的。
一般情况下，若c＞2n，程序6-2的复杂性为W (n2n)。

作为一个思考题，希望大家考虑，假如c特别大，为消除反复计算二引入的二维数组f就特别大，如何来解决存储问题，一般可以考虑不用数组，采用其他手段来存，如前面所述的三元组或者稀疏矩阵等。
2.2图像(字节流)压缩
下列方法实际也用于其他数据的压缩。
数字化图像是m×m的像素阵列。假定每个像素有一个0~255的灰度值。因此存储一个像素至多需8位。若每个像素存储都用最大位8位，则总的存储空间为8m2位。为了减少存储空间，我们将采用变长模式（variable bit scheme），即不同像素用不同位数来存储。
像素值为0和1时只需1位存储空间；值2、3各需2位；值4，5，6和7各需3位；以此类推，使用变长模式的环节如下：
l 图像线性化 根据图6-3a中的折线将m×m维图像转换为1×m2维矩阵。

图6-3 数字图像
a)“之”形的行主顺序 b)灰度值

l 分段 将像素组提成若干个段，分段原则是：每段中的像素位数相同。每个段是相邻像素的集合且每段最多含256个像素，因此，若相同位数的像素超过256个的话，则用两个以上的段表达。
l 创建文献 涉及三个部分：SegmentLength,BitsPerPixel和Pixels。第一个部分为所建的段的长度-1，其中各项均为8位长。第二部分：BitsPerPixel为各段中每个像素的存储位数-1，其中各项均为3位。第三部分Pixels则是以变长格式存储的像素的二进制串。
l 压缩文献 压缩所建立的文献，以减少空间需求。
上述压缩方法的效率（用所得压缩率表达）很大限度上取决于长段的出现频率。
例6-8考察图6-3b的4×4图像。按照蛇形的行主顺序，灰度值依次为10，9，12，40，50，35，15，12，8，10，9，15，11，130，160和240。各像素所需的位数分别为4，4，4，6，6，6，4，4，4，4，4，4，4，8，8，8，按等长的条件将像素分段，可以得到4个段[10，9，12]、[40，50，35]、[15,12,8,10,9,15,11]和［130,160,240］。
因此，第一部分SegmentLength为2，2，6，2；第二部分BitsPerSegment的内容为3，5，3，7；第三部分Pixels包含了按蛇形行主顺序排列的16个灰度值，其中头三个各用4位存储，接下来三个各用6位，再接下来的七个各用4位，最后三个各用8位存储。因此存储单元中前30位存储了前六个像素：
1010 1001 1100 111000 110010 100011
这三个部分需要的存储空间分别为：第一部分SegmentLength需32位；BitsPerSegment需12位；Pixels需82位，共需126位。而假如每个像素都用8位存储，则存储空间需8×16=128位，因而在本例图像中，节省了2位的空间。我们的问题不在这种方法能压缩多少
而是由此引入问题：
能否在此基础上，合并不同的段，从而减少所需的位数

进而，找出最优合并

假设在2)之后，产生了n个段。段标题（segment header）用于存储段的长度以及该段中每个像素所占用的位数。每个段标题需11位。现假设li和bi分别表达第i段的段长和该段每个像素的长度，则存储第i段像素所需要的空间为li*bi。
在2)中所得的三个部分的总存储空间为11n+ålibi。
可通过将某些相邻段合并的方式来减少空间消耗。如当段i和i+1被合并时，合并后的段长应为li+li+1。此时每个像素的存储位数为max{bi，bi+1}位。尽管这种技术增长了Pixels部分的空间消耗，但同时也减少了一个段标题的空间。

例6-9假如将例6-8中的第1段和第2段合并，合并后，SegmentLength部分变为5，6，2，BitsPerSegment部分变为5，3，7。而Pixels部分的前36位存储的是合并后的第一段：
001010 001001 001100 111000 110010 100011
其余的像素（例6-8第3段）没有改变。由于减少了1个段标题， SegmentLength部分和BitsPerPixel部分的空间消耗共减少了11位，而Pixels部分的空间增长6位，因此总共节约的空间为5位，空间总消耗为121位。

我们希望能设计一种算法，使得在产生n个段之后，能对相邻段进行合并，以便产生一个具有最小空间需求的新的段集合。
令sq为前q个段的最优合并所需要的空间。定义s0=0。考虑第i段(i＞0)，假如在最优合并C中，第i段与第i-1，i-2，.，i-r+1段相合并，而不涉及第i-r段。合并C所需要的空间消耗等于：
第1段到第i-r段所需空间+lsum(i-r+1,i)*bmax(i-r+1,i)+11
其中lsum(a,b)=，bmax(a,b)=max{ba,...,bb}。
假如在C中第1段到第i-r段的合并不是最优合并，那么需要对合并进行修改，以使其具有更小的空间需求。因此还必须对段1到段i-r进行最优合并，也即保证最优原则得以维持。故C的空间消耗为：
si=si-r+lsum（i-r+1,i）*bmax（i-r+1,i）+11
r的值介于1到i之间，其中规定lsum不超过256(由于段长限制在256之内)。尽管我们不知道如何选择r，但我们知道，由于C具有最小的空间需求，因此在所有选择中，r必须产生最小的空间需求。因此可得递归式：

假定kayi表达取得最小值时k的值，sn为n段的最优合并所需要的空间，因而一个最优合并可用kay的值构造出来。

例6-10假定得到五个段，它们的长度为[6，3，10，2，3]，总计24个像素（24BYTE），像素位数为[1，2，3，2，1]，要用公式计算sn，必须先求出sn-1，…，s0的值。s0为0，现计算s1：
s1=s0+l1*b1 +11=17
kay1=1
s2由下式得出：
s2=min{s1+l2b2,s0+(l1+l2)*max{b1,b2}}+11=min{17+6,0+9*2}+11=29
kay2=2
以此类推，可得s1…s5=[17，29，67，73，82]，kay1…kay5=[1，2，2，3，4]。

把s3写出来，可以看得更加清楚些：
s3=min{s2+l3b3 ,s1+(l2+l3)*max{b2,b3}, s0+(l1+l2+l3)*max{b1,b2,b3}}+11

由于s5=82，所以最优空间合并需82位的空间。可由kay5导出本合并的方式，由于kay5=4，表达其最小值是在：
s5 = s1 + lsum（2,5）*bmax（2,5）+11
时取得的，可以得出合并的方案，即1，2-5两段。

1）递归方法
用递归式可以递归地算出si和kayi。程序6-3为递归式的计算代码。l，b，和kay是一维的全局整型数组，L是段长限制（256），header为段标题所需的空间(11)。调用S(n)返回sn的值且同时得出kay值。调用Traceback(kay,n)可得到最优合并。
现给出程序6-3的复杂性：S的复杂性为t(n)=(2n)。Traceback的复杂性为(n)。

程序6-3递归计算s,kay及最优合并
int S(int i)
{//返回S(i)并计算kay[i]
if(i==0) return 0;
//k=1时,根据公式（6-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 n)
{//合并段
if(n==0)return;
Traceback(n-kay[n]);
cout<<"New segment begins at"<<(n-kay[n]+1)<<endl;
}

2）无反复计算的递归方法
通过避免反复计算si，可将函数S的复杂性减少到(n)。注意这里只有n个不同的si。
例6-11再考察例6-10中五个段的例子。当计算s5时，先通过递归调用来计算s4，…，s0。计算s4时，通过递归调用计算s3，…，s0，因此s4只计算了一次，而s3计算了两次，每一次计算s3要计算一次s2，因此s2共计算了四次，而s1反复计算了16次！
可运用一个数组s来保存先前计算过的si以避免反复计算。改善后的代码见程序6-4，其中s为初值为0的全局整型数组。

程序6-4避免反复计算的递归算法
int S(int i)
{//计算S(i)和kay[i]
//避免反复计算
if(i==0) return 0;
if(s[i]>0) return s[i];//已计算完
//计算s[i]
//一方面根据公式（6-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];
}

可以证明程序6-4的时间复杂性为(n)。

3）迭代方法
倘若用式（6-3）依序计算s1,…,sn，便可得到一个复杂性为(n)的迭代方法。在该方法中，在si计算之前，sj必须已计算好。该方法的代码见程序6-5，其中仍运用函数Traceback（见程序6-3）来获得最优合并。
程序6-5迭代计算s和kay
void Vbits(int n)
{//计算s[i]和kay[i]
int *s = new int[n+1];
int L=256,header=11;
s[0]=0;
//根据式（6-3）计算s[i]
for(int i=1;i<=n;i++)
{
//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];
if(s[i]>s[i-k]+lsum*bmax)
{
s[i]=s[i-k]+lsum*bmax;
kay[i]=k;
}
}
s[i]+=header;
}
}

2.3矩阵乘法链
m×n矩阵A与n×p矩阵B相乘需花费(mnp)的时间。我们把mnp作为两个矩阵相乘所需时间的测量值。现假定要计算三个矩阵A、B和C的乘积，有两种方式计算此乘积。在第一种方式中，先用A乘以B得到矩阵D，然后D乘以C得到最终结果，这种乘法的顺序可写为(A*B)*C。第二种方式写为A*(B*C),道理同上。尽管这两种不同的计算顺序所得的结果相同，但时间消耗会有很大的差距。

例6-12 一个较为极端的例子，假定A为100×1矩阵，B为1×100矩阵，C为100×1矩阵，则A*B的时间花费为10000，得到的结果D为100×100矩阵，再与C相乘所需的时间花费为1000000，因此计算(A*B)*C的总时间为1010000。B*C的时间花费为10000，得到的中间矩阵为1×1矩阵，再与A相乘的时间消耗为100，因而计算A*（B*C）的时间花费竟只有10100！并且，计算（A*B）*C时，还需10000个单元来存储A*B，而A*（B*C）计算过程中，只需用1个单元来存储B*C。

下面举一个得益于选择合适顺序计算A*B*C矩阵的实例：考虑两个3维图像的匹配。图像匹配问题的规定是，拟定一个图像需旋转、平移和缩放多少次才干逼近另一个图像。实现匹配的方法之一便是执行约100次迭代计算，每次迭代需计算一个12×1向量T：
T=åA(x,y,z)*B(x,y,z)*C(x,y,z)
其中A，B和C分别为12×3，3×3和3×1矩阵。(x,y,z)为矩阵中向量的坐标。设t表达计算A(x,y,z)*B(x,y,z)*C(x,y,z)的计算量。假定此图像含256×256×256个向量，在此条件中，这100个迭代所需的总计算量近似为100*2563*t≈1.7*109t。
若三个矩阵是按由左向右的顺序相乘的，则t=12*3*3+12*3*1=144；但假如从右向左相乘，t=3*3*1+12*3*1=45。由左至右计算约需2.4*1011个操作，而由右至左计算大约只需7.5*1010个操作。假如使用一个每秒可执行1000亿次操作的计算机，由左至右需2.4秒钟，而由右至左只需0.75秒钟。
在计算矩阵运算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)。用Mij表达链Mi×…×Mj（i≤j）的乘积。设c(i,j)为用最优法计算Mij的消耗，kay(i,j)为用最优法计算Mij的最后一步Mik×Mkj的k值。因此Mij的最优算法涉及如何用最优算法计算Mik和Mkj以及计算Mik×Mkj。根据最优原理，可得到如下的动态规划递归式我们用矩阵的规模近似地当作运算的复杂性

：

以上求c的递归式可用递归或迭代的方法来求解。c(1，q)为用最优法计算矩阵链的消耗，kay(1,q)为最后一步的消耗。其余的乘积可由kay值来拟定。

1）递归方法
与求解0/1背包及图像压缩问题同样，本递归方法也须避免反复计算c(i,j)和kay(i,j)，否则算法的复杂性将会非常高。
例6-13设q=5和r=（10,5,1,10,2,10），即：
[10*5][5*1][1*10][10*2][2*10]五个矩阵相乘。
由动态规划的递归式得：

式中待求的c中有四个c的s=0或1，因此用动态规划方法可立即求得它们的值：c(1,1)=c(5,5)=0;c(1,2)=50;c(4,5)=200。
现计算C(2，5)：
c(2,5)=min{ c(2,2)+c(3,5)+50,c(2,3)+c(4,5)+500,
c(2,4)+c(5,5)+100 }
其中c(2,2)=c(5,5)=0；c(2,3)=50；c(4,5)=200。再用递归式计算c(3,5)及c(2,4):
c(3,5)=min{ c(3,3)+c(4,5)+100,c(3,4)+c(5,5)+20 }
=min{ 0+200+100,20+0+20 } = 40
c(2,4)=min{ c(2,2)+c(3,4)+10,c(2,3)+c(4,4)+100 }
=min{ 0+20+10,50+10+20 }= 30
由以上计算还可得kay(3,5)=4，kay(2,4)=2。现在,计算c(2,5)所需的所有中间值都已求得，将它们代入式（6-5）得：
c(2,5)=min{ 0+40+50,50+200+500,30+0+100 }=90
且kay(2,5)=2。再用式（6-4）计算c(1,5)，在此之前必须算出c(3,5)、c(1,3)和c(1,4)。同上述过程，亦可计算出它们的值分别为40、150和90，相应的kay值分别为4、2和2。代入式（6-4）得：
c(1,5)=min{ 0+90+500,50+40+100,150+200+1000,90+0+200 }
=190且kay(1,5)=2
此最优乘法算法的消耗为190，由kay(1,5)值可推出该算法的最后一步，kay(1,5)等于2，因此最后一步为M12×M35，而M12和M35都是用最优法计算而来。由kay(1,2)=1知M12等于M11×M22，同理由kay(3,5)=4得知M35由M34×M55算出。依此类推，M34由M33×M44得出。因而此最优乘法算法的环节为：
M11×M22=M12
M33×M44=M34
M34×M55=M35
M12×M35=M15
计算c(i,j)和kay(i,j)的递归代码见程序6-6。在函数C中，r为全局一维数组变量，kay是全局二维数组变量，函数C返回c(ij)之值且置kay[a][b]=kay(a,b)(对于任何a,b)，其中c(a,b)在计算c(i,j)时皆已算出。函数Traceback运用函数C中已算出的kay值来推导出最优乘法算法的环节。

设t(q)为函数C的复杂性，其中q=j-i+1（即Mij是q个矩阵运算的结果）。当q为1或2时，t(q)=d，其中d为一常数；而q>2时，，其中e是一个常量。因此当q＞2时，t(q)＞2t(q-1)+e，所以t(q)=W (2q)。函数Traceback的复杂性为(q)。
程序6-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[i+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<j;k++)
{
int t=C(i,k)+C(k+1,j)+r[i]*r[k+1]*r[j+1];
if(t<u)
{//小于最小值的情形
u=t;
kay[i][j]=k;
}
}
return u;
}
void Traceback(int i,int j,int **kay)
{//输出计算Mij的最优方法
if(i==j)return;
Traceback(i,kay[i][j],kay);
Traceback(kay[i][j]+1,j,kay);
cout<<"MultiplyM"<<i<<","<<kay[i][j];
cout<<"andM"<<(kay[i][j]+1)<<","<<j<<end1;
}

2）无反复计算的递归方法
若避免再次计算前面已经计算过的c（及相应的kay），可将复杂性减少到（q3）。而为了避免反复计算，需用一个全局数组c[][]存储c(i,j)值，该数组初始值为0。函数C的新代码见程序6-7：
程序6-7无反复计算的c(i,j)计算方法
int C(int i,int j)
{//返回c(i,j)并计算kay(i,j)=kay[I][j]
//避免反复计算
//检查是否已计算过
i