关于动态规划解决矩阵连乘问题

关于动态规划解决矩阵连乘问题

以《计算机算法设计与分析 第5版》教材中的例子为例。

动态规划解决具体问题分为4步：

1. 分析最优解结构。

最优子结构是：问题的最优解包含着子问题的最优解。

2. 建立递归关系。

递归的定义最优值。

3. 计算最优值。

根据递归式，写出递归算法。

4. 构造最优解。

根据算法获取到的信息构造最优解。

空谈比较空泛，现在利用动态规划寻找矩阵连乘的最优顺序：
例：计算矩阵连乘积A₁A₂A₃A₄A₅A₆的最优计算顺序。
其中各矩阵维度分别是
A₁:30×35
A₂:35×15
A₃:15×5
A₄: 5×10
A₅:10×20
A₆:20×25

分析：如果直接莽，从左乘到右，将要计算30×35×15×5×10×20×25次乘法。（3.9亿多次）

再稍微讲下这个矩阵连乘的乘法次数的计算，例如：A₁×A₂，其中A₁是p×r的矩阵，A₂是r×q的矩阵。则结果A₀应是一个p×q的矩阵。
一共要进行p×r×q次乘法计算，因为A₀的每一项都要进行r次乘法，例如A₀[a][b]=${\sum }_{i=1}^n$A₁[a,i]×A₂[i,b]。
而A₀共有p×q项，所以需要计算p×r×q次乘法。

所以不能直接莽，要靠动态规划。
方便起见，将矩阵连乘A_iA_i+1…A_j简记为A[i:j]，则原问题就是计算A[1:n]的最优解。
将矩阵维度存储在p[7]数组中，其中p₀=A₁的行，p₁=A₁的列，p₂=A₂的列，p₃=A₃的列……p₆=A₆的列。
①A[1:n]的最优乘法次序包含了A[1:k]与A[k+1:n]的最优次序，因为矩阵乘法是不能交换位置的（不满足交换律），所以矩阵连乘满足最优子结构性质。

②设一个数组m[6][6] (这里不是从0到5，而是从1到6) 用m[i][j]表示A[i:j]的最少乘法次数。
根据最优子结构和矩阵乘法的乘法次数计算方法，可以得出下列m[i][j]的递归定义：
当i=j时，m[i][j]=0.这是很显然的，单一矩阵无需计算。
当i≠j时，m[i][j]=min{ m[i][k]+m[k+1][j]+p_i-1p_kp_j }.这是把m[i][j]分为左右两个矩阵连乘，再用左右两个矩阵计算总的乘法次数。p_i-1p_kp_j就是左边矩阵的行×左边矩阵的列×右边矩阵的列。

③递归关系已经建立，接下来就是设计算法计算最优值，也就是计算m[6][6]数组。在写代码之前，先要弄明白m[][]数组的计算规则
根据递推式，应当沿对角线层层推进式的计算：（网上找的一张同题的图）把 a[][]数组看为m[][]数组就可以了。

例如m[2][5]的计算：
m[2][5]有3种划分方式，

• ①划分为m[2][2]与m[3][5]。
• ②划分为m[2][3]与m[4][5]。
• ③划分为m[2][4]与m[5][5]。

按照递推式分别计算3种情况，取最小值。

递推式为：
当i=j时，m[i][j]=0.
当i≠j时，m[i][j]=min{ m[i][k]+m[k+1][j]+p_i-1p_kp_j }.
p₀=30，p₁=35，p₂=15，p₃=5，p₄=10，p₅=20，p₆=25

第①种情况：m[2][5]=m[2][2]+m[3][5]+p₁p₂p₅=0+2500+35×15×20=13000
第②种情况：m[2][5]=m[2][3]+m[4][5]+p₁p₃p₅=2625+1000+35×5×20=7125
第③种情况：m[2][5]=m[2][4]+m[5][5]+p₁p₄p₅=4375+0+35×10×20=11375

很显然，第②种情况的m[2][5]最低，所以应当将m[2][5]赋值为7125。
第②种情况是将m[2][5]分为m[2][3]与m[4][5]，断点k=3,可以另起一个数组s[][]存储断点，那么s[2][5]赋值为3。

好了，计算规则清楚了，再来看书上的代码：

void MatrixChain(int *p, int n, int **m, int **s){
	for(int i=1, i<=n; i++)
		m[i][i]=0;//这个循环就是把对角线置0
	for(int r=2; r<=n ;r++){//r是斜线的层数，对角线就是r=1
		for(int i=1; i<=n-r+1; i++){//i为行号，因为随着斜线层数r增加，i的上限逐渐降低。例如当r=2时，i只需要从1到n-1。
			int j=i+r-1;//j为列号，当r与i确立的时候，实际上坐标已经定位了，j只是把它变为直角坐标系。
						//随着斜线层数r的增加，j的起始位置右移，随着i的增加，j也沿着斜线向右增加。
			m[i][j]=m[i+1][j]+p[i-1]*p[i]*p[j];
			//其实这里省略了一个m[i][i]，但因为m[i][i]==0所以省略了。
			//这里先默认“按照第一种情况将m[i][j]划分为m[i][i]与m[i+1][j]是最优的”下面再用一层循环来计算其他几种情况的乘法次数
			s[i][j]=i;//这里也是默认第一种情况，默认断点在i处
			for(int k=i+1; k<j; k++){//这一层循环就是寻找其他情况的代价，将所有可以将m[i][j]切开的断点k都尝试一遍，然后计算乘法次数t。
				int t=m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];
				if(t<m[i][j]){//最后比较，如果t小于所存储的乘法次数，则更新乘法次数并且更新断点k
					m[i][j] = t;
					s[i][j] = k;
				}
			}
		}
	}
}

④构造最优解，根据算法计算出了最优值后，我们就可着手构造最优解了，最优解的信息存在于s[][]数组中。

若s[i][j]==k，那么这意味着A[i:j]的最优断点在k处，即应当将A[i][j]分为A[i][k]与A[k+1][j]两部分，再分别探讨s[i][k]与s[k+1][j]，寻找A[i][k]与A[k+1][j]的最优断点，一直递归下去直到s[m][n]0（等价于mn）为止。

有了s[][]数组，原理就很简单了，书中代码如下：

void Traceback(int i, int j, int **s){
	if(i==j)
		return;
	Traceback(i, s[i][j], s);//断点为s[i][j]，左边是i到s[i][j]
	Traceback(s[i][j]+1, j, s);//右边是s[i][j]+1到j
	cout<<"Multiply A"<< i <<"," <<s[i][j];
	cout<<" and A "<<(s[i][j]+1)<<","<<j<<endl;
	//最后两句话就是输出A[i:s[i][j]]×A[(s[i][j]+1):j]
}