【动态规划入门】从矩阵链乘法到最优二叉查找树
闲话说在前头
进入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\)作为树的根，则我们有
注意到\(w[i][r-1]+p_r+w[r+1][j]=w[i][j]\)，因此我们可以得到下面的关系：
可以验证，这个关系对第二层也是成立的。因此我们可以得到这样的状态转移方程：
我们发现，\(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;
}
感谢你看到这里~
