算法-动态规划

完整可执行代码见 sun2ot/algorithm-learning

基本思想

与分治法类似，其基本思想也是将待求解问题分解成若干个子问题，先求解子问题，然后从这些子问题的解得到原问题的解。

与分治法不同的是，适合于用动态规划法求解的问题，经分解得到的子问题往往不是互相独立的，有些子问题被重复计算了许多次。如果我们能够保存已解决的子问题的答案，而在需要时再找出已求得的答案，这样就可以避免大量的重复计算，从而得到多项式时间算法。

为了达到此目的，可以用一个表来记录所有已解决的子问题的答案。不管该子问题以后是否被用到，只要它被计算过，就将其结果填入表中。这就是动态规划法的基本思想。

矩阵连乘问题

问题定义

矩阵连乘算法问题是一个经典的优化问题，其目标是找到一种最优的计算顺序，以最小化矩阵连乘操作的总次数。

假设有一系列矩阵 $A_1, A_2, A_3, ..., A_n$ ，其中矩阵 $A_i$ 的维度为 $p_{i-1} \times p_i$ ，我们想要计算它们的连乘积 $A_1 \cdot A_2 \cdot A_3 \cdot ... \cdot A_n$ 。注意，矩阵连乘是一个不可交换的操作，即矩阵的乘法顺序会影响最终结果。

矩阵连乘问题的目标是找到一种最优的计算顺序，使得计算矩阵连乘的总次数最小。这个问题可以使用动态规划的方法来解决，其中关键是定义一个递归的状态和状态转移方程。

分析最优子结构

将矩阵连乘积 $A_iA_{i+1}...A_j$ 记为 $A[i:j]$ 。假设最优计算次序是在 $A_k$ 和 $A_{k+1}$ 之间断开，则原问题的解即为 $A[1:k]$ 的计算量加上 $A[k+1:n]$ 的计算量再加上 $A[1: k]\cdot A[k+1:n]$ 的计算量。

那么这个拆解过程存在一个前提，就是计算 $A[1:n]$ 的最优次序所包含的矩阵子链 $A[1:k]$ 和 $A[k+1:n]$ 的次序也是最优的。我们可以用反证法来证明这一点是成立的。

假设当 $A[1:k], A[k+1:n]$ 这种计算次序最优时， $A[1:k]$ 的计算次序不是最优的，即存在 $k'$ ，使得 $A[1:k']$ 的计算量小于 $A[1:k]$ 。那么显然， $A[1: k'], A[k'+1:n]$ 这种拆解方式的计算量一定小于原先的 $A[1:k]$ 和 $A[k+1:n]$ ，即 $A[1:k], A[k+1:n]$ 的计算次序不是最优，这与假设矛盾。

同理可证 $A[1: n]$ 所包含的矩阵子链 $A[k+1:n]$ 也是最优的。

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

问题的最优子结构性质是该问题可用动态规划算法求解的显著特征。

算法思想

在确定了原问题具有最优子结构性质后，我们需要递归地定义最优值。假设计算 $A[i:j]$ 所需的最少数乘次数为 $m[i][j]$ ，则原问题 $A[1:n]$ 的最优值为 $m[1][n]$ 。

• 当 $i=j$ 时， $A[i:j]$ 就是一个单独的矩阵，所以 $m[i][j]=0$
• 当 $i \lt j$ 时， $m[i][j]=m[i][k]+m[k+1][j]+p_{i-1}p_kp_j$

$i<j$ 时的解释
关于 $p$ ，我们将矩阵 $A_i$ 的维数记为 $p_{i-1} \times pi$ 。此处有一个问题，就是对于两个矩阵相乘，总共的数乘次数为多少？若 $A$ 是 $p \times q$ 矩阵， $B$ 是 $q\times r$ 矩阵，则计算 $AB$ 总共需要 $pqr$ 次数乘。

p 行，每行都要乘 r 列，一共就是 pr 次乘法。每一次乘法都是 q 个数和 q 个数相乘，因此总共为 prq=pqr 次。

