C语言矩阵连乘 (动态规划)详解

算法原理

矩阵乘法不满足交换律和结合律，因此计算矩阵连乘的顺序会影响计算时间。即使只有6个矩阵相乘，也有可能有超过百万种计算次序。因此需要通过算法来优化时间复杂度。动态规划是一种可用于求解最优化问题的算法，它将原问题分解为子问题求解，并将每个子问题的最优解存储在表格中，以便在较大的子问题中简化计算。

设矩阵 $A_{1 \times 5}$、$A_{5 \times 10}$、$A_{10 \times 3}$ 和 $A_{3 \times 12}$ 依次相乘，分别用 $A_1$、$A_2$、$A_3$ 和 $A_4$ 表示。根据动态规划的思想，将矩阵连乘问题分解为子问题，即：

$M(1,4)=min{M(1,k)+M(k+1,4)+p_{0}p_{1}p_{4}}$

其中 $p_0$、$p_1$ 和 $p_4$ 分别表示 $A_1$ 的列数、$A_2$ 的列数和 $A_4$ 的行数，$M(1,4)$ 表示将矩阵 $A_1$ 到 $A_4$ 连乘所需要的最少次数，$k$ 表示从 $A_1$ 到 $A_4$ 中的任意一矩阵，即矩阵乘法的括号位置。

算法实现

步骤一：分配空间

由于要计算从 $A_1$ 到 $A_n$ 的所有矩阵连乘次序的最小值，因此需要分配一个 $n \times n$ 的 int 类型数组存储最小次数。

int m[MAXSIZE][MAXSIZE];

同样地，还需要分配一个 $n \times n$ 的 int 类型数组存储最优解对应的括号位置。

int s[MAXSIZE][MAXSIZE];

步骤二：计算最小次数和最优解括号位置

根据矩阵连乘算法原理，可以用递归方法暴力枚举所有可能次序，并计算出最小次数和最优解括号位置，时间复杂度为$O(2^n)$，其中 $n$ 表示矩阵个数。

为了提高速度，使用动态规划算法。将矩阵乘法的次序由2个矩阵相乘逐渐扩大，直到求解所有 $n$ 个矩阵相乘。假设现在要求解从 $A_i$ 到 $A_j$ 连乘所需的最小次数，设 $k$ 为括号位置，则可以得到递推式：

$$ M(i,j) = \min{M(i,k) + M(k+1,j) + p_{i-1}p_kp_j} $$

其中 $p_i$ 表示第$i$个矩阵的行数和第$i+1$个矩阵的列数，$M(i,j)$ 表示将矩阵 $A_i$ 到 $A_j$ 连乘所需要的最少次数。为了实现上表达式中的$M(i,k)$和$M(k+1,j)$，需要额外存储一个 $n \times n$ 的 int 类型数组存储中间过程。

我们只需要枚举括号的位置 $k$，得到最小的连乘次数 $M(i,j)$。同时，将括号位置 $k$ 存储到 $s(i,j)$ 中。动态规划问题的最终结果为 $M(1,n)$，此时需要用 $s$ 数组中存储的括号位置信息来生成最优解。

void matrixChainOrder() {
    for (int i = 1; i <= n; i++) {
        m[i][i] = 0;
    }
    for (int l = 2; l <= n; l++) {
        for (int i = 1; i <= n - l + 1; i++) {
            int j = i + l - 1;
            m[i][j] = INT_MAX;
            for (int k = i; k < j; k++) {
                int q = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];
                if (q < m[i][j]) {
                    m[i][j] = q;
                    s[i][j] = k;
                }
            }
        }
    }
}

步骤三：输出最优解

根据存储的括号位置信息 $s$，可以生成最优解的字符串，即括号序列。下面给出一个示例输出。

假设有矩阵 $A_1$、$A_2$、$A_3$ 和 $A_4$，形状分别为 $1 \times 5$、$5 \times 10$、$10 \times 3$ 和 $3 \times 12$。将这四个矩阵连乘所需要的最少次数为 620。最优解为 $((A_1(A_2A_3))A_4)$。

代码如下：

void printOptimalParens(int i, int j) {
    if (i == j) {
        printf("A%d", i);
    } else {
        printf("(");
        printOptimalParens(i, s[i][j]);
        printOptimalParens(s[i][j] + 1, j);
        printf(")");
    }
}

示例说明

假设有四个矩阵相乘，形状分别为 $A_{1 \times 5}$、$A_{5 \times 10}$、$A_{10 \times 3}$ 和 $A_{3 \times 12}$。代码实现如下：

#include <stdio.h>
#include <limits.h>

#define MAXSIZE 100

int n = 4;
int p[MAXSIZE + 1] = { 1, 5, 10, 3, 12 };
int m[MAXSIZE][MAXSIZE];
int s[MAXSIZE][MAXSIZE];

void matrixChainOrder() {
    for (int i = 1; i <= n; i++) {
        m[i][i] = 0;
    }
    for (int l = 2; l <= n; l++) {
        for (int i = 1; i <= n - l + 1; i++) {
            int j = i + l - 1;
            m[i][j] = INT_MAX;
            for (int k = i; k < j; k++) {
                int q = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];
                if (q < m[i][j]) {
                    m[i][j] = q;
                    s[i][j] = k;
                }
            }
        }
    }
}

void printOptimalParens(int i, int j) {
    if (i == j) {
        printf("A%d", i);
    } else {
        printf("(");
        printOptimalParens(i, s[i][j]);
        printOptimalParens(s[i][j] + 1, j);
        printf(")");
    }
}

int main() {
    matrixChainOrder();
    printf("最少次数：%d\n", m[1][n]);
    printf("最优解：%s\n", "");
    printOptimalParens(1, n);
    printf("\n");
    return 0;
}

运行结果如下：

最少次数：620
最优解：((A1(A2A3))A4)

在这个示例中，我们将 $A_1$ 到 $A_4$ 连乘，记为 $((A_1A_2)A_3)A_4$ 和 $A_1(A_2(A_3A_4))$，其中前者需要计算 $90+500+1500=2090$ 次，后者需要计算 $5+180+4320=4510$ 次。由此可见，采用动态规划算法可以得到最优解，而暴力枚举算法时间复杂度太高，无法得到最优解。