【动态规划】矩阵连乘问题 C++（附代码实例和复杂度分析）

问题描述

给定n个矩阵{A1,A2,…An}，其中Ai和Ai+1是可乘的。对于这n个矩阵的连乘积，可能有不同的计算次序，如何确定计算矩阵连乘积的计算次序，使得需要的数乘次序最少

问题分析

由于可以在矩阵链的任意位置加括号改变矩阵连乘的运算次序，使用穷举法会面临指数级的复杂度，不能算一个有效的算法。究其根本，是在穷举过程中进行了大量的重复运算，比如前面已经算过A1xA2了，后面计算A1xA2xA3又要计算一次A1xA2，这些重复计算耗费了大量时间，对此可以采用动态规划的方法

使用动态规划算法可按以下几个步骤进行

1. 分析最优解的结构

将矩阵A1,…,An的连乘积简记为A[1:n]，则该问题可拆分成两个子问题，即求A[1:k]和A[k+1:n]的解，最后再将结果相乘。如果A[1:n]为最优解，那么它的两个子问题A[1:k]和A[k+1:n]也一定是最优解，否则用它们的最优解替换当前解，能够得到更优的A[1:n]，这显然与A[1:n]为最优解冲突

因此，矩阵连乘问题的最优解包含其子问题的最优解，这种性质称为最优子结构性质。

2. 建立递推方程

设计算A[i:j]的最少数乘次数为m[i:j]

可以得到如下递归方程
$$m[i][j]=\{\begin{array}{ll}0& \text{if }i=j,\\ {\min }_{i\le k<j}\left(m[i][k]+m[k+1][j]+p_{i-1}\times p_k\times p_j\right)& \text{if }i<j.\end{array}$$

这个方程的意思是，将矩阵 Ai,…,Aj 在Ak处将矩阵链一分为二，左边进行连乘得到的m[i][k]加上右边连乘得到的m[k+1][j]，最后加上$p_{i-1}\times p_k\times p_j$（左边连乘得到的矩阵和右边连乘得到的矩阵相乘所需要的乘法次数）就是整个连乘的总次数。

注：对于维度分别为$p_i\times p_j$和$p_j\times p_k$的两个矩阵相乘，所需要的乘法次数是$p_i\times p_j\times p_k$次。

如果将对应m[i:j]断开的位置记为s[i:j]，在计算出最优值m[i:j]后，可以递归的由s[i:j]构造最优解

3. 计算最优值

现在考虑不重复的矩阵连乘子问题有多少种情况
1个矩阵连乘，有n种情况
2个矩阵连乘，有n-1种情况（12、23、34、…(n-1)n）
3个矩阵连乘，有n-2种情况（123、234、345、…、(n-2)(n-1)n）
4个矩阵连乘，有n-3种情况
…
n-1个矩阵连乘，有2种情况
n个矩阵连乘，有1种情况
所以不重复的连乘情况有$n^2$种

以如下6个矩阵连乘为例

• 1个矩阵连乘

  • m[1]=0
  • m[2]=0
  • m[3]=0
  • m[4]=0
  • m[5]=0
  • m[6]=0

• 2个矩阵连乘

  • m[1:2]=m[1]+m[2]+p[0]*p[1]*p[2]=15750
  • m[2:3]=m[2]+m[3]+p[1]*p[2]*p[3]=2625
  • m[3:4]=m[3]+m[4]+p[2]*p[3]*p[4]=750
  • m[4:5]=m[4]+m[5]+p[3]*p[4]*p[5]=1000
  • m[5:6]=m[5]+m[6]+p[4]*p[5]*p[6]=5000

• 3个矩阵连乘

  • m[1:3]=min(m[1]+m[2:3]+p[0]*p[1]*p[3]，m[1:2]+m[3]+p[0]*p[2]*p[3])=7875
  • m[2:4]=min(m[2]+m[3:4]+p[1]*p[2]*p[4]，m[2:3]+m[4]+p[1]*p[3]*p[4])=4375
  • m[3:5]=2500
  • m[4:6]=3500

4、5、6个矩阵连乘同理

• 4个矩阵连乘

  • m[1:4]=9375
  • m[2:5]=7125
  • m[3:6]=5375

• 5个矩阵连乘

  • m[1:5]=11875
  • m[2:6]=10500

• 6个矩阵连乘

  • m[1:6]=15125

代码
整体思路是利用循环去填充表中的每一斜行，最后得到的m[1][n]就是最优解
在MatrixChain函数中，首先将m[i][i]位置都置0，然后分别利用k控制外层循环，j控制内层循环

• 外层循环：控制当前i和j的差值，即当前几个矩阵做连乘（k+1个）
• 内层循环：遍历当前连乘个数下的所有可能情况（填充表中的斜行）

find_optimal_parens函数用于寻找当前加括号方法的最优解（min），并将这个方法存储到s数组中

void MatrixChain(int *p, int n) {
	for(int i=1;i<=n;i++)
		m[i][i] = 0;
	for (int k = 1; k < n; k++) {
		for (int j = k + 1, i = j - k; j <= n; j++,i++) {
			find_optimal_parens(i, j);
		}
	}
}

void find_optimal_parens(int i, int j) {
	for (int k = i; k < j; k++) {//在k处断开
		if (k==i || m[i][j] > m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j]) {
			m[i][j] = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];
			s[i][j] = k;
		}
	}
	return;
}

4. 构造最优解

前面已经计算出最优解了，并且将最优情况记录在了s数组里，s数组如下

只要递归的查找s数组，就能构造出最优情况

比如，s[1][6]=3，说明A[1:6]是在A3处断开 (A1A2A3)(A4A5A6);
s[1][3]=1，说明A[1:3]是在A1处断开(A1)(A2A3);
s[4][6]=5，说明A[4:6]是在A5处断开(A4A5)(A6)

合起来就得到了最终的结果（(A1)(A2A3))((A4A5)(A6)）

代码

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

完整实例

#define _CRT_SECURE_NO_WARNINGS
#include<iostream>
using namespace std;
int m[7][7];
int s[7][7];
int p[7] = { 30,35,15,5,10,20,25 };
void find_optimal_parens(int i, int j) {
	for (int k = i; k < j; k++) {//在k处断开
		if (k==i || m[i][j] > m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j]) {
			m[i][j] = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];
			s[i][j] = k;
		}
	}
	return;
}
void MatrixChain(int* p, int n) {
	for (int i = 1; i <= n; i++)
		m[i][i] = 0;
	for (int k = 1; k < n; k++) {
		for (int j = k + 1, i = j - k; j <= n; j++, i++) {
			find_optimal_parens(i, j);
		}
	}
}
void Traceback(int i, int j) {
	if (i == j) {
		cout << "A" << i;
		return;
	}
	cout << "(";
	Traceback(i, s[i][j]);
	Traceback(s[i][j] + 1, j);
	cout << ")";
}
int main() {
	int n = 6;
	MatrixChain(p, n);
	cout << m[1][n] << endl;
	Traceback(1, n);
	return 0;
}

运行结果

复杂度分析

算法MatrixChain的主要计算量取决于三重循环。循环体内的计算量为O(1)，而三重循环的总次数为$O(n^3)$，因此该算法的时间上界为$O(n^3)$。占用的空间为$O(n^2)$。显然比穷举法更加有效