因此， $m[i][k]$ 表示从 $A_i$ 乘到 $A_k$ ，其中 $A_i$ 是 $(i-1) \times i$ 矩阵， $A_k$ 是 $(k-1) \times k$ 矩阵，由于矩阵连乘的前提一定是矩阵链是可乘的，因此 $A[i:k]$ 的结果是一个 $(i-1) \times k$ 矩阵。同理， $A[k+1:j]$ 的结果是 $k \times j$ 矩阵。所以这两个矩阵相乘的数乘次数为 $p_{i-1}p_kp_j$ 。

那么 $k$ 取多少呢？我们无法确定，因为 $i \leq k \lt j$ 的位置有 $i,i+1,...,j-1$ 共 $j-i$ 种可能，因此 $m[i][j]$ 递归地定义为：

$$m[i][j] = \begin{cases} 0 & i=j \\ \min\limits_{i \leq k \lt j}\{m[i][j]=m[i][k]+m[k+1][j]+p_{i-1}p_kp_j\} & i \lt j \end{cases}$$

上式给出了两个信息，一是最优值 $m$ ，即矩阵连乘所需的最小次数；二是最优解 $k$ ，即应该怎么划分才能获得最优值。

代码实现

/**
 * 使用动态规划解决矩阵连乘问题，并返回最小计算次数和最优解
 * 
 * @param dimensions 矩阵的维度数组
 * @param n 矩阵的数量
 * @param dp 存储最小计算次数的二维数组
 * @param split 存储最优划分点的二维数组
 * @return 最小计算次数
 */
int matrixChainOrder(int dimensions[], int n, int dp[MAX_MATRICES][MAX_MATRICES], int split[MAX_MATRICES][MAX_MATRICES]) {
    // 初始化dp数组和split数组
    for (int i = 1; i <= n; i++) {
        dp[i][i] = 0;
    }

    // 计算子问题的最优解和划分点
    for (int len = 2; len <= n; len++) {  // 子问题规模
        for (int i = 1; i <= n - len + 1; i++) {  // 子问题的起始位置
            int j = i + len - 1;  // 子问题的结束位置
            dp[i][j] = INT_MAX;  // 初始化最优值为最大值

            for (int k = i; k < j; k++) {  // 子问题的划分点
                int cost = dp[i][k] + dp[k + 1][j] + dimensions[i - 1] * dimensions[k] * dimensions[j];
                if (cost < dp[i][j]) {  
                    dp[i][j] = cost;  // 更新最优值
                    split[i][j] = k; // 记录最优划分点
                }
            }
        }
    }

    // 返回最小计算次数
    return dp[1][n];
}


/**
 * 输出最优解
 * 
 * @param split 分割矩阵的位置
 * @param i 矩阵的起始位置
 * @param j 矩阵的结束位置
 */
void printOptimalParenthesis(int split[MAX_MATRICES][MAX_MATRICES], int i, int j) {
    if (i == j) {  // 矩阵只有一个时，直接输出
        printf("A%d", i);
    } else {  // 矩阵有多个时，递归输出
        printf("(");  // 输出左括号
        printOptimalParenthesis(split, i, split[i][j]);  // 输出左边的矩阵链的
        printOptimalParenthesis(split, split[i][j] + 1, j);  // 输出右边的矩阵链
        printf(")");  // 输出右括号
    }
}

在这段代码中，三重循环的循环变量分别代表以下内容：

1. 外层循环变量 len：表示子问题的规模，从2开始逐渐增加，直到达到矩阵的大小 n。这个循环用于控制子问题的规模。
2. 中间循环变量 i：表示子问题的起始位置，从1开始逐渐增加，直到 n - len + 1。这个循环用于遍历所有可能的子问题起始位置。
3. 内层循环变量 k：表示子问题的划分点，从 i 开始逐渐增加，直到 j。这个循环用于遍历所有可能的划分点，计算出最优的划分点。

