矩阵连乘——动态规划算法
目录
预备知识
了解矩阵连乘相关算法的前提是了解矩阵相乘的规则,此处应该询问度娘。
什么是矩阵连乘的括号化方案?举个例子来理解。
设有四个矩阵A、B、C、D,它们的维数分别是:A=50*10、B=10*40、C=40*30、D=30*5,总共有五种完全加括号的方式,如下:(A((BC)D))、(A(B(CD)))、((AB)(CD))、(((AB)C)D)、((A(BC))D). 如果看到这里还不懂矩阵相乘的相关知识,那么你可以放弃了,这个算法不适合你。
前言
对于矩阵连乘问题可以有多种解决办法,例如穷举法、递归与分治法、动态规划等。顾名思义,穷举法即把所有情况一一列举出来然后进行比较,当矩阵数目越来越大时,显然该方法不现实。采用递归与分治的思想比穷举法的可行性更强,但是递归的同时会重复计算同一个子问题,因此该算法依然不够理想。该篇主要讲解用动态规划算法如何解决矩阵连乘问题。
问题
给定n个矩阵{A1,A2,…,An },其中Ai 与Ai+1 是可乘的,i=1,2,3,···,n-1. 如何确定计算矩阵连乘积的计算次序,使得依次次序计算矩阵连乘积需要的数乘次数最少?
应用动态规划方法的步骤
-
刻划一个最优解的结构特征;
-
递归地定义最优解的值;
-
计算最优解的值,通常采用自底向上的方法;
-
利用计算出的信息构造一个最优解。
按照步骤分析问题
步骤1:最优括号化方案的结构特征
我们用符号Ai·j (i<=j)表示AiAi+1···Aj 乘积的结果矩阵,如果问题是非凡的,即i<j,那么为了对AiAi+1···Aj 进行括号化,我们就必须在某个Ak 和Ak+1之间将矩阵链划分开(k为i<=k<j之间的整数),也就是说我们首先计算矩阵Ai···k 和Ak+1···j ,然后再计算他们的乘积得到最终结果Ai···j 。现在问题就变成了对矩阵Ai···k 和Ak+1···j 进行独立求解,并且要满足最优括号化方案,最后将子问题的最优解合并就可以得到问题的最优解。
步骤2:一个递归求解方案
令m[i,j]表示计算矩阵Ai···j 数乘次数的最小值,那么原问题的最优解——计算A1···n 所需的最低代价就是m[1,n]。i=j时不需要做任何标量乘法运算,所以m[i,i]=0.若i<j,利用步骤1中得到的最优子结构来计算m[i,j]。设矩阵Ai···j 的大小为pi-1×pi ,那么Ai···k 和Ak+1···j 相乘的代价为pi-1*pk*pj 。因此可得.
此递归公式假定最优分割点k是已知的,但实际上我们并不知道。不过,k只有j-i种可能的取值,即k=i,i+1,···,j-1.由于最优分割点必其中,我们只需要检查所有可能的情况,推到最优者即可。因此,AiAi+1···Aj 最小代价括号化方案的递归求解公式变为
步骤3:计算最优代价
注意到,我们需要求解的不同子问题的数目是相对较少的,这种子问题重叠的性质是应用动态规划的另一个标识。假定矩阵Ai 的规模为pi-1×pi 。它的输入是一个序列p=<p0,p1,···,pn >,其长度为p.length=n+1.过程用一个辅助表来保存m[i,j],用另一个辅助表s[1···n-1, 2···n]记录最优值相应的分割点k,这样我们就可以利用表s构造最优解。
步骤4:构造最优解
虽然步骤3中求出了计算矩阵链乘积所需的最少标量乘法运算次数,但并未直接指出如何进行这种最优代价的矩阵链乘法计算。表s记录了最优解所需的信息,即记录了所有的k值,通过查询s表即可得到最终的答案。
算法的核心代码
void MatrixChain(int *p, int n, int**m ,int**s)
{
for(int i=1; i<= n; i++) m[i][i]=0; //最小的子问题,矩阵个数为1
for(int r= 2; r<= n; r++) //r代表子问题的矩阵个数,从2到n
{
for(int i=1;i<=n-r+1; i++){ //分析n-r+1种情况
int j=i+r-1;
m[i][j]=m[i+1][j]+p[i-1]*p[i]*p[j];//假设每次都从第i个矩阵后分开
s[i][j]=i;
for(int k=i+1;k<j;k++){ //讨论k的所有可能情况
int t= m[i][k]+m[k+1][j]+ p[i-1]*p[k]*p[j];
if(t<m[i][j]) {m[i][j]=t; s[i][j]=k;}}}//选择数乘次数最小的情况,并将相关数据覆盖于表中的相应位置
}
}
举例验证算法
A1 | A2 | A3 | A4 | A5 | A6 |
30*35 | 35*15 | 15*5 | 5*10 | 10*20 | 20*25 |
例如:
算法复杂度分析
算法的主要计算量取决于算法中对r、i和k的三重循环。循环体内的计算量为O(1),而三重循环的总次数为O(n3 ),因此算法的计算时间上界为O(n3 ),算法所占用的空间显然为O(n2 ).
算法的基本要素
1、最优子结构
矩阵连乘计算次序问题的最优解包含着其子问题的最优解。这种性质称为最优子结构性质。
利用问题的最优子结构性质,以自底向上的方式递归地从子问题的最优解逐步构造出整个问题的最优解。最优子结构是问题能用动态规划算法求解的前提。
同一个问题可以有多种方式刻划它的最优子结构,有些表示方法的求解速度更快(空间占用小,问题的维度低)。
2、重叠子问题
递归算法求解问题时,每次产生的子问题并不总是新问题,有些子问题被反复计算多次。这种性质称为子问题的重叠性质。
动态规划算法,对每一个子问题只解一次,而后将其解保存在一个表格中,当再次需要解此子问题时,只是简单地用常数时间查看一下结果。
3、备忘录方法
矩阵连乘的备忘录方法
备忘录方法为每个已经计算的子问题建立备忘录, 即保存子问题的计算结果以备需要时引用,从而避免了相同子问题的重复求解。
算法MemoizedMatrixChain是解矩阵连乘最优计算次序问题的备忘录方法。
#include <iostream>
using namespace std;
int MemoizedMatrixChain(int n, int** &, int ** &,int *);
int LookupChain(int ,int ,int *,int** &, int ** &);
int MemoizedMatrixChain(int n, int** &m, int ** &s,int *p)
{
for ( int i=0; i <= n-1; i++){ //初始化表格
for ( int j = 0; j <= n-1; j++)
{
m[i][j] =0;
s[i][j] =0;
}
}
return LookupChain(1,n,p,m,s);
};
int LookupChain(int i,int j,int *p,int** &m, int ** &s) // 动态规划算法的核心代码,与前面提到的MatrixChain相似
{
if (i == j) return 0;
int u = LookupChain(i,i,p,m,s) + LookupChain(i+1,j,p,m,s) + p[i-1]*p[i]*p[j];
s[i][j] = i;
for (int k = i+1; k < j; k++) {
int t = LookupChain(i,k,p,m,s) + LookupChain(k+1,j,p,m,s) + p[i-1]*p[k]*p[j];
if (t < u) { u = t; s[i][j] = k;}
}
m[i][j] = u;
return u;
};
int main()
{
int n=6;
int p[]={30,35,15,5,10,20,25};
int** m;
int** s;
m=new int*[n]; //声明m表
for (int zz=0;zz<n;zz++)
m[zz]=new int[n];
s=new int*[n]; //声明s表
for (int ii=0; ii<n; ii++)
s[ii]=new int[n];
cout<<MemoizedMatrixChain(n,m,s,p)<<endl; //输出最小数乘次数
return 0;
}
运行结果:
算法复杂度分析:与动态规划算法一样,备忘录算法耗时O(n3 )。m[i][j]这些记录项的初始化耗费O(n2 )时间,每个记录项只填入一次。每次填入时,不包括填入其他记录项的时间,共耗费O(n)时间。因此通过备忘录技术,直接递归算法的计算时间从O(2n )降低到O(n3 )。
总之,为求解矩阵链乘法问题,我们既可以用带备忘录的自顶向下动态规划算法,也可以用自底向上的动态规划算法,时间复杂度均为O(n3 )。两种方法都利用了重叠子问题性质,通常情况下,自底向上动态规划算法会比自顶向下备忘算法快(都是O(n3 )时间,相差一个常量系数),因为自底向上算法没有递归调用的开销,表的维护开销也更小。
闲聊时刻
前几天看毕淑敏的书,在此跟大家分享一段话。
人不是不可以怯懦和懒惰,但他不能把这些陋习伪装成高风亮节,不能由于自己做不到高尚,就诋毁所有做到了这些的人是伪善。你可以跪在泥里,但你不可以把污泥抹在整个世界的胸膛,并因此煞有介事地说到处都是污垢。