上一篇介绍了动态规划的基本概念,又介绍了一个典型的例子链条分割,这一篇来介绍另一个例子-----矩阵链乘法
给定一个n个矩阵的序列,我们希望计算它的乘积 A1A2A3...An,我们可以用括号标定计算次序,然后再利用标准的矩阵乘法算法来计算,由于矩阵乘法满足结合律,所以任何一种计算顺序最后的计算结果都是相同的。
完全括号化的矩阵乘积链:它是单一矩阵,或者两个完全括号化的矩阵乘积链的积,并且已经外加括号;例如矩阵链(A1,A2,A3,A4),则有下边5种完全括号化的矩阵链:
对矩阵链加括号虽然对计算结果没影响,但是计算过程中所花费的代价却是影响巨大。我们先来做简单的分析,首先A,B矩阵能相乘的条件是两个矩阵相容(即A矩阵的列数等于B矩阵的行数),设A是n*m阶,B是m*r阶,那么要进行n*m*r次标量乘法运算;举个例子:设A1,A2,A3三个矩阵分别为10x100,100x5,5x50阶,若按(A1(A2A3))顺序计算,那么总共进行(100x5x50)+(10x100x50)=75000次标量乘法运算;若按((A1A2)A3)顺序,那么总共进行(10x100x5)+(10x5x50)=7500次标量乘法运算,第二种顺序比第一种运算快10倍。
引入矩阵乘法链问题:给定n个矩阵的链(A1,A2,...An),矩阵Ai的规模pi-1xpi,求完全括号化方案,使得A1A2...An的乘积所需的标量乘法次数最小。
在此之前先来考虑一下完全括号化方案的总数,用P(n)表示方案总数,当n=1时,只有一个矩阵,只有一种方案;n>=2时,这样考虑,完全括号化的矩阵乘积可以看成两个完全括号化的部分积的形式,总结成下边的式子:
下边根据动态规划的4个步骤来设计算法:
1、最优括号化方案的结构特征:
用Ai...j表示AiAi+1...Aj的乘积,为了进行完全括号化,必须在Ak,Ak+1之间将链断开,即对一个整数k,我们先计算Ai...k,Ai+1...j,再计算两者的乘积,此方案的代价是Ai...k,Ak+1...j的计算代价,加上两者相乘的计算代价。有一个重要的结论:在原问题Ai...Aj的最优括号化方案中,对子链Ai...Ak,Ak+1...Aj括号化的方案就是原问题的括号化方案。
2、一个递归求解方案:
用m[i,j]表示计算Ai...Aj所需的标量乘法次数的最小值,那么m[1,n]就表示整个问题的解。这样来定义递归式,当i=j时,输入只有一个矩阵,即m[i,i]=0,当i<j时,可以在Ak和Ak+1之间划分,那么m[i,j]就是Ai...Ak的次数和Ak+1...Aj的次数,再加上两者相乘次数的最小值。
最后再加一个量s[i,j]保存划分位置。
3、计算代价,给出实现:
自底向上表格法
//自底向上表格法
void Min_ChainOrder_down2up(DataType p[],int n,DataType m[][N],DataType s[][N])
{
for(int l=2;l<=n;l++)
{
for(int i=1;i<=n-l+1;i++)
{
int j=l+i-1;
m[i][j]=MAX_LIMITS;
for(int k=i;k<=j-1;k++)
{
int 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;
}
}
}
}
}
算法时间复杂度O(n^3),并且需要O(n^2)的内存来维护m和s这两个表
4、构造最优解:void print_result(DataType s[][N],int i,int j)
{
if(i==j)
cout<<"A"<<i;
else
{
cout<<"(";
print_result(s,i,s[i][j]);
print_result(s,s[i][j]+1,j);
cout<<")";
}
}
下边给出规模为6的测试代码,关于此算法的详细分析和理解另开一帖详细解释
#include <iostream>
#define N 7
#define MAX_LIMITS 100000
typedef int DataType;
using namespace std;
//自底向上表格法
void Min_ChainOrder_down2up(DataType p[],int n,DataType m[][N],DataType s[][N])
{
for(int l=2;l<=n;l++)
{
for(int i=1;i<=n-l+1;i++)
{
int j=l+i-1;
m[i][j]=MAX_LIMITS;
for(int k=i;k<=j-1;k++)
{
int 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;
}
}
}
}
}
void print_result(DataType s[][N],int i,int j)
{
if(i==j)
cout<<"A"<<i;
else
{
cout<<"(";
print_result(s,i,s[i][j]);
print_result(s,s[i][j]+1,j);
cout<<")";
}
}
int main()
{
//6个矩阵,规模为6,N=7
DataType p[N]={30,35,15,5,10,20,25};
int n=N-1;
//
DataType m[N][N]={0};
DataType s[N][N]={0};
Min_ChainOrder_down2up(p,n,m,s);
for(int i=0;i<N;i++)
{
for(int j=0;j<N;j++)
{
cout<<"m["<<i<<"]["<<j<<"]="<<m[i][j]<<endl;
cout<<"s["<<i<<"]["<<j<<"]="<<s[i][j]<<endl;
}
}
cout<<"result:"<<endl;
print_result(s,1,n);
return 0;
}