矩阵连乘

参考资料：《计算机算法设计与分析》第3版：Page49，王晓东 著

一、问题描述：

给定n个矩阵{A₁, A₂, ..., A_n}，其中A_i与A_i+1是可乘的，i=1, 2, ..., n-1。考察n个矩阵的连乘积。

二、完全加括号：

矩阵乘法满足结合律，故计算矩阵的连乘积可以有很多不同的计算次序，这种计算次序可以用加括号的方式来确定。若一个矩阵的连乘积的计算次序完全确定，也就是说该连乘积已经完全加括号，则可以依此次序反复调用2个矩阵相乘的标准算法计算出矩阵连乘积。完全加括号的矩阵连乘积可递归的定义为：

1. 单个矩阵是完全加括号的；

2. 矩阵连乘积A是完全加括号的，则A可以表示为2个完全加括号的矩阵的连乘积B和C的乘积并加括号，即A=(BC)。

每一种完全加括号方式对应于一种矩阵连乘积的计算次序，而矩阵连乘积的计算次序与其计算量有密切关系。

三.两个矩阵乘积所需的计算量：A_p×q, B_q×r, C_p×r=AB: pqr

两个矩阵乘积的标准算法如下，其中 ra, ca 和 rb, cb 分别表示矩阵A和B的行数和列数。

void matrix_multiply(int **a,int **b,int **c,int ra,int ca,int rb,int cb){
	if(ca!=rb){
		cout<<"can not multiply";
		return;
	}
	int i,j,k;
	for(i=1;i<=ra;++i){
		for(j=1;j<=cb;++j){
			c[i][j]=0;
			for(k=1;k<=ca;++k){
				c[i][j]+=a[i][k]*b[k][j];
			}
		}
	}
	return;
}

穷举法不同计算次序（指数级）：

P(n)=1, n=1;

P(n)=P(1)P(n-1)+P(2)P(n-2)+...+P(n-1)P(1)+, n>1;

四.最优子结构：

记矩阵连乘积A_iA_i+1...A_j简记为A[i:j]，考察A[1:n]的最优计算次序。设这个计算次序在矩阵A_k与A_k+1之间将矩阵链断开，1<=k<n，则相应加括号方式为((A_i...A_k)(A_k+1...A_n))。依次次序，先计算A[1:k]和A[k+1:n]，然后将计算结果相乘得A[1:n]，依此次序总计算量为A[1:k]和A[k+1:n]计算量之和，再加上A[1:k]和A[k+1:n]相乘的计算量。此问题的关键特征是：计算A[1:n]的最优次序所包含的计算矩阵子链A[1:k]和A[k+1:n]的次序也是最优的。

五.递归关系：

设计算A[1:n]，1<=i<=j<=n，所需的最小数乘次数为m[i][j]，则原问题的最优值为m[1][n]。

m[i][j]=0, i=j;

m[i][j]=min{m[i][k]+m[k][j]+p_i-1p_kp_j} (i<=k<j), i<j;

m[i][j]给出了最优值，即计算A[i:j]所需的最少数乘次数，同时还确定了计算A[i:j]的最优次序中的断开位置k，将该位置记为s[i][j]，在计算出最优值m[i][j]后，可递归的由s[i][j]构造出相应的最优解。对n个矩阵连乘问题，其不同子问题的个数为C(n,2)+n=n(n-1)/2+n。

六.计算最优值：

自底向上计算。输入参数{p₀, p₁, ..., p_n}存储于数组p中，n记录矩阵连乘链长，记录最优值的数组m，记录最优断开位置的数组s。要计算最优值只需调用matrix_chain(p,n,m,s)。

void matrix_chain(int *p,int n,int **m,int **s){
	int i,j,k,r,t;
	for(i=1;i<=n;++i){
		m[i][i]=0;
	}
	for(r=2;r<=n;++r){
		for(i=1;i<=n-r+1;++i){
			j=i+r-1;
			m[i][j]=m[i+1][j]+p[i-1]*p[i]*p[j];
			s[i][j]=i;
			for(k=i+1;k<j;++k){
				t=m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];
				if(t<m[i][j]){
					m[i][j]=t;
					s[i][j]=k;
				}
			}
		}
		cout<<endl;
	}
}

七.构造最优解：

按求解最优值时记录的最优断点矩阵s指示的加括号方式输出计算A[i:j]的最优计算次序。要得到A[1:n]的最优计算次序只需调用matrix_solution(1,n,s)。

void matrix_solution(int i,int j,int **s){
	if(i==j){
		cout<<"A"<<i;
		return;
	}
	cout<<"(";
	matrix_solution(i,s[i][j],s);
	matrix_solution(s[i][j]+1,j,s);
	cout<<")";
}

八.测试代码：

#include<iostream>

using namespace std;
void matrix_chain(int *p,int n,int **m,int **s);
void matrix_solution(int i,int j,int **s);
int main(int argc,char *argv[]){
    int *p;
    int **m;
    int **s;
    int n;
    
    cin>>n;
    p=new int[n+1];
    m=new int*[n+1];
    s=new int*[n+1];
    for(int i=1;i<=n;++i){
        m[i]=new int[n+1];
        s[i]=new int[n+1];
    }
    
    for(int i=0;i<=n;++i){
        cin>>p[i];
    }

    matrix_chain(p,n,m,s);
    cout<<"minmum multiply operation counts:"<<m[1][n]<<endl;
    matrix_solution(1,n,s);
    
    for(int i=1;i<=n;++i){
        delete[] m[i];
        delete[] s[i];
    }
    delete[] p;
    delete[] m;
    delete[] s;
    return 0;
}