【动态规划入门】从矩阵链乘法到最优二叉查找树

闲话说在前头

进入ACM队以后，作为三人团队里最拉跨的一个，我默默地开启了天坑修补工作。最近在配合着《挑战程序设计竞赛2：算法与数据结构》刷AIZU的相关专题，个人感觉2比1通俗易懂多了，《挑战1》感觉啥都讲，但是啥都是点到为止，《挑战2》就逮着基础数据结构算法详细展开，感觉舒服多了。话说回来，最近推到了基础篇的动态规划专题，碰到了两个非常有意思的问题。记性差的我决定记录下来。

这两题在AIZU上的ID是 ALDS1_10_B和 ALDS_10_D。在这里我放了VJudge的链接。那我们就一个个来看吧。

先从简单的开始

（1）矩阵链乘法

矩阵链乘法问题的目标是找到最有效的方法来乘以给定的 $n$ 个矩阵 $M_1,M_2,M_3,...,M_n$。
编写一个程序，读取 $M_i$ 的维度，并找到最小次数标量乘法以计算 $M_1M_2...M_n$。

首先由线性代数的知识知道，矩阵乘法满足结合律，不满足交换律，一个 $p\times q$ 的矩阵和一个 $q \times r$ 的矩阵相乘需要 $p \times r \times q$ 次标量乘法运算，因为得到的是一个$p \times r$ 的矩阵，而这个矩阵每一个元素的计算都经历了$q$次乘法与$q-1$次加法。我们希望通过改变运算顺序，让总的乘法计算次数最小。自然，这种情况下加法的次数也会变小。

我们设这个矩阵链的第一个矩阵$M_1$大小为$p_0 \times p_1$， 第二个矩阵$M_2$大小$p_1\times p_2$，以此类推。

先从最简单的来看。很显然一个矩阵需要计算0次乘法。两个相邻矩阵$M_iM_{i+1}$只有一种计算方法，乘法次数为$p_{i-1} \times p_i \times p_{i+1}$。

那么三个矩阵$M_iM_{i+1}M_{i+2}$呢？由结合律，我们可以计算$(M_iM_{i+1})M_{i+2}$，也可以计算$M_i(M_{i+1}M_{i+2})$。前者需要计算$p_{i-1} \times p_i \times p_{i+1} + 0 + p_{i-1} \times p_{i+1} \times p_{i+2}$次，后者需要计算$0 + p_i \times p_{i+1} \times p_{i+2} + p_{i-1} \times p_{i} \times p_{i+2}$ 次。是不是出现了一些熟悉的身影？对了，三个矩阵组成的链可以拆分成前两个和后一个，或者前一个和后两个。而一个矩阵和两个矩阵的情况是确定的，因此只要两种情况取最小就行。

四个矩阵呢？就拆分成前一后三，前二后二，前三后一的情况来考虑，然后取最小值就行。一种、二种、三种的情况可以提前算出来，这就完成了已经计算的结果的再利用。接下来，一步步往上递推吧。当然这个递推式子在动态规划里面叫做“状态转移方程”。这个问题的状态转移方程用数学语言描述如下：

