基于动态规划的矩阵连乘最优方法

问题描述：

在科学计算中经常要计算矩阵的乘积。矩阵A和B可乘的条件是矩阵A的列数等于矩阵B的行数。若A是一个p×q的矩阵，B是一个q×r的矩阵，则其乘积C=AB是一个p×r的矩阵。其标准计算公式为：

由该公式知计算C=AB总共需要pqr次的数乘。 为了说明在计算矩阵连乘积时加括号方式对整个计算量的影响，我们来看一个计算3个矩阵{A1，A2，A3}的连乘积的例子。设这3个矩阵的维数分别为10×100，100×5和5×50。若按第一种加括号方式((A1A2)A3)来计算，总共需要10×100×5+10×5×50=7500次的数乘。若按第二种加括号方式(A1(A2A3))来计算，则需要的数乘次数为100×5×50+10×100×50=75000。第二种加括号方式的计算量是第一种加括号方式的计算量的10倍。由此可见，在计算矩阵连乘积时，加括号方式，即计算次序对计算量有很大影响。

所以，对于给定的相继n个矩阵{A1,A2,…,An}(其中Ai的维数为pi-1×pi ，i=1,2,…,n)，如何确定计算矩阵连乘积A1A2…An的一个计算次序(完全加括号方式)，使得依此次序计算矩阵连乘积需要的数乘次数最少，即所谓的矩阵连乘积的最优计算次序问题。

思路：

采用动态规划：动态规划算法的基本思想是将待求解问题分成若干个子问题，先求解子问题，然后从这些子问题的解得到原问题的解。与分治法不同的是，动态规划法经分解得到的子问题往往不是相互独立的，前一子问题的解为后一子问题的解提供有用的信息，可以用一个表来记录所有已解决的子问题的答案，不管该子问题以后是否被用到，只要它被计算过，就将其结果填入表中。

设a[n][m]为第n个矩阵到第m个矩阵连乘的最小乘法数（n, m >= 1），如果n=m，则a[n][m]这段中就一个矩阵，需要计算的次数为0。b[i-1], b[i]分别为第i个矩阵的行数和列数（i >=1），那么：

1. a[n][n + 1]易求，为相邻两个矩阵相乘的乘法数，即b[n-1] * b[n] * b[n + 1]；

2. An ~ Am可以任意拆分为An ~ Ai及Ai+1 ~ Am两部分相乘（n <= i <= m），则a[n][m]为所有拆分情况中乘法次数最少的一种，即a[n][m] = min(a[n][i] + a[i+1][m] + b[n - 1] * b[i] * b[m])，相当于将其分成两个分部，用三个步骤完成矩阵连乘。a[n][i]表示所拆分的第1部分的最少的乘法次数，a[i+1][m]为第2部分最少的乘法次数，b[n - 1] * b[i] * b[m]则是上述两者结果的乘法次数。

综上可以构造递归关系如下：

代码实现时需要注意的问题：计算顺序！！！  因为要保证在计算a[n][m]查找a[n][i]和a[i+1][m]的时候，a[n][i]和a[i+1][m]已经计算出来了。

所以坐标的关系图：

从上述的递归式可以发现解的过程中会有很多重叠子问题，可以用一个nXn维的辅助表m[n][n] s[n][n]分别表示最优乘积代价及其分割位置k 。 辅助表s[n][n]可以由2种方法构造，一种是自底向上填表构建，该方法要求按照递增的方式逐步填写子问题的解，也就是先计算长度为2的所有矩阵链的解，然后计算长度3的矩阵链，直到长度n；另一种是自顶向下填表的备忘录法，该方法将表的每个元素初始化为某特殊值(本问题中可以将最优乘积代价设置为一极大值)，以表示待计算，在递归的过程中逐个填入遇到的子问题的解。

对于 p={30 35 15 5 10 20 25}

第一种：自底向上法

自底向上法过程：

由p可以知道各个矩阵的行列情况：

计算过程：

通过上述表格可以获得下面的结果：

以下进行C++的编程实现：

#include<iostream>
using namespace std;

