算法：矩阵连乘求的最小乘法次数

2018-12-29 本文已影响24人柳年思水

今天来讨论一道算法题，这道算法题我在做的时候真的是没什么思路，甚至看到解法之后依然想了好久才想明白，好久没看过算法了，本来对动态规划这块就不是很熟，这里记录一下。

题目

给定一个 n 的矩阵序列，我们希望计算它们的乘积：

$A_1·A_2·A_3····A_n$

其中， $A_i$ 是 $a_i \times a_{i+1}$ 矩阵

注意，这里不是要计算乘积，而是希望找到一个明确的计算次序，使得这个矩阵连乘的乘法次数最少，并求这个最小的乘法次数（ $m(1, n)$ ，这个值表示第 1 个到第 n 个矩阵相乘的最小乘法次数）。下面举举几个例子：

当 $n=1$ 时， $m(1,1) = 0$ ；
当 $n=2$ 时， $m(1,2) = a_1 * a_2 * a_3$ ；
当 $n=3$ 时，有两种情况：

$((A_1 · A_2) · A_3)$ ，这乘法次数为： $a_1*a_2*a_3 + a_1*a_3*a_4$ ；
$(A_1 · (A_2 · A_3))$ ，这乘法次数为： $a_2*a_3*a_4 + a_1*a_2*a_4$ ；

而， $m(1,3) = min \begin{Bmatrix} a_1*a_2*a_3 + a_1*a_3*a_4, a_2*a_3*a_4 + a_1*a_2*a_4 \end{Bmatrix}$

简单来说，这道题的目的，就是在计算矩阵连乘时，求出一种方案，使得计算所需的乘法次数最小，输出的结果是这个最小的乘法次数。

暴力穷尽解法

能想到的第一种方案就是，穷举法，这里不讲述穷举法如何解决上面的问题，主要来看下穷举法的复杂度，假设 $P(n)$ 表示对于一个 n 个矩阵连乘时可选的方案数（穷举法，需要计算每个方案最后的结果，然后选择最小的一个），那么其实是可以得到下面的递推公式的：

$P(n)= \begin{cases} 1 &\mbox{n=1}\\\ \sum_{k=1}^{n-1} P(k)·P(n-k) &\mbox{$n \geq 2$} \end{cases}$

对于一个 n 个矩阵连乘情况，我们在矩阵中选择一个 k 点，将一个大矩阵分成两个小矩阵，然后分别求其可选的方案数，结果相乘，就是此时的总方案数。那么这种暴力穷尽解法的时间复杂度是多少呢？

这个时间复杂度其实是： $\Omega(2^n)$ ，我们来证明一下：

假设 $P(k) \geq c·2^k , c \leq \frac{1}{2}$

当 $n=1$ 时， $P(1) = 1 \geq c*2^1$ ；
当 $n\geq1$ 时， $P(n) = \sum_{k=1}^{n-1} { P(k)P(n-k) } \geq c^2·2^1·2^{n-1} + c^2·2^2·2^{n-2} + ... + c^2·2^{n-1}·2^1 \geq c^2·(n-1)·2^n$ ；

所以，当 $c(n-1) \geq 1$ ，即 $n \geq \frac{1}{c} +1$ 时， $P(n) \geq c·2^n$ ；

因此， $P(n) = \Omega(2^n)$ ，这种暴力穷尽解法的时间复杂度是 $\Omega(2^n)$ 。

时间复杂度是与 n 呈指数关系，这种方案明显是一个糟糕的选择。

动态规划解法

前面的穷举法时间度过高，明显是不可取，那么应该怎么去优化这个问题呢？这里介绍一种使用动态规划的解法。

子问题拆分

其实，在前面的穷尽法中有一个可取的地方，那就是子问题的拆分，在 $\begin{Bmatrix} A_1 ... A_n\end{Bmatrix}$ 中，先选择一个矩阵 $A_k$ ，这样的话就将一个大矩阵拆分为两个小矩阵，分别求这两个小矩阵的最小乘法次数，然后再将 i 从 $0...n$ 遍历一边，取一个最小值，就可以得到我们想要的结果。

