矩阵连乘

矩阵连乘

题目：

​ 给定n个矩阵，确定计算矩阵连乘积的计算次序（完全加括号的方式），使得依此计算次序计算矩阵连乘积需要的 乘数 次数最少。

分析：

​ 矩阵乘法属于最优解问题，第一反应是是否能用贪心来解决。不能。不会数学证明，但可以举例子，参考链接：

明显贪心算法所得结果不是最优的，故不能用贪心来求解。解决这道问题的正确方法是采用动态规划算法：

1. 带备忘录的自顶向下法
2. 自底向上法

实现一：动态规划（自底向上）

1. 分析最优解的结构

​ 计算 A[i:j] 的最优次序，即 m[i][j] 的值，一定存在一个分界k使得

m[i][j]最小，用数学公式表示为：

m[i][j] = min{ m[i][k] + m[k+1][j] + p[i-1]p[k]p[j] }

​ 证明：m[i][j] 最优时，子结构 m[i][k]、m[k+1][j] 也是最优的。

​ 反证法：当 m[i][j] 最优时，p[i-1]p[k]p[j] 是固定值（常量），倘若存在m[i][k] 或 m[k+1][j]不是最优的，则他们之和得出的m[i][j]也不是最优的，故矛盾，可证明 m[i][j] 最优时，子结构 m[i][k]、m[k+1][j] 也是最优的。

2. 建立递归关系

m [ i ] [ j ] = { 0 , i = j m i n {   m [ i ] [ k ] + m [ k + 1 ] [ j ] + p [ i − 1 ] p [ k ] p [ j ]   } , i < j m[i][j]= \begin{cases} 0, \quad i=j \\ min\{ \space m[i][k] + m[k+1][j] + p[i-1]p[k]p[j] \space \}, \quad i<j \end{cases} m[i][j]={0,i=jmin{ m[i][k]+m[k+1][j]+p[i−1]p[k]p[j] },i<j​

3. 计算最优值

​ 用二维数组中的某个元素表示对应 A[i:j] 的乘法次数，自底向上计算最优值。规模最小的是自己乘自己时间开销为0，表现在矩阵中就是主对角线位置全为0。接着计算距离为1，表现在矩阵上就是主对角线左上方平行的斜线，如图（红色的线）：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-rFoUfNna-1668480940633)(C:\Users\肥肥龙\AppData\Roaming\Typora\typora-user-images\image-20220921111807401.png)]

依次往上，每次填写的都是平行的斜线。在代码方面用r表示斜线序号，i控制行，j控制列，k表示分界点位置，代码如下：

// 矩阵连乘
/*
p: 矩阵的行列
n: 矩阵个数
m: m[i][j] 表示矩阵i到j乘积最小开销
s: s[i][j] 最小断开处
*/
//void matrixChain(int *p, int n, int **m, int **s) {
void matrixChain(int p[], int n, int m[][n], int s[][n]) {
	// m[i][i]处全为0
	for (int i = 0; i < n; i++) {
		m[i][i] = 0;
	}

	// 动态规划，自底向上，先算距离为1的，后依次递增（从主对角线开始斜边向上）
	for	(int r = 1; r < n; r++) { // 斜边
		for (int i = 0; i < n - r; i++) { // 行
			int j = i + r; // 列
			int minCost = INT_MAX;
			int t;  // 临时值
			int minPivot;  // 分界点

			for (int k = i; k < j; k++) { // 寻找最优分界点
				t = m[i][k] + m[k + 1][j] + p[i] * p[k + 1] * p[j + 1];

				if (minCost > t) {
					minCost = t;
					minPivot = k;
				}
			}

			m[i][j] = minCost;
			s[i][j] = minPivot + 1;  // 加1是为了表示矩阵的编号

		}
	}
}

用教材上的例子：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-roZazI31-1668480940635)(C:\Users\肥肥龙\AppData\Roaming\Typora\typora-user-images\image-20220921113642511.png)]

