动态规划之矩阵连乘

以下内容参考（摘抄）《算法设计与分析》，王晓东编著，清华大学出版社2003年1月第1版。

给定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*
;

for( int i=0;i<n;i++)

m[i]=new int
;

s=new int*
;

for(int i=0;i<n;i++)

s[i]=new int
;

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。

本文来自CSDN博客，转载请标明出处：http://blog.csdn.net/liguisen/archive/2008/03/08/2158127.aspx