【算法设计】矩阵乘法

研究生课程系列文章参见索引《在信科的那些课

题目

设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://blog.csdn.net/xiaowei_cqu 未经允许请勿用于商业用途)





没有更多推荐了,返回首页