1 动态规划的定义
动态规划通常用来解决最优化问题,在这类问题中,通过做出一组选择来达到最优解。在做出每个选择的同时,通常会生成与原问题形式相同的子问题。当多于一个选择子集都生成相同的子问题时,动态规划技术通常就会很有效,其关键技术就是对每个这样的子问题都保存其解,当其重复出现时即可避免重复求解。
动态规划与分治方法相似,都是通过组合子问题的解来求解原问题。分治方法将问题划分为互不相交的子问题,递归地求解子问题,再将它们地解结合起来,求出原问题的解。与之相反,动态规划应用于子问题重叠的情况,即不同的子问题具有公共的子子问题(子问题的求解是递归进行的,将其划分为更小的子子问题)。在这种情况下,分治算法会做很多不必要的工作,它会反复地求解那些公共子子问题。而动态规划算法对每个子子问题只求解一次,将其解保存在一个表格中,从而无需每次求解一个子子问题时都重新计算,避免了这种不必要的计算工作。
动态规划方法通常用来求解最优化问题。这类问题可以有很多可行解。每个解都有一个值,希望寻找具有最优解的解。我们称这样的解为问题的一个最优解。
通常按如下4个步骤来设计一个动态规划:
- 刻画一个最优解的结构特征
- 递归地定义最优解的值
- 计算最优解的值,通常采用自底向上的方法
- 利用计算出的信息构造一个最优解
步骤1~3时动态规划算法求解问题的基础。如果仅仅需要一个最优解的值,而非解本身,可以忽略步骤4。如果确实要做步骤4,有时就需要在执行步骤3的过程中维护一些额外信息,以便用来构造一个最优解。
2 问题求解
2.1 钢条切割
Serling公司购买长钢条,将其切割为短钢条出售。切割工序本身没有成本支出。公司管理层希望知道最佳的切割方案。假定Serling公司出售一段长度为i英寸的钢条的价格为\(p_i(i=1, 2, \dots ,\) 单位为美元)。钢条的长度均为整英寸。下图是对应价格表:
钢条切割问题是这样的:给定一段长度为n英寸的钢条和一个价格表 \(p_i(i=1, 2, \dots, n)\) , 求切割钢条方案,使得销售收益\(r_n\)最大。注意,如果长度为n英寸的钢条的价格\(p_n\)足够大,最优解可能就是完全不需要切割。
考虑 n = 4 的情况。下图给出了 4 英寸钢条所有可能的切割方案,包括根本不切割的方案。发现,将一段长度为4英寸的钢条切割为两段各长2 英寸的钢条,将产生\(p_1+p_2=5+5=10\)的收益,为最优解。
长度为n英寸的钢条共有 \(2^{n-1}\) 种不同的切割方案,因为在距离钢条左端 \(i(i=1, 2, \dots, n-1)\) 英 寸处,总是可以选择切割或不切割。用普通的加法符号表示切割方案,因此 7=2+2+3 表示将长度为 7 英寸的钢条切割为三段——两段长度为 2 英寸、一段长度为 3 英寸。 如果一个最优解将钢条切割为 k 段(对某个 \(1 \leq k \leq n \),那么最优切割方案
\(n = i_1 + i_2 + \dots + i_k\)
将钢条切割为长度分别为 \(i_1, i_2, \dots, i_k\)的小段,得到最大收益
\(r_n = p_{i_1} + p_{i_2} + \dots + p_{i_k}\)
对于上述价格表样例,我们可以观察所有最优收益值\(r_i(i=1, 2, \dots, 10)\)及对应的最优切割方案:
- \(r_1 = 1\),切割方案1 = 1(无切割)
- \(r_2 = 5\),切割方案2 = 2(无切割)
- \(r_3 = 8\),切割方案3 = 3(无切割)
- \(r_4 = 10\),切割方案4 = 2+2
- \(r_5 = 13\),切割方案5 = 2+3
- \(r_6 = 17\),切割方案6 = 6(无切割)
- \(r_7 = 18\),切割方案7 = 1+6 或 7 = 2+2+3
- \(r_8 = 22\),切割方案8 = 2+6
- \(r_9 = 25\),切割方案9 = 3+6
- \(r_10 = 30\),切割方案10 = 10(无切割)
更一般地,对于\(r_n(n \geq 1 )\),可以用更短的钢条的最优切割收益来描述它:
\(r_n = max(p_n, r_1+r_{n-1}, r_2+r_{n-2}, \dots, r_{n-1}+r_1)\)
第一个参数\(p_n\)对应不切割,直接出售长度为n英寸的钢条的方案。其他 n-1 个参数对应另外 n-1 种方案:对每个\(i = 1, 2, \dots, n-1\),首先将钢条切割为长度为 i 和 n-i 的两段,接着求解这两段的最优切割收益\(r_i\)和\(r_{n-i}\)(每种方案的最优收益为两段的最优收益之和)。由于无法预知哪种方案会获得最优收益,必须考察所有可能的 i ,选取其中收益最大者。如果直接出售原钢条会获得最大收益,当然可以选择不做任何切割。
注意到,为了求解规模为 n 的原问题,先求解形式完全一样,但规模更小的子问题。即当完成首次切割后,将两段钢条看成两个独立的钢条切割问题实例(这个过程类似于分成左右两边后,对左右两边分别进行递归,处理像归并排序一样)。通过组合两个相关子问题的最优解,并在所有可能的两段切割方案中选取组合收益最大者,构成原问题的最优解。称钢条切割问题满足最优子结构性质:问题的最优解由相关子问题的最优解组合而成,而这些子问题可以独立求解。
除了上述求解方法外,钢条切割问题还存在一种相似的但更为简单的递归求解方法:我们将钢条从左边切割下长度为 i 的一段,只对右边剩下的长度为 n-i 的一段继续进行切割(递归求解), 对左边的一段则不再进行切割。即问题分解的方式为:将长度为 n 的钢条分解为左边开始一段, 以及剩余部分继续分解的结果。这样,不做任何切割的方案就可以描述为:第一段的长度为n,收 益为\(p_n\), 剩余部分长度为 0 , 对应的收益为 \(r_0 = 0\)。于是可以得到收益公式的简化版本:
\(r_n = max_{1 \leq i \leq n}(p_i+r_{n-i})\)
在此公式中,原问题的最优解只包含一个相关子问题(右端剩余部分)的解,而不是两个。
自顶向下递归实现,以下代码实现了公式\(r_n = max_{1 \leq i \leq n}(p_i+r_{n-i})\)。
/*
p[]代表价格数组
n代表长度为n的钢条
*/
int CutRod(int p[],int n){
if(n==0) // 若n为0,则收益也对应0,所以返回0
return 0;
int q=-2e9;
for(int i=1;i<=n;i++){
q=max(q,p[i]+CutRod(p,n-i));
}
return q;
}
递归的效率非常差,因为CutRod函数反复地用相同的参数值对自身进行递归调用,即它反复求解相同的子问题,下图显示了n=4时的调用过程:
递归调用的时间复杂度为\(O(2^n)\)。
使用动态规划方法求解最优钢条切割问题
朴素递归算法之所以效率很低,是因为它反复求解相同的子问题。因此,动态规划方法仔细安排求解顺序,对每个子问题只求解一次,并将结果保存下来。如果随后再次需要此子问题的解,只需查找保存的结果,而不必重新计算。因此,动态规划方法是付出额外的内存空间来节省计算时间,是典型的时空权衡的例子。而时间上的节省可能是非常巨大的:可能将一个指数时间的解转化为一个多项式时间的解。如果子问题的数量是输入规模的多项式函数,而我们可以在多项式时间内求解出每个子问题,那么动态规划方法的总运行时间就是多项式阶的。
动态规划有两种等价的实现方法:
第一种方法称为带备忘的自顶向下法。此方法仍按自然的递归形式编写过程,但过程会保存每个子问题的解(通常保存在一个数组或散列表中)。当需要一个子问题的解时,过程首先检查是否已经保存过此解。如果是,则直接返回保存的值,从而节省了计算时间;否则,按通常方式计算这个子问题。我们称这个递归过程是带备忘的, 因为它”记住”了之前已经计算出的结果。
第二种方法称为自底向上法。这种方法一般需要恰当定义子问题“规模”的概念,使得任何子问题的求解都只依赖于“更小的”子问题的求解。因而我们可以将子问题按规模排序,按由小至大的顺序进行求解。当求解某个子问题时,它所依赖的那些更小的子问题都已求解完毕,结果已经保存。每个子问题只需求解一次,当我们求解它(也是第一次遇到它)时,它的所有前提子问题都已求解完成。
两种方法得到的算法具有相同的渐近运行时间。仅有的差异是在某些特殊情况下,自顶向下方法并未真正递归地考察所有可能的子问题。由于没有频繁的递归函数调用的开销,自底向上方法的时间复杂性函数通常具有更小的系数。
下面是自顶向下的代码,加入了备忘机制:
/*
p[]代表价格数组
n代表长度为n的钢条
*/
int Memoized_CutRod(int p[], int n){
int r[n+1]; // r代表不同尺寸的钢条最大收益的备忘数组
for(int i=0;i<=n;i++){
r[i]=-2e9; // 意为 初始化为无穷小
}
return memorized_CutRod_Aux(p, n, r);
}
int Memoized_CutRod_Aux(int p[], int n, int r[]){
if(r[n]>=0)return r[n]; // 如果备忘数组存在钢条长度为n的最大收益,那么直接返回。
int q;
if(n==0) // 若n为0,则收益也对应0,所以q=0
q = 0;
else{ // 否则的话,则计算最大收益(通过递归的方式)
q=-2e9;
for(int i=1;i<=n;i++){
q=max(q,p[i]+Memoized_CutRod_Aux(p,n-i,r));
}
}
r[n]=q; // 记录计算出来的钢条长度为n的最大收益,防止递归重复计算子问题,降低时间复杂度。
return q;
}
自底向上版本更为简单:
/*
p[]代表价格数组
n代表长度为n的钢条
*/
int Bottom_Up_CutRod(int p[], int n){
int r[n+1]; // r代表不同尺寸的钢条最大收益的备忘数组
r[0]=0; // 将长度为0的钢条的最大收益初始化为0
for(int i=1;i<=n;i++){ // 向上递推长度为i的钢条的最大收益,类似斐波那契
int q=-2e9;
for(int j=1;j<=i;j++){ // 长度为i的钢条,可以在j=1,2,...,i的位置进行切割,右边切割部分调用其最大收益r[i-j]
q=max(q,p[i]+r[i-j]);
}
r[i]=q; // 记录长度为i的钢条的最大收益
}
return r[n];
}
自底向上版本Bottom_Up_CutRod采用子问题的自然顺序,过程依次求解规模为i=0, 1, …, n的子问题。
自底向上算法和自顶向下算法的时间复杂度都为\(O(n^2)\)。
子问题图
下图显示了 n = 4 时钢条切割问题的子问题图。
它是一个有向图,每个顶点唯一地对应一个子问题。 若求子问题 x 的最优解时需要直接用到子问题 y 的最优解,那么在子问题图中就会有一条从子问题 x 的顶点到子问题 y 的顶点的有向边。例如,如果自顶向下过程在求解 x 时需要直接递归调用自身来求解 y 那么子问题图就包含从 x 到 y 的一条有向边。可以将子问题图看做自顶向下递归调用树的”简化版”或 “收缩版”,因为树中所有对应相同子问题的结点合并为图中的单一顶点,相关的所有边都从父结点指向子结点。
自底向上的动态规划方法处理子问题图中顶点的顺序为:对于一个给定的子问题 x ,在求解它之前求解邻接至它的子问题 y。自底向上动态规划算法是按”逆拓扑序”或”反序的拓扑序”来处理子问题图中的顶点。换句话说,对于任何子问题,直至它依赖的所有子问题均已求解完成,才会求解它。
子问题图 G=(V, E) 的规模可以帮助我们确定动态规划算法的运行时间。由于每个子问题只求解一次,因此算法运行时间等于每个子问题求解时间之和。通常,一个子问题的求解时间与 子问题图中对应顶点的度(出射边的数目)成正比,而子问题的数目等于子问题图的顶点数。因此,通常情况下,动态规划算法的运行时间与顶点和边的数量呈线性关系。
重构解
前文给出的钢条切割问题的动态规划算法返回最优解的收益值,但并未返回解本身(一个长度列表,给出切割后每段钢条的长度)。可以扩展动态规划算法,使之对每个子问题不仅保存最优收益值,还保存对应的切割方案。利用这些信息,就能输出最优解。
下面给出的是Bottom_Up_CutRod的扩展版本,它对长度为 j 的钢条不仅计算最大收
益值\(r_j\),还保存最优解对应的第一段钢条的切割长度\(s_j\):
int s[N]; // 记录不同尺寸的钢条的最优切割位置
int r[N]; // 记录不同尺寸的钢条的最优收益
void Extended_Bottom_Up_CutRod(int p[],int n){
r[0]=0;
for(int i=1;i<=n;i++){
int q=-2e9;
for(int j=1;j<=i;j++){
if(q<p[i]+r[i-j]){
q=p[i]+r[i-j];
s[i]=j; // 代表长度为i的钢条在j的位置进行了切割
}
}
r[i]=q;
}
return r;
}
// 打印最优切割方案
void Print_CutRod_Solution(int p[],int n){
Extended_Bottom_Up_CutRod(p,n);
while(n>0){
cout<<s[n]; // 长度为n的钢条切割在s[n]的位置
n=n-s[n]; // s[n]位置的右侧钢条又进行了切割,所以n=n-s[n],再次进行循环。
}
}
2.2 矩阵链乘法
这个例子是求解矩阵链相乘问题的动态规划算法。给定一个n个矩阵的序列(矩阵链)\(<A_1, A_2, \dots, A_n>\),希望计算它们的乘积
\(<A_1A_2 \dots A_n>\)
为了计算上式,可以先用括号明确计算次序,然后利用标准的矩阵相乘算法进行计算。由于矩阵乘法满足结合律,因此任何加括号的方法都会得到相同的计算结果。我们称有如下性质的矩阵乘积链为完全括号化的:它是单一矩阵,或者是两个完全括号化的矩阵乘积链的积,且已外加括号。例如,如果矩阵链为\(<A_1, A_2, A_3, A_4>\),则共有 5 种画完全括号化的矩阵乘积链:
\((A_1 ( A_2 ( A_3 A_4 )))\)
\((A_1 (( A_2 A_3 ) A_4 ))\)
\(((A_1 A_2 ) (A_3 A_4 ))\)
\(((A_1 ( A_2 A_3 )) A_4 )\)
\((((A_1 A_2 ) A_3 ) A_4 )\)
对矩阵链加括号的方式会对乘积运算的代价产生巨大影响。先来分析两个矩阵相乘的代价。下面的代码给出了两个矩阵相乘的标准算法。
int C[N][N]; // 存储矩阵相乘的结果
bool flag=true; // flag为true代表两矩阵可以合法相乘,否则不能
/*
A, B代表两矩阵相乘的数组
rowA, rowB分别代表A, B数组的行
columnA, columnB分别代表A, B数组的列
*/
int Matrix_Multiply(int A[][], int B[][], rowA, rowB, columnA, columnB){
if(columnA != rowB){ // 说明不满足两矩阵相乘的条件 即第一个数组的列要等于第二个数组的行才能相乘。
cout<<"incompatible dimensions";
}
else{
// 矩阵相乘的计算过程
for(int i=1;i<=rowA;i++){
for(int j=1;j<=columnB;j++){
C[i][j]=0;
for(int k=1;k<=columnA;k++){
C[i][j]+=A[i][k]*B[k][j];
}
}
}
}
}
两个矩阵 A 和 B 只有相容,即 A 的列数等于 B 的行数时,才能相乘。如果 A 是 \(p \times q\) 的矩阵,B是\(p \times q\)的矩阵,那么乘积 C 是\(p \times q\)的矩阵。计算 C 所需时间由第8 行的标量乘法的次数决定,即 pqr ,下文中将用标量乘法的次数来表示计算代价。
我们以矩阵链\(<A_1, A_2, A_3>\)相乘为例,来说明不同的加括号方式会导致不同的计算代价。 假设三个矩阵的规模分别为\( 10 \times 100、100 \times 5和 5 \times 50\)。如果按\(((A_1 A_2)A_3)\)的顺序计算,为计算\(A_1 A_2\)( 规模 \(10 \times 5\)),需要做\(10 \times 100 \times 5 =5 000\)次标量乘法,再与\(A_3\)相乘又需要做\(10ᅠ\times 5 \times 50 = 2 500\)次标量乘法,共需7 500次标量乘法。如果按\((A_1 (A_2A_3))\)的顺序,计算\(A_2A_3\)(规 模 \(100 \times 50\)),需\(100 \times 5 \times 50=25 000\)次标量乘法,\(A_1\) 再与之相乘又需\(10ᅠ\times 100ᅠ\times 50=50 000\)次标量乘法, 共需75 000次标量乘法。因此,按第一种顺序计算矩阵链乘积要比第二种顺序快10倍。
矩阵链乘法问题可描述如下:给 定 n 个矩阵的链\(〈A_1 ,A_2 ,\times ,A_n 〉\),矩阵\( A_i \)的规模为\(p_{i-1} \times p_i(1 \leq i \leq n)\),求完全括号化方案,使得计算乘积\(A_1A_2 \cdots A_n\)所需标量乘法次数最少。
注意,求解矩阵链乘法问题并不是要真正进行矩阵相乘运算,目标只是确定代价最低的计算顺序。确定最优计算顺序所花费的时间通常要比随后真正进行矩阵相乘所节省的时间 (例如仅进行7 500次标量乘法而不是75 000次)要少。
计算括号化方案的数量
在用动态规划方法求解矩阵链乘法问题之前,先来说服自己——穷举所有可能的括号化方案不会产生一个高效的算法。对一个n个矩阵的链,令 P(n) 表示可供选择的括号化方案的数量。当n = 1 时,由于只有一个矩阵,因此只有一种完全括号化方案。当\(n \geq 2 \) 时,完全括号化 的矩阵乘积可描述为两个完全括号化的部分积相乘的形式,而两个部分积的划分点在第k 个矩阵 和第k + 1个矩阵之间,k 为 1,2 , … ,n-1 中的任意一个值。因此,我们可以得到如下递归 公式:
\(P(n)= \begin{cases} 1 & 如果n=1 \\ \sum_{k=1}^{n-1} P(k)P(n-k) & 如果 n \geq 2 \end{cases}\)
递归公式的时间复杂度为\(O(2^n)\), 括号化方案的数量与n呈指数关系,通过暴力搜索穷尽所有可能的括号化方案来寻找最优方案,是一个糟糕的策略。
应用动态规划方法
下面用动态规划方法来求解矩阵链的最优括号化方案,还是按照本章开头提出的4个步骤进行:
- 刻画一个最优解的结构特征。
- 递归地定义最优解的值。
- 计算最优解的值,通常采用自底向上的方法。
- 利用计算出的信息构造一个最优解。
按顺序进行这几个步骤,清楚地展示针对本问题每个步骤应如何做。
步骤1:最优括号化方案的结构特征
动态规划方法的第一步是寻找最优子结构,然后就可以利用这种子结构从子问题的最优解构造出原问题的最优解。在矩阵链乘法问题中,此步骤的做法如下所述。为方便起见,用符号\(A_{i,j}(i \leq j)\)表示\(A_iA_{i+1} \cdots A_j\)乘积的结果矩阵。可以看出,如果问题是非平凡的,即 i < j,那么为了对\(A_iA_{i+1} \cdots A_j\)进行括号化,就必须在某个\(A_k\)和\(A_{k+1}\)之间将矩阵链划分开(k为\(i \leq k < j\)间的整数)。也就是说,对某个整数k,我们首先计算矩阵\(A_{i..k}\)和\(A_{k+1..j}\),然后再计算它们的乘积得到最终结果\(A_{i..j}\)。此方案的计算代价等于矩\(A_{i..k}\)的计算代价,加上矩阵\(A_{k+1..j}\)的计算代价, 再加上两者相乘的计算代价。
假设 \(A_iA_{i+1} \cdots A_j\) 的最优括号化方案的分割点在 \(A_k\) 和 \(A_{k+1}\)之间。那么,继续对”前缀”子链 \(A_iA_{i+1} \cdots A_k\) 进行括号化时,我们应该直接采用独立求解它时所得的最优方案。问题的最优解由相关子问题的最优解组合而成,而这些子问题可以独立求解,所以矩阵链乘法是满足最优子结构的。
一个非平凡的矩阵链乘法问题实例的任何解都需要划分链,而任何最优解都是由子问题实例的最优解构成的。因此,为了构造一个矩阵链乘法问题实例的最优解,我们可以将问题划分为两个子问题(\(A_iA_{i+1} \cdots A_k\)和 \(A_{k+1}A_{k+2} \cdots A_{j} \) 的最优括号化问题),求出子问题实例的最优解,然后将子问题的最优解组合起来。必须保证在确定分割点时,已经考察了所有可能的划分点, 这样就可以保证不会遗漏最优解。
步骤2: 一个递归求解方案
下面用子问题的最优解来递归地定义原问题最优解的代价。对矩阵链乘法问题,我们可以将对所有\(1 \leq i \leq j \leq n\)确定\(A_iA_{i+1} \cdots A_j\)的最小代价括号化方案作为子问题。令m[i, j]表示计算矩阵\(A_{i..j}\)所需标量乘法次数的最小值,那么,原问题的最优解—-计算\( A_{1..n} \)所需的最低代价就是m[1, n]。
我们可以递归定义m[i, j]如下。对于i=j时的平凡问题,矩阵链只包含唯一的矩阵\(A_{i..i} = A_{j}\),因此不需要做任何标量乘法运算。所以,对所有i = 1,2,… ,n,m[i, i] = 0。若\(i < j \),我 们利用步骤1中得到的最优子结构来计算m[i, j]。我们假设\(A_iA_{i+1} \cdots A_j\)的最优括号化方案的分割点在矩阵\(A_k\)和\(A_{k+1}\)之间,其中\(i \leq k < j\)。那么,m[i, j]就等于计算\(A_{i..k}\)和\(A_{k+1..j}\)的代价加上两者相乘的代价的最小值。由于矩阵\(A_i\)的大小为\(p_{i-1} \times p_i[latex],易知[latex]A_{i..k}\)与\(A_{k+1..j}\)相乘的代价为\(p_{i-1}p_kp_j\)次标量乘法运算。因此,得到
\(m[i, j]=m[i, k]+m[k+1,j]+p_{i-1}p_kp_j\)
此递归公式假定最优分割点 k 是已知的,但实际上我们是不知道的。不过,k只有 j-i 种可能的取值,即k = i, i+1, …, j-1。由于最优分割点必在其中,只需检查所有可能情况, 找到最优者即可。因此, \(A_iA_{i+1} /cdots A_j\)最小代价括号化方案的递归求解公式变为:
\(m[i,j] = \begin{cases} 0 & 如果 i = j \\ min_{i \leq k <j} \{ m[i,k] + m[k+1,j] + p_{i-1}p_kp_j \} & 如果 i < j \end{cases}\)
m[i,j] 的值给出了子问题最优解的代价,但它并未提供足够的信息来构造最优解。为此,
我们用 s[i,j] 保存\(A_iA_{i+1} /cdots A_j\)最优括号化方案的分割点位置k,即使得\(m[i, j]=m[i, k]+m[k+1,j]+p_{i-1}p_kp_j\)成立的k 值。
步骤3:计算最优代价
递归算法会在递归调用树的不同分支中多次遇到同一个子问题。这种子问题重叠的性质是应用动态规划的另一个标识(第一个标识是最优子结构)。
采用自底向上表格法代替基于上面的的递归求解公式来计算最优代价。下面给出的过程Matrix_Chain_Order实现了自底向上表格法。此过程假定矩阵\(A_i\)的规模为\(p_{i-1} \times p_i(i=1, 2, \cdots, n) \)。它的输入是一个序列\(p=<p_0, p_1, \cdots, p_n>\),其长度为p.length = n+1。过程用一个辅助表m[1..n, 1..n]来保存代价m[i, j], 用另一个辅助表s[1..n-1, 2..n]记录最优值m[i, j]对应的分割点k。我们就可以 利用表s构造最优解。
为了实现自底向上方法,我们必须确定计算m[i, j]时需要访问哪些其他表项显示,j-i+1 个矩阵链相乘的最优计算代价m[i, j]只依赖于那些少于j-i+1个矩阵链相乘的最优计算代价。也就是说,对k= i, i+1, …, j-1,矩阵\(A_{i..k}\)是k-i+1 < j-i+1个矩阵的积,矩阵\(A_{k+1..j}\)是j-k<j-i+1个矩阵的积。因此,算法应该按长度递增的顺序求解矩阵链括号化问题,并按对应的顺序填写表m。对矩阵链 \( A_i A_{i+1} \cdots A_j \) 最优括号化的子问题,我们认为其规模为链的长度j-i+1。
int m[N][N],s[N][N];
void Matrix_Chain_Order(int p[], int n){ // p是矩阵规模数组, n是p的尺寸
n=n-1;
for(int i=1;i<=n;i++){
m[i][i]=0;
}
for(int l=2;l<=n;l++){ // 第一重循环遍历矩阵链的长度
for(int i=1;i<=n-l+1;i++){ // 第二重循环遍历矩阵链的起点
j=i+l-1; // 通过矩阵链长度和矩阵链的起点计算矩阵链的终点
m[i][i]=2e9; // 初始化矩阵相乘代表为无穷大,方便后续比较替换。
for(int k=i;k<=j-1;k++){ // 第三重循环遍历矩阵链分割点
q=m[i][k]+m[k+1][j]+p[i-1]p[k]p[j];
if(q<m[i][j]){
m[i][j]=q;
s[i][j]=k; // 存储最优分割位置
}
}
}
}
}
算法首先对所有 i = 1 ,2 , … ,n 计算m[i, i] = 0 (长度为1 的链的最小计算代价)。接着在第一重for循环的第一个循环步中,利用递归公式对所有 i =1,2 , …, n-1 计算m[i, i+1](长度 l = 2 的链的最小计算代价)。在第二个循环步中,算法对所有i=1, 2 , …,n-2计算m[i, i+2](长度l= 3的链的最小计算代价),依此类推。在每个循环步中,计算代价m[i, j]时仅依赖于已经计算出的表项m[i, k]和m[k+1, j]。
下图展示了对一个长度为6 的矩阵链执行此算法的过程。由于定义m[i, j]仅在\(i \leq j\)时有意义,因此表m只使用主对角线之上的部分。图中的表是经过旋转的,主对角线已经旋转到了水平方向。矩阵链的规模列在了图的下方。在这种布局中,可以看到子矩阵链\( A_i A_{i+1} \cdots A_j \)相乘的代价m[i, j]恰好位于始于\(A_i\)的东北至西南方向的直线与始于\(A_j\) 的西北至东南方向的直线的交点上。表中同一行中的表项都对应长度相同的矩阵链。Matrix_Chain_Order按自下而上、自左至右的顺序计算所有行。当计算表项m[i, j]时,会用到乘积\(p_{i-1}p_kp_j(k=i, i+1, \cdots , j-1) \),以及m[i, j]西南方向和东南方向上的所有表项。
简单分析Matrix_Chain_Order的嵌套循环结构,可以看到算法的运行时间为\(O(n^3)\)。 循环嵌套的深度为三层,每层的循环变量(l、i和k)最多取n一1个值。算法还需要\(O(n^2)\)的内存空间来保存表m和s。 因此,Matrix_Chain_Order比起穷举所有可能的括号化方案来寻找最优解的指数阶算法要高效得多。
步骤4:构造最优解
虽然Matrix_Chain_Order求出了计算矩阵链乘积所需的最少标量乘法运算次数,但它并未直接指出如何进行这种最优代价的矩阵链乘法计算。表s[1..n-1, 2..n]记录了构造最优解所需的信息。每个表项s[i, j]记录了一个k值,指出\(A_iA_{i+1} \cdots A_j \)的最优括号化方案的分割点应在\(A_k\) 和\(A_{k+1}\) 之间。因此,我们知道\(A_{1..n}\)的最优计算方案中最后一次矩阵乘法运算应该是\(A_{1..s[1,n]}A_{s[1,n]+1..n}\)。我们可以用相同的方法递归地求出更早的矩阵乘法的具体计算过程,因为式s[1, s[1,n]]指出了计算A_{1..s[1,n]}时应进行的最后一次矩阵乘法运算;s[s[1,n]+1, n]指出了计算\(A_{s[1,n]+1..n}\)时应进行的最后一次矩阵乘法运算。下面给出的递归过程可以输出\(<A_i, A_{i+1}, \cdots, A_j>\)的最优括号化方案,其输入为Matrix_Chain_Order得到的表s及下标 i 和 j 。调用Print_Optimal_Parens(s, 1, n)即可输出\(<A_1, A_2, \cdots, A_n>\)的最优括号化方案。
void Print_Optimal_Parens(s, i, j){
if(i==j){
cout<<"A"<<i;
}
else{
cout<<"(";
Print_Optimal_Parens(s,i,s[i,j]);
Print_Optimal_Parens(s,s[i,j]+1,j);
cout<<")";
}
}