首页 > 代码库 > 动态规划之矩阵连乘讲解
动态规划之矩阵连乘讲解
注：当时在学矩阵连乘的时候，在网上发现这篇文章总结的不错，便摘抄了下来，今天与大家共享。同时望原创不要见怪！
动态规划之矩阵连乘
给定n个矩阵{A1,A2,…,An}，其中Ai与Ai+1是可乘的，i=1,2,…,n-1。考察这n个矩阵的连乘积A1A2…An。由于矩阵乘法满足结合律，故计算矩阵的连乘积可以有许多不同的计算次序，这种计算次序可以用加括号的方式来确定。若一个矩阵连乘积的计算次序完全确定，则可以依此次序反复调用2个矩阵相乘的标准算法（有改进的方法，这里不考虑）计算出矩阵连乘积。若A是一个p×q矩阵，B是一个q×r矩阵，则计算其乘积C=AB的标准算法中，需要进行pqr次数乘。
矩阵连乘积的计算次序不同，计算量也不同，举例如下：
先考察3个矩阵{A1,A2,A3}连乘，设这三个矩阵的维数分别为10×100，100×5，5×50。若按（（A1A2）A3）方式需要的数乘次数为10×100×5＋10×5×50＝7500，若按（A1（A2A3））方式需要的数乘次数为100×5×50＋10×100×50＝75000。
下面使用动态规划法找出矩阵连乘积的最优计算次序。
1， 设矩阵连乘积AiAi+1…Aj简记为A[i:j]，设最优计算次序在Ak和Ak+1之间断开，则加括号方式为：
（（AiAi+1…Ak）（Ak+1…Aj））
则依照这个次序，先计算A[i:k]和A[K+1:j]然后再将计算结果相乘，计算量是：
A[i:k]的计算量加上A[K+1:j]的计算量再加上它们相乘的计算量。
问题的一个关键是：计算A[i:j]的最优次序所包含的两个子过程（计算A[i:k]和A[K+1:j]）也是最优次序。
2， 设计算A[i:j]所需的最少数乘次数为m[i][j]。
i=j时为单一矩阵，则m[i][i]=0，
i<j时，设最优计算次序在Ak和Ak+1之间断开，则m[i][j]=m[i][k]+m[k+1][j]+pipk+1pj+1，其中p表示数组的维数，例如A0到A5共6个数组（为了C语言的描述方便，下标从0开始），他们表示如下：
//p[0]:第一个矩阵的行数
//p[1]:第一个矩阵的列数，第二个矩阵的行数
//p[2]:第二个矩阵的列数，第三个矩阵的行数
k此时并未确定，需要从i到j-1遍历以寻找一个最小的m[i][j]。我们把这个最小的k放在s[i][j]。
以下是完整实现代码，以一个具体的例子实现，稍加修改即可通用。
#include <iostream>
using namespace std;
//////////////////////////////////////////////////////////////////////////////
//MatrixChain计算m[i][j]所需的最少数乘次数
//并记录断开位置s[i][j]
//////////////////////////////////////////////////////////////////////////////
void MatrixChain(int *p,int n,int **m,int **s)
{
for(int i=0;i<n;i++)
m[i][i]=0;//单个矩阵相乘，所需数乘次数为0
//以下两个循环是关键之一，以6个矩阵为例(为描述方便，m[i][j]用ij代替)
//需按照如下次序计算
//01 12 23 34 45
//02 13 24 35
//03 14 25
//04 15
//05
//下面行的计算结果将会直接用到上面的结果。例如要计算14，就会用到12，24；或者13，34等等
for(int r=1;r<n;r++)
{
for(int i=0;i<n-r;i++)
{
int j=i+r;
//首先在i断开，即（Ai*（Ai+1...Aj））
m[i][j]=m[i][i]+m[i+1][j]+p[i]*p[i+1]*p[j+1];
s[i][j]=i;
for(int k=i+1;k<j;k++)
{
//然后在k（从i+1开始遍历到j-1）断开，即（（Ai...Ak)*（Ak+1...Aj））
int t=m[i][k]+m[k+1][j]+p[i]*p[k+1]*p[j+1];
if(t<m[i][j])//找到更好的断开方法
{
m[i][j]=t;//记录最少数乘次数
s[i][j]=k;//记录断开位置
}
}
}
}
//如果使用下面注释的循环，则是按照如下次序计算
//01 02 03 04 05
//12 13 14 15
//23 24 25
//34 35
//45
//当要计算时14，会用到12，24，而此时24并没有被计算出来。
/*
for(int i=0;i<n;i++)
{
for( int j=i+1;j<n;j++)
{
m[i][j]=m[i][i]+m[i+1][j]+p[i]*p[i+1]*p[j+1];
s[i][j]=i;
for(int k=i+1;k<j;k++)
{
int t=m[i][k]+m[k+1][j]+p[i]*p[k+1]*p[j+1];
if(t<m[i][j])
{
m[i][j]=t;
s[i][j]=k;
}
}
}
}
*/
}
//////////////////////////////////////////////////////////////////////////////
//Traceback打印A[i:j]的加括号方式
//////////////////////////////////////////////////////////////////////////////
void Traceback(int i,int j,int **s)
{
//s[i][j]记录了断开的位置，即计算A[i:j]的加括号方式为：
//(A[i:s[i][j]])*(A[s[i][j]+1:j])
if(i==j)return;
Traceback(i,s[i][j],s);//递归打印A[i:s[i][j]]的加括号方式
Traceback(s[i][j]+1,j,s);//递归打印A[s[i][j]+1:j]的加括号方式
//能走到这里说明i等于s[i][j]，s[i][j]+1等于j
//也就是说这里其实只剩下两个矩阵，不必再分了
cout<<"A"<<i<<"和A"<<(s[i][j]+1)<<"相乘"<<endl;
}
int _tmain(int argc, _TCHAR* argv[])
{
int n=6;//矩阵的个数
int *p=new int[n+1];
//p[0]:第一个矩阵的行数
//p[1]:第一个矩阵的列数，第二个矩阵的行数
//p[2]:第二个矩阵的列数，第三个矩阵的行数
p[0]=30;
p[1]=35;
p[2]=15;
p[3]=5;
p[4]=10;
p[5]=20;
p[6]=25;
int **m,**s;
m=new int*[n];
for( int i=0;i<n;i++)
m[i]=new int[n];
s=new int*[n];
for(int i=0;i<n;i++)
s[i]=new int[n];
MatrixChain(p,n,m,s);
Traceback(0,n-1,s);
for(int i=0;i<n;i++)
{
delete []m[i];
m[i]=NULL;
delete []s[i];
s[i]=NULL;
}
delete []m;
m=NULL;
delete []s;
s = NULL;
delete []p;
p = NULL;
return 0;
}
打印结果是：
A1和A2相乘
A0和A1相乘
A3和A4相乘
A3和A5相乘
A0和A3相乘
实际上要表达的是如下加括号方式：
（（A0（A1A2））（（A3A4）A5））
加了括号之后用第一个来代替，例如（A1A2）可看作A1，这个结果的数乘次数是15125。
具体实例见：
NYOJ 460 项链
NYOJ 536 开心的mdd