最开始考虑这个问题时，当时没想明白，一个大矩阵的最优解怎么拆分为两个小矩阵的最优解，忘记了最外围还有一个遍历，然后再取最优解。

这是一种非常有用算法思想，将一个大问题拆分为一些子（小）问题，先解决这些小问题，最后大问题就迎刃而解了，其实跟递归、分治这些思想都有一些类似之处。

递推公式推导

在 $\begin{Bmatrix} A_1 ... A_n\end{Bmatrix}$ 中，选择一个矩阵 $A_k ( 1 \leq k < n)$ 将一个大矩阵拆分为 $\begin{Bmatrix} A_1 ... A_k\end{Bmatrix}$ 和 $\begin{Bmatrix} A_{k+1} ... A_n\end{Bmatrix}$ 的问题，最后把这两个做乘就可以得到最后的结果，，因此，可以得到（k 为最佳切分点时）：

$m(1,n) = m(1, k)+m(k+1, n) + a_1*a_{k+1}*a_n$

上面的公式是假设 $A_k$ 为最佳的切分点，但是实际上这个是无法判断的，对于
$\begin{Bmatrix} A_1 ... A_n\end{Bmatrix}$ ，k 是有 $n-1$ 种可能的值，最优的切分点必然在这其中，只需要遍历求最优解即可，最后的递推公式为：

$m(i,j)= \begin{cases} 0 &\mbox{i=j}\\\ min_{i \leq k < j} \begin{Bmatrix} m(i, k)+m(k+1, j) + a_i*a_{k+1}*a_{j+1} \end{Bmatrix} &\mbox{$ i < j$} \end{cases}$

到这里，最优解的递推公式就推导完成了，那么如何求解呢？

计算最小的乘法次数

有了前面的递推公式，下面就看如何根据这个递推去求最优解？

递归算法

如果直接按照递推公式去设计程序，其实现如下：

public int countMatricsChain(int[] a, int i, int j, int[][] m) {
    if (i == j) {
        return 0;
    }
    m[i][j] = Integer.MAX_VALUE;
    for (int k = 0; k < j - 1; k++) {
        int q =
            countMatricsChain(a, i, k, m) + countMatricsChain(a, k + 1, j, m) + a[i] * a[k
                + 1] * a[j+1];
        if (q < m[i][j]) {
            m[i][j] = q;
        }
    }
    return m[i][j];
}

这时的时间复杂度也将是指数时间，与暴力穷尽解法区别并不大。其时间复杂度的分析如下：

这里假设 $T(n)$ 为计算 n 个矩阵相乘的最优解的所花费的时间，那么得到下面的公式：

$T(i) \geq 1$ ；
$T(n) \geq 1+ \sum_{k=1}^{n-1} {(T(k) + T(n-k) + 1)}，n >1$ ；

上面的公式，可以简化为：

$T(n) \geq 2·\sum_{k=1}^{n-1} {T(k)} + n$

下面还是根据代入法证明： $\Omega(2^n)$

这里先根据数据归纳法证明：对于所有的 $n\geq 1, T(n) \geq 2^{n-1}$ 都成立， $n=1$ 的时候很简单，当 $n\geq2$ 时，有：

$T(n) \geq 2·\sum_{k=1}^{n-1} {2^{k-1}} + n = 2·\sum_{k=0}^{n-2} {2^k} + n = 2(2^{n-1}-1)+n =2^n-2+n \geq 2^{n-1}$

这里可以看到，即使使用前面分析的递推公式去做，最后解法的时间复杂度依然是 $\Omega(2^n)$ ，那么有没有更好的解法呢？

动态规划求解