$dp[i][j]= \left \{ \begin{aligned} 0 ,i=j \\ p_{i-1} \times p_i \times p_{i+1} , j=i+1 \\ \min \{ dp[i][k]+dp[k+1][j]+p_{i-1} \times p_k \times p_j \enspace | \enspace k=i , i+1 , \dots , j-1 \} , others \end{aligned} \right.$

我们设$dp[i][j]$表示第$i$个矩阵到第$j$个矩阵相乘的最小乘法次数。最终我们要通过递推关系计算得$dp[1][n]$的值。根据这道题的样例输入，我们可以这样读入矩阵的尺寸：

#define rep(i,a,b) for(int i=a;i<=b;++i)
#define IN2(a,b) scanf("%d %d",&a,&b)
rep(i, 1, n) {
	IN2(p[i - 1], p[i]);
	dp[i][i] = 0;
}

并且将$dp[i][i]$赋值为0。

接下来，分别计算矩阵链里连续2个、3个、……、 n个矩阵相乘的最小标量乘法次数。设长度为len，则若是第 $i$ 个矩阵开始的连续len个矩阵，那么最后一个矩阵就是第$i+len-1$个。我们就可以写出以下关键代码：

#define min(a,b) ((a)<(b)?(a):(b))
rep(len, 2, n) {
	rep(i, 1, n - len + 1) {
		int j = i + len - 1;
		dp[i][j] = INT_MAX;
		rep(k, i, j - 1) {
			dp[i][j] = min(dp[i][j], dp[i][k] + dp[k + 1][j] + p[i - 1] * p[k] * p[j]);
		}
	}
}

接下来上完整代码：（学校的OJ系统这个专题只允许交C语言代码，这里就放C版本的了）

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <stdbool.h>
#include <limits.h>
#define IN1(a) scanf("%d",&a)
#define IN2(a,b) scanf("%d %d",&a,&b)
#define IN3(a,b,c) scanf("%d %d %d",&a,&b,&c)
#define GRP int T;scanf("%d",&T);rep(C,1,T)
#define etr putchar('\n')
#define rep(i,a,b) for(int i=a;i<=b;++i)
#define rrep(i,a,b) for(int i=a;i>=b;--i)
#define mem(arr,val) memset(arr,val,sizeof(arr))
#define min(a,b) ((a)<(b)?(a):(b))
typedef long long ll;
typedef unsigned long long ull;

int n;
int p[110], dp[110][110];

int main(void) {
	IN1(n);
	rep(i, 1, n) {
		IN2(p[i - 1], p[i]);
		dp[i][i] = 0;
	}
	rep(len, 2, n) {
		rep(i, 1, n - len + 1) {
			int j = i + len - 1;
			dp[i][j] = INT_MAX;
			rep(k, i, j - 1) {
				dp[i][j] = min(dp[i][j], dp[i][k] + dp[k + 1][j] + p[i - 1] * p[k] * p[j]);
			}
		}
	}
	printf("%d\n", dp[1][n]);
	return 0;
}

提高了一些难度

（2）最优二叉查找树

给定n个互异的关键字组成的序列$K=<k_1,k_2,...,k_n>$，且关键字有序$(k_1<k_2<...<k_n)$，我们想用这些关键字构造一棵二叉查找树。对每个关键字$k_i$，一次搜索搜索到的概率为$p_i$。可能有一些搜索的值不在$K$内，因此还有$n+1$个“失败节点”$d_0,d_1,...,d_n$作为虚拟键，他们代表不在$K$内的值。具体：$d_0$代表所有小于$k_1$的值，$d_n$代表所有大于$k_n$的值。而对于$i = 1,2,...,n-1$,虚拟键$d_i$代表所有位于$k_i$和$k_{i+1}$之间的值。对于每个虚拟键，一次搜索对应于$d_i$的概率为$q_i$。要使得查找一个节点的代价（代价定义为：从根节点到目标节点的路径上的节点数目）期望值最小，就需要建立一棵最优二叉查找树。

在这里我们首先进行反推：如果一棵树是最优二叉查找树，那么它的两个子树（如果有的话）也一定是最优二叉查找树。否则子树还可以再继续优化，使得局部期望值减小，进而总期望值减小，不符合假设条件。

由题意，失败节点一定是叶节点。我们按照查找树的层数对这个问题进行考虑，$n$个关键字组成的查找树最高可能有$n$层。当处于最底层的时候，子树中只有一个节点，按照题意一定是失败节点，期望代价就是$q_i$。

我们往上爬一层。二叉查找树中序遍历得到的序列一定有序，因此如果根节点是$k_i$，那么它的左右儿子一定是$d_{i-1}$和$d_i$。对子节点而言，深度增加了1，因此代价为$2 \times ( q_{i-1} + q_i ) + p_i$。

接下来爬第三层。假设根节点$k_i$，那么可能的情况有两种，需要分别计算代价并取最小值：

• 左子树$d_{i-1}$，右子树为$k_{i+1} , d_i , d_{i+1}$组成的子树。此时的代价为$2 \times q_{i-1} + [ 2 \times p_{i+1} + 3 \times ( q_i + q_{i+1} ) ] +p_{i}$。
• 左子树为$k_{i-1},d_{i-2},d_{i-1}$组成的子树，右子树为$d_i$。此时的代价为$[ 2 \times p_{i-1} + 3 \times ( q_{i-2} + q_{i-1} ) ] + 2 \times q_{i} +p_{i}$。

我们对第一种情况得到的式子做一个变形，得到$(q_{i-1} + p_{i-1} + q_i + q_{i+1} ) + \{1 \times q_{i-1} + [ 1 \times p_{i+1} + 2 \times ( q_i + q_{i+1} ) ] \} +p_{i}$。我们发现，第一个括号里面的是根节点所有子树取到的概率（权值）的总和，第二个括号内是左右子树的平均查找代价，最后再加上根节点的概率。其实这也很好理解：往上把两颗树挂在一个节点上的时候，子树的深度都加了1，所以需要把子树的概率乘以增加的深度再进行求和，这就是增加的查找代价。对于第二种情况，我们也可以做同样的分析。看到这个挂节点的过程，是不是和矩阵链乘法的推理过程非常像？

接下来，我们仿照矩阵链乘法选取一个区间进行分析，定义几个变量。令$w[i][j]$表示包含关键字 $k_i,k_{i+1},\dots ,k_{j}$的二叉查找树的概率总和，$e[i][j]$表示这段关键字构成的最优二叉查找树的代价期望。最终，我们想要得到$e[1][n]$。若在其中选择$k_r$作为树的根，则我们有

$$e[i][j]=p_r+(e[i][r-1]+w[i][r-1])+(e[r+1][j]+w[r+1][j])$$

注意到$w[i][r-1]+p_r+w[r+1][j]=w[i][j]$，因此我们可以得到下面的关系：

$$e[i][j]=e[i][r-1]+e[r+1][j]+w[i][j]$$

可以验证，这个关系对第二层也是成立的。因此我们可以得到这样的状态转移方程：

$$e[i][j]= \left \{ \begin{aligned} q_{i-1} \enspace , j=i-1 \\ \min \{ e[i][r-1]+e[r+1][j]+w[i][j] \enspace | \enspace r=i , i+1 , \dots , j \} , others \end{aligned} \right.$$

我们发现，$w[i][j]$也可以通过动态规划的思想进行计算。首先我们初始化数组：

rep(i, 1, n + 1) {
	w[i][i - 1] = q[i - 1];
	e[i][i - 1] = q[i - 1];
}

接下来DP过程中的关键代码和矩阵链乘法的关键代码一模一样，这里不再赘述。这里我们用到了式子$w[i][j] = w[i][j - 1] + p_j + q_j$。

rep(len, 1, n) {
	rep(i, 1, n - len + 1) {
		int j = i + len - 1;
		e[i][j] = INF;
		w[i][j] = w[i][j - 1] + p[j] + q[j];
		rep(k, i, j) {
			e[i][j] = min(e[i][j], e[i][k - 1] + e[k + 1][j] + w[i][j]);
		}
	}
}

接下来上完整代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <stdbool.h>
#define IN1(a) scanf("%d",&a)
#define IN2(a,b) scanf("%d %d",&a,&b)
#define IN3(a,b,c) scanf("%d %d %d",&a,&b,&c)
#define GRP int T;scanf("%d",&T);rep(C,1,T)
#define etr putchar('\n')
#define rep(i,a,b) for(int i=a;i<=b;++i)
#define rrep(i,a,b) for(int i=a;i>=b;--i)
#define elif else if
#define mem(arr,val) memset(arr,val,sizeof(arr))
#define min(a,b) ((a)<(b)?(a):(b))
typedef long long ll;
typedef unsigned long long ull;

const double INF = 1e9;
double w[510][510], e[510][510];
double p[510], q[510];
int n;

int main(void) {

	IN1(n);
	rep(i, 1, n) {
		scanf("%lf", &p[i]);
	}
	rep(i, 0, n) {
		scanf("%lf", &q[i]);
	}

	rep(i, 1, n + 1) {
		w[i][i - 1] = q[i - 1];
		e[i][i - 1] = q[i - 1];
	}
	rep(len, 1, n) {
		rep(i, 1, n - len + 1) {
			int j = i + len - 1;
			e[i][j] = INF;
			w[i][j] = w[i][j - 1] + p[j] + q[j];
			rep(k, i, j) {
				e[i][j] = min(e[i][j], e[i][k - 1] + e[k + 1][j] + w[i][j]);
			}
		}
	}
	printf("%.6f", e[1][n]);

	return 0;
}

感谢你看到这里~