【算法设计】矩阵乘法

题目

设A1,A2,…,An为矩阵序列，Ai为Pi-1×Pi阶矩阵，i = 1,2,…,n. 确定乘法顺序使得元素相乘的总次数最少.
输入：向量P = <P0, P1, … , Pn>

实例：

P = <10, 100, 5, 50> A1: 10 × 100, A2: 100 × 5, A3: 5 × 50

乘法次序：

(A1 A2)A3: 10 × 100 × 5 + 10 ×5 × 50 = 7500

A1(A2 A3): 10 × 100 × 50 + 100 × 5 × 50 = 75000

 

 

搜索空间的规模

先将矩阵链加括号分为两部分，即P=A1*A2*...*An=（A1*A2...*Ak）*（Ak+1*...An），则有f(n)=f(1)*f(n-1)+f(2)*f(n-2)+...+f(n-1)*f(1)种方法。f(n)为一个Catalan数，所以一般的方法要计算种。

动态规划算法

输入P=< P0, P1, …, Pn>，Ai..j 表示乘积 AiAi+1…Aj 的结果，其最后一次相乘是：

m[i,j] 表示得到Ai..j的最少的相乘次数。

递推方程：

为了确定加括号的次序，设计表s[i,j],记录求得最优时最一位置。

算法递归实现

由上面的递归公式，很容易得到算法的递归实现：

const int N=5;
int m[N][N]; //m[i][j]存储Ai到Aj的最小乘法次数
int s[N][N];//s[i][j]存储Ai到Aj之间加括号的位置

int RecurMatrixChain(int P[],int i,int j)
{
	m[i][j]=100000;
	s[i][j]=i;
	if(i==j)
		m[i][j]=0;
	else{
		for(int k=i;k<j;k++){
			int q=RecurMatrixChain(P,i,k)+RecurMatrixChain(P,k+1,j)+P[i]*P[k+1]*P[j+1];
			if(q<m[i][j]){
				m[i][j]=q;
				s[i][j]=k;
			}
		}
	}
	return m[i][j];
}

int main()
{
	int P[N+1]={30,35,15,5,10,20};
	for(int i=0;i<N;i++)
		m[i][i]=0;
	m[0][N-1]=RecurMatrixChain(P,0,N-1);
	return 0;
}

递归实现的复杂性

复杂性满足递推关系：

由数学归纳法可得：

可见递归实现的复杂性虽然较一般算法有改进，但还是较高。分析原因，主要是子问题重复程度高。如下图所示：

1..4表示计算Ai..j中i=1,j=4的子问题，其子问题包括A1..1，而A1..2，A1..3中都包括子问题A1..1，所以很多子问题被重复计算了多次。

于是，我们想到用自底向上的迭代实现。

 

算法迭代实现

迭代实现主要思想是子问题由小到大，每个子问题只计算一次，并且把结果保存起来，后来用到这个子问题时，直接代入。

void MatrixChain(int P[],int n)
{
	int r,i,j,k,t;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			m[i][j]=0;
	//r为当前计算的链长（子问题规模）
	for(r=2;r<=n;r++){  
		//n-r+1为最后一个r链的前边界
		for(i=0;i<n-r+1;i++){
			//计算前边界为r，链长为r的链的后边界
			j=i+r-1;
			//将链ij划分为A(i) * ( (A(i+1) ... A(j) )
			m[i][j]=m[i+1][j]+P[i]*P[i+1]*P[j+1];
			//记录分割位置
			s[i][j]=i;
			for( k=i+1;k<j-1;k++){
				//将链ij划分为( A(i)...A(k) )* ( (A(k+1) ... A(j) )
				t=m[i][k]+m[k+1][j]+P[i]*P[i+1]*P[j+1];
				if(t<m[i][j]){
					m[i][j]=t;
					s[i][j]=k;
				}
			}
		}
	}
}

int main()
{
	int P[N+1]={30,35,15,5,10,20};
	MatrixChain(P,N);
}

迭代实现的复杂性

行7,9,16的循环为O(n)，外层循环为O(1)，所以算法复杂度W(n)=O(n^3)

迭代过程的一个实例

子问题由小到大的计算过程如下图所示：

结果打印

再写一个打印结果，以及打印优化函数备忘录m和标记函数的s的函数：

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

void PrintMS(int m[][N],int s[][N],int N)
{
	for(int r=0;r<N;r++){
		for(int i=0;i<N-r;i++){
			int j=i+r;
			cout<<"m["<<i+1<<","<<j+1<<"]="<<m[i][j]<<"\t";
		}
		cout<<endl;
	}
	for(int r=1;r<5;r++){
		for(int i=0;i<N-r;i++){
			int j=i+r;
			cout<<"s["<<i+1<<","<<j+1<<"]="<<s[i][j]+1<<"\t";
		}
		cout<<endl;
	}
}

一个简单的测试实例

用一个N=5，P=<30,35,15,5,10,20>的简单实例，运行上述代码：

*注意：上述代码与解释中的下标不同，即代码中s[i-1][j-1]表示实际中的s[i,j]

代码下载：http://download.csdn.net/detail/xiaowei_cqu/4764899

参考资料：屈婉玲 刘田等 《算法设计与分析》

（转载请注明作者和出处：http://blog.csdn.net/xiaowei_cqu未经允许请勿用于商业用途）