其实，关于动态规划问题，有一个非常明显的特点就是重叠子问题，上面的解法之所以时间复杂度那么高，一个重要的原因就是在使用递归方法时，有很多重复的计算，比如对于 $m(1,2)$ 这个值，其实 $m(1,3), m(1,4), m(1,5)...$ 都会用到，如何防止重复计算将是这个问题优化的一个重点，容易想到的方法就是使用空间换时间，在计算的过程中，额外申请一个二维数组（ $n*n$ ）保存 $m(i,j)$ 这个值，避免重复计算。

这个二维数据，也类似一个表格，下面的问题就是怎么填充这个表格（这个解法在算法导论上叫做自底向上的表格填充法）

关于表格这个说法，之前没有听说过，这是第一次听说，当时在做这道题的时候真的是有点脑子转不过来。

这里，我们举个例子，有一个表格如下，首先，我们有：

$m(i,j)=0, 如果 i=j$

所以表格对角线上的值均为0，因为是要求 $i \leq j$ 的，所以这个表格对角线下面的空格的值是不需要去填充的。

假设这里，我们要去求 $m(1,6)$ （图中红色的空格），根据 $min_{i \leq k < j} \begin{Bmatrix} m(i, k)+m(k+1, j) + a_i*a_{k+1}*a_{j+1} \end{Bmatrix}$ 这个公式，其实是需要下面几个空格的值（图中黄色空格的数值）：

$m(1,1), m(2,6)$ ；
$m(1,2), m(3,6)$ ；
$m(1,3), m(4,6)$ ；
$m(1,4), m(5,6)$ ；
$m(1,5), m(6,6)$ ；

这里是有一个规律的，那就是上面的值均在红线以下。而且最开始的对角线的值是有的，所以将对角线依次往右上方平移去计算，这样的话，在求 $m(i,j)$ 时，其所需要的值都是已知的。

表格填充法

右上角的深红色的点，就是我们要求解的最优解。

这部分的代码实现如下：

// a[0] 不使用，使用的是，1到 n+1
public int countMatricsChain(int[] a, int n) {
    int[][] m = new int[n + 1][n + 1];
    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] = Integer.MAX_VALUE;
            for (int k = i; k < j; k++) {
                int q = m[i][k] + m[k + 1][j] + a[i] * a[k + 1] * a[j + 1];
                if (q < m[i][j]) {
                    m[i][j] = q;
                }
            }
        }
    }
    return m[1][n];
}

关于这个时间复杂度比较容易求解，是 $\Omega(n^3)$ 。

平衡树描述问题

这里介绍另一种思路去理解这里问题，就是平衡树的思想（其实用平衡树去理解反而有些难度，不过明白这道题的解法之后再看就清晰很多），我们可以把上面的矩阵想象成一棵树最下面的子节点，这棵树的深度是 $n$ ，在处理矩阵相乘时，其实每一层都会合并一颗子树（且只会合并一棵树），过程如下图：

用平衡树的思想

其中， $A_5, A_3, A_1, A_7, B_3, B_1, D_1$ 分别为最佳的计算点。其实，上面的图只是其中一个切分过程，因为最佳计算点的值无法判断的。

自己在做的时候，开始有一个地方没明白：那就是将一颗树分为两棵树求解时，一棵树的最优解怎么转化为两个树最优解求解，因为当时想着还有最后一步合并，所以还会有一个变量，没有理解上面的拆分怎么成立，知道了最后的解法再去看就容易理解了，只要保证两颗子树是最优解，然后再加上那个变量（最后一步的矩阵相乘），遍历整个拆分点即可选择最小值即可，当时这一步没有转过来。

总结

其实这道题，算是一道常规题，如果对动态规划熟悉的，是很容易解决的，动态规划的题是有一套专门的解题方式/思想，说到底，还是自己的算法基础比较弱，之前并没有仔细看过动态规划的内容，希望自己能按照耗子叔的 ARTS 计划，保持每周学习一道算法题，坚持下去，来弥补这块的短板。

参考：

《算法导论》；