关于 len

循环从 len = 2 开始，是因为在矩阵乘法中，我们需要先计算较小的子矩阵的乘积，然后逐渐增加子矩阵的大小。通过从较小的子矩阵开始，我们可以利用已经计算出的乘积结果来计算更大的子矩阵的乘积。

例如，对于一个 4x4 的矩阵，我们首先计算 2x2 的子矩阵的乘积，然后计算 3x3 的子矩阵的乘积，最后计算整个 4x4 矩阵的乘积。通过从较小的子矩阵开始，我们可以利用已经计算出的 2x2 子矩阵的乘积来计算 3x3 子矩阵的乘积，从而减少计算量。

因此，循环从 len = 2 开始，以便从较小的子矩阵开始计算乘积，并逐渐增加子矩阵的大小，直到计算整个矩阵的乘积。

为什么 $i \leq n-len+1$ ？

假设矩阵下标从 1 开始到 n。当问题规模为 len 时，即 $[1:len],\ [2:len+1], ...$ 这样逐个划分计算，其对应的起始位置就是 $1,2,...$ 。那么，我们要求起始位置的上界，对应的极端情况显然是长度为 len 的矩阵链在整个矩阵连乘的式子中尽可能地靠右，这样得到的起始位置就是尽可能大的。

所以问题就变成了，已知最后（最右的）一个矩阵为 n，从它开始向左推一个长度为 len 的矩阵链，这个链的起始位置是谁? 显然，答案是 $n-len+1$ 。

矩阵的维度有两个值，为什么可以通过一维数组 dimensions 来表示？

参见上文两矩阵的数乘次数问题，我们可知，对于第 i 个矩阵， dimensions[i-1] 是它的行数， dimensions[i] 是它的列数。这是因为，如果我们按照顺序进行矩阵链乘法，那么第 i 个矩阵的列数必须等于第 i+1 个矩阵的行数（矩阵可乘，左列等于右行），所以我们可以用一个连续的数组来存储所有矩阵的行数和列数。

最长公共子序列

问题定义

给定两个序列 $X=\{x_1, x_2,..., x_m\},\, Y=\{y_1, y_2, ..., y_n\}$ ，找出最长的、公共的子序列。

分析最优子结构

最长公共子序列问题有没有最优子结构性质，说白了就是，如果有最长公共子序列 $Z_k=\{z_1,z_2,...,z_k\}$ ， $Z_{k-1}$ 是不是第二长的？（当然是的）下面给出该问题的最优子结构性质定义：

设序列 $X_m=\{x_1, x_2,..., x_m\},\, Y_n=\{y_1, y_2, ..., y_n\}$ 的最长公共子序列为 $Z_k=\{z_1,z_2,...,z_k\}$ ，则

• 若 $x_m=y_n$ ，则 $z_k=x_m=y_n$ ，且 $Z_{k-1}$ 是 $X_{m-1}$ 和 $Y_{n-1}$ 的最长公共子序列
• 若 $x_m\neq y_n$ 且 $z_k\neq x_m$ ，则 $Z_k$ 是 $X_{m-1}$ 和 $Y$ 的最长公共子序列
• 若 $x_m\neq y_n$ 且 $z_k\neq y_n$ ，则 $Z_k$ 是 $X$ 和 $Y_{n-1}$ 的最长公共子序列

下面分别针对这三点给出证明

若 $x_m=y_n$ 时，$z_k\neq x_m$，那意味着最长公共子序列 $Z_k$ 不是最长的。因为 $\{z_1,z_2,...z_k,x_m\}$ 显然也是公共子序列，且长度为 $k+1$，比 $Z_k$ 还长，构成矛盾。所以当 $x_m=y_n$ 时，必有 $z_k=x_m=y_n$，同时 $Z_{k-1}$ 是 $X_{m-1}$ 和 $Y_{n-1}$ 的公共子序列。