const int N=7;
//p为矩阵链，p[0],p[1]代表第一个矩阵，p[1],p[2]代表第二个矩阵，length为p的长度
//所以如果有六个矩阵，length=7，m为存储最优结果的二维矩阵，s为存储选择最优结果路线的
//二维矩阵
//nXn维的辅助表m[n][n] s[n][n]分别表示最优乘积代价及其分割位置k
void MatrixChainOrder(int *p,int m[N][N],int s[N][N],int length)
{
	int n=length-1;//矩阵个数
	int l,i,j,k,q=0;
	//m[i][i]只有一个矩阵，所以相乘次数为0，即m[i][i]=0;
	for(i=1;i<length;i++)
	{
		m[i][i]=0;//对角线位置
	}
	//l表示矩阵链的长度
	// l=2时，计算 m[i,i+1],i=1,2,...,n-1 (长度l=2的链的最小代价)
	for(l=2;l<=n;l++)//l=2为对角线紧邻的右上位置，随着l的增大，即沿着右上角的方向递增，扩大链的长度
	{
		for(i=1;i<=n-l+1;i++)
		{
			j=i+l-1; //以i为起始位置，j为长度为l的链的末位，
			m[i][j]=0x7fffffff;
			//k从i到j-1,以k为位置划分
			for(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;
				}
			}
		}
	}
	cout << m[1][N-1] << endl;
}
//以递归形式调用，添加括号，注意理解！
void PrintAnswer(int s[N][N],int i,int j)
{
	if(i==j)//只有一个矩阵，直接输出
	{
		cout<<"A"<<i;
	}
	/*else if(i+1==j)
	{
		cout<<"(A"<<i<<"A"<<j<<")";
	}*/
	else//是否需要再添加一种情况：两个矩阵，加括号输出？？
	{
		cout<<"(";
		PrintAnswer(s,i,s[i][j]);//递归，从得到最优解的地方开始s[i][j]处断开
		PrintAnswer(s,s[i][j]+1,j);
		cout<<")";
	}
}
int main()
{
	int p[N]={30,35,15,5,10,20,25};
	int m[N][N],s[N][N];
	MatrixChainOrder(p,m,s,N);
	PrintAnswer(s,1,N-1);
	return 0;
}

运行结果如下：

第二种：自顶向下 （备忘录）法

#include<iostream>
using namespace std;

const int N=7;
//备忘录法（自顶向下 ）
int MatrixChainOrder2(int *p,int m[N][N],int s[N][N],int i, int j)
{

	if (i == j)
	{
		return 0;
	}
	if (m[i][j] < 0x7fffffff)//是否已经计算，以计算后的值对初始值进行替换
	{
		return m[i][j];
	}

	for (int k=i; k<j; ++k)
	{

		int tmp = MatrixChainOrder2(p,m,s,i,k) + MatrixChainOrder2(p,m,s,k+1,j) + p[i-1]*p[k]*p[j];
		if (tmp < m[i][j])
		{
			m[i][j] = tmp;//从k处断开，分别求得每次的数乘次数
			s[i][j] = k;
		}
	}
	return m[i][j];
}


void PrintAnswer(int s[N][N],int i,int j)
{
	if(i==j)
	{
		cout<<"A"<<i;
	}
	else
	{
		cout<<"(";
		PrintAnswer(s,i,s[i][j]);
		PrintAnswer(s,s[i][j]+1,j);
		cout<<")";
	}
}


int main()
{
	int p[N]={30,35,15,5,10,20,25};
	int m[N][N],s[N][N];

	for (int i=0; i<N; ++i)
	{
		for (int j=0; j<N; ++j)
			m[i][j] = 0x7fffffff;//初始化为最大值
	}
	cout << MatrixChainOrder2(p,m,s,1,N-1)<<endl;
	PrintAnswer(s,1,N-1);
	return 0;
}

运行结果是一样的。

参考：

http://blog.csdn.net/sixtyfour/article/details/12250415

http://www.cnblogs.com/liushang0419/archive/2011/04/27/2030970.html

http://www.cnblogs.com/scarecrow-blog/p/3712580.html

http://www.warmtoyou.com/articles/detailArticlePub.htm?articleId=81