得出的m和s矩阵如下图：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-H0niLkgY-1668480940635)(C:\Users\肥肥龙\AppData\Roaming\Typora\typora-user-images\image-20220921114058245.png)]

4. 构造最优解

​ 根据计算出来的s矩阵递归计算最优解，递归方程如下*（可能有误，请指正）*：
$$A[i:j]=\{\begin{array}{l}A[i],i=j\\ A[i:s[i][j]]\ast A[s[i][j]+1:j],i<j\end{array}$$
代码如下：

// 递归分治的方法找到最优结合方法
/*
	start: 开始矩阵
	end: 结束矩阵
	s: 动态规划中求得的分界点矩阵
*/
void traceBack(int start, int end, int s[][6]) {
	if (start == end) {
		printf("%d", start);
		return;
	}

	int pivot = s[start - 1][end - 1];
	printf(" ( ");
	traceBack(start, pivot, s);  // 递归左边
	printf(" ) ");
	printf(" ( ");
	traceBack(pivot + 1, end, s); // 递归右边
	printf(" ) ");
}

结果如下：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-bM40gV6I-1668480940636)(C:\Users\肥肥龙\AppData\Roaming\Typora\typora-user-images\image-20220921115828044.png)]

5. 完整代码

#include <stdio.h>
#define INT_MAX 2147483647

// 矩阵连乘
/*
p: 矩阵的行列
n: 矩阵个数
m: m[i][j] 表示矩阵i到j乘积最小开销
s: s[i][j] 最小断开处
*/
//void matrixChain(int *p, int n, int **m, int **s) {
void matrixChain(int p[], int n, int m[][n], int s[][n]) {
	// m[i][i]处全为0
	for (int i = 0; i < n; i++) {
		m[i][i] = 0;
	}
	// 动态规划，自底向上，先算距离为1的，后依次递增（从主对角线开始斜边向上）
	for	(int r = 1; r < n; r++) { // 斜边
		for (int i = 0; i < n - r; i++) { // 行
			int j = i + r; // 列
			int minCost = INT_MAX;
			int t;  // 临时值
			int minPivot;  // 分界点
			for (int k = i; k < j; k++) { // 寻找最优分界点
				t = m[i][k] + m[k + 1][j] + p[i] * p[k + 1] * p[j + 1];
				if (minCost > t) {
					minCost = t;
					minPivot = k;
				}
			}
			m[i][j] = minCost;
			s[i][j] = minPivot + 1;  // 加1是为了表示矩阵的编号
		}
	}
}
// 递归分治的方法找到最优结合方法
/*
	start: 开始矩阵
	end: 结束矩阵
	s: 动态规划中求得的分界点矩阵
*/
void traceBack(int start, int end, int s[][6]) {
	if (start == end) {
		printf("A%d", start);
		return;
	}
	int pivot = s[start - 1][end - 1];
	printf(" ( ");
	traceBack(start, pivot, s);  // 递归左边
	printf(" ) ");
	printf(" ( ");
	traceBack(pivot + 1, end, s); // 递归右边
	printf(" ) ");
}
void main() {
	int p[] = {30, 35, 15, 5, 10, 20, 25};  // 矩阵行列
	int n = sizeof(p) / sizeof(int) - 1;
	int m[n][n];
	memset(m, 0, sizeof(m));
	int s[n][n];
	memset(s, 0, sizeof(s));
	matrixChain(p, n, m, s);
	printf("m:                                        s:\n");
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			printf("%6d", m[i][j]);
		}
		printf("      ");
		for (int j = 0; j < n; j++) {
			printf("%6d", s[i][j]);
		}
		printf("\n");
	}
	printf("结合方法:\n");
	traceBack(1, n, s);
}

n; j++) {
printf(“%6d”, m[i][j]);
}
printf(" “);
for (int j = 0; j < n; j++) {
printf(”%6d", s[i][j]);
}
printf(“\n”);
}
printf(“结合方法:\n”);
traceBack(1, n, s);
}