接下来证明 $Z_{k-1}$ 是最长的。假设 $X_{m-1}$ 和 $Y_{n-1}$ 有长度 $\gt k-1$ 的公共子序列 $W$，那就八达鸟了。因为 $W$ 的长度至少是 $k$，再把两个父串共有的 $x_m$ 加上，$W$ 的长度直接达到 $k+1$，此之谓无中生有，故矛盾。

---

若 $x_m\neq y_n$ 且 $z_k\neq x_m$ ，则 $Z_k$ 肯定是 $X_{m-1}$ 和 $Y$ 的公共子序列，不过是不是最长的呢？还是反证法，我们假设 $X_{m-1}$ 和 $Y$ 有长度大于 $k$ 的公共子序列 $W$，显然 $W$ 也是 $X$ 和 $Y$ 的公共子序列。但 $X$ 和 $Y$ 的最长公共子序列为 $Z_k=\{z_1,z_2,...,z_k\}$ 长度也只有 $k$，故矛盾。

---

第三点证明与二同理，略。

综上，最长公共子序列问题具有最优子结构性质。

算法思想

根据最优子结构性质可知，

• 任一串为空时，最长公共子序列长度为 0
• 当 $x_m=y_n$ 时，找出 $X_{m-1}$ 和 $Y_{n-1}$ 的最长公共子序列，然后尾部加上 $x_m\ \text{or}\ y_n$ 即为 $X$ 和 $Y$ 的最长公共子序列。
• 当 $x_m\neq y_n$ 时，取 $X_{m-1}$ 和 $Y$、$X$ 和 $Y_{n-1}$ 的最长公共子序列中更长的一个

将序列 $X_i$ 和 $Y_j$ 的最长公共子序列长度记为 $c[i][j]$，根据以上内容，构建递推关系式如下：

$$c[i][j] = \begin{cases} 0 & i=0\ \text{or}\ j=0 \\ c[i-1][j-1] + 1 & i,j \gt 0;\ x_i = y_j \\ \max\{c[i][j-1],c[i-1][j]\} & i,j \gt 0;\ x_i \neq y_j \end{cases}$$

代码实现

/**
 * 使用动态规划解决最长公共子序列问题
 * 
 * @param str1 字符串1
 * @param str2 字符串2
 */
void longestCommonSubsequence(char* str1, char* str2) {
    int m = strlen(str1);
    int n = strlen(str2);

    int dp[m + 1][n + 1];

    // 初始化dp数组的第一行和第一列
    for (int i = 0; i <= m; i++) {  
        dp[i][0] = 0;  
    }
    for (int j = 0; j <= n; j++) {
        dp[0][j] = 0;
    }

    // 计算dp数组的其他元素
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= n; j++) {
            if (str1[i - 1] == str2[j - 1]) {  // 两个子串的最后一个字符相同
                dp[i][j] = dp[i - 1][j - 1] + 1;
            } else {
                dp[i][j] = (dp[i - 1][j] > dp[i][j - 1]) ? dp[i - 1][j] : dp[i][j - 1];
            }
        }
    }

    // 构造最长公共子序列
    int lcsLength = dp[m][n];
    char lcs[lcsLength + 1];
    lcs[lcsLength] = '\0';

    int i = m, j = n;
    while (i > 0 && j > 0) {
        // 串的扫描和LCS的构建都是从右往左，倒着来
        if (str1[i - 1] == str2[j - 1]) {  // 两个串的尾部字符相等，说明该字符属于LCS
            lcs[lcsLength - 1] = str1[i - 1];
            i--;
            j--;
            lcsLength--;
        } else if (dp[i - 1][j] > dp[i][j - 1]) {  // 删去str1最后一个字符得到的LCS比删去str2最后一个字符得到的LCS长，因此i向前回溯
            i--;
        } else {
            j--;
        }
    }
    printf("LAS Length is: %d\n", dp[m][n]);
    printf("LCS is: %s\n", lcs);  // 注意，最长公共子序列并不唯一
}