算法矩阵列乘法 | 动态规划

Posted 2021-04-18 鸡翅编程

矩阵链乘法问题：给定 n 个矩阵的链<A1, A2, A3, ..., An>，矩阵 Ai 的规模为 。求完全括号化方案，使得计算乘积 所需要标量的乘法次数最少。

两个矩阵 A 和 B 只有相容，即 A 的列数等于 B 的行数时，才能相乘。如果 A 是 矩阵 ，B 是 矩阵，那么乘积 C 是 矩阵。计算 C 的进行了 次乘法运算。我们采用乘法进行的次数来表示计算的代价。

对于多个矩阵连乘的问题，我们称作矩阵链乘问题。矩阵的乘法满足结合律，因此任何加括号的方法并不会影响多个矩阵相乘的结果。 但不同的加括号的方案（改变计算顺序）会影响矩阵链乘法计算的代价：

---

一、算法剖析

1、最优括号化方案的结构特征

动态规划方法的第一步是寻找最优子结构，然后就可以利用这种子结构，从子问题的最优解构造出原问题的最优解。

（为方便起见，我们用 表示 的乘积结果矩阵）

• 从 整体往局部看：对于每一个矩阵链 ，如果对其括号化，我们应该在某两个矩阵 与 之间划分开，然后我们再讨论 和 的括号化。
• 再从 局部往整体看：计算 的代价 = 计算 的代价 + 计算 的代价 + 计算 的代价。

因此，为了构造一个矩阵链 乘法问题的最优解，我们可以将问题分解成两个子问题: 和 最优化括号问题。求出子问题的最优解后，再将子问题的最优解组合起来。

2、 一个递归的求解方案

• 矩阵的规模对应存储在一维数组
  p[0..n] 中，矩阵 的规模为 。
• 计算的代价对应储存在二维数组
  ans[1..n, 1..n] 中，我们用  ans[i, j] 表示矩阵 的最小计算代价，那么原问题的最小计算代价最终就在 ans[1, n] 内。
• 最优分割点对应储存在二维数组 divide[1..n, 1..n] 中，我们用 divide[i, j] 记录最优代价 ans[i, j] 对应的最优分割点 d。（此数组在打印最优结构时十分有用）

根据上述最优括号化方案的结构特征，我们可以列出状态转移方程：

• 当 i = j 时：
  
• 当 i < j 时：
  
• ps：
  在定义时将  ans 数组全部初始化为0。

---

二、算法实现

1、计算最优代价

从算法剖析中已知，本问题需要自底向上的方法：先处理好子矩阵链的最优问题，再总合起来。

那么我们实现时：应该从长度为2的矩阵链开始讨论，然后再在此基础上讨论长度为3的矩阵链...最终讨论长度的n的矩阵链，即得到答案。在某个已确定长度的讨论中，不可以漏去情况，不然会对最优值有影响。

那么我们的循环结构应该是这样：

• 第一层循环：长度
• 第二层循环：讨论每个长度为 l 的矩阵链 ：
• 第三层循环：确定最优分割点

/* 计算n元矩阵链p的最优代价
 * 动态规划的每一步结果储存在 ans与 divide中 */
void MatrixChainOrder(int p[], int n) {
    /* 矩阵链长度为2到n */
    for (int l = 2; l <= n; l++) {
        /* 讨论长度为l的矩阵链A[i..j] */
        for (int i = 1; i <= n - l + 1; i++) {
            int j = i + l - 1;
            ans[i][j] = INT_MAX;
            /* 依次讨论每一个分割点d:将矩阵链A[i..j]分成A[i..d]和 A[d+1..j] */
            for (int temp, d = i; d < j; d++) {
                temp = ans[i][d] + ans[d + 1][j] + p[i - 1] * p[d] * p[j];
                /* 记录下矩阵链A[i..j]最小的情况 */
                if (temp < ans[i][j]) {
                    ans[i][j] = temp;
                    divide[i][j] = d;
                }
            }
        }
    }
}

2、构造最优解

我们在二维数组 divide[i, j] 中已经记录了 的最优分界点，由于矩阵链的最优划分可以分解为子链的最优划分，故我们可以递归地打印出来。

/* 根据divide数组，输出A[i..j]的最优情况 */
void PrintDivide(int i, int j) {
    if (i == j)
        printf("A%d", i);
    else {
        int d = divide[i][j];  //找到分界点
        printf("(");
        PrintDivide(i, d);
        PrintDivide(d + 1, j);
        printf(")");
    }
}

---

三、例题

矩阵链乘问题

输入：
        共两行
第一行 N （ 1<=N<=100 ），代表矩阵个数。
第二行有 N+1 个数，分别为 A1 、 A2 ...... An+1 （ 1<=Ak<=2000 ）， Ak 和 Ak+1 代表第 k 个矩阵是个 Ak X Ak+1 形的。
        输出：
        共两行
第一行 M ，为最优代价。注：测试用例中 M 值保证小于
第二行为最优顺序。如 (A1((A2A3)A4)) ，最外层也加括号。注意:测试用例已经保证了输出结果唯一,所以没有AAA的情况.

---

具体方法上面已经详细讲解啦~

主要注意：当矩阵链中矩阵数目为1的输出要加上括号，不然会wa一个。

完整AC代码：

//
// Created by LittleCat on 2020/3/10.
//
 
#include <cstdio>
#include <climits>
 
using namespace std;
#define MAX 2010
 
int ans[MAX][MAX] = {0}; // 保存矩阵链 A[i..j]的最小代价
int divide[MAX][MAX] = {0};   // 记录最小代价ans[i,j]对应的分割点
 
/* 计算n元矩阵链p的最优代价
 * 动态规划的每一步结果储存在 ans与 divide中 */
void MatrixChainOrder(int p[], int n) {
    /* 矩阵链长度为2到n */
    for (int l = 2; l <= n; l++) {
        /* 讨论长度为l的矩阵链A[i..j] */
        for (int i = 1; i <= n - l + 1; i++) {
            int j = i + l - 1;
            ans[i][j] = INT_MAX;
            /* 依次讨论每一个分割点d:将矩阵链A[i..j]分成A[i..d]和 A[d+1..j] */
            for (int temp, d = i; d < j; d++) {
                temp = ans[i][d] + ans[d + 1][j] + p[i - 1] * p[d] * p[j];
                /* 记录下矩阵链A[i..j]最小的情况 */
                if (temp < ans[i][j]) {
                    ans[i][j] = temp;
                    divide[i][j] = d;
                }
            }
        }
    }
}
 
/* 根据divide数组，输出A[i..j]的最优情况 */
void PrintDivide(int i, int j) {
    if (i == j)
        printf("A%d", i);
    else {
        int d = divide[i][j];  //找到分界点
        printf("(");
        PrintDivide(i, d);
        PrintDivide(d + 1, j);
        printf(")");
    }
}
 
int main() {
    int n;
    scanf("%d", &n);
    int p[MAX];
    for (int temp, i = 0; i <= n; i++)
        scanf("%d", &p[i]);
 
    MatrixChainOrder(p, n);
    printf("%d
", ans[1][n]);
    if (n == 1)
        printf("(A1)");
    else
        PrintDivide(1, n);
    printf("
");
}
 

---

你AC了吗！