算法设计与分析——动态规划（三）：矩阵链乘法

分类目录：《算法设计与分析》总目录
相关文章：
· 动态规划（一）：基础知识
· 动态规划（二）：钢条切割
· 动态规划（三）：矩阵链乘法
· 动态规划（四）：动态规划详解
· 动态规划（五）：最长公共子序列

---

这个例子是求解矩阵链相乘问题的动态规划算法。给定一个n个矩阵的序列（矩阵链）$(A_1.A_2.\cdots ,A_n)$，我们希望计算它们的乘积$A_1{A}_2\cdots A_n$其中，我们可以先用括号明确计算次序，然后利用标准的矩阵相乘算法进行计算。由于矩阵乘法满足结合律，因此任何加括号的方法都会得到相同的计算结果。我们称有如下性质的矩阵乘积链为完全括号化的：它是单一矩阵，或者是两个完全括号化的矩阵乘积链的积，且已外加括号。例如，如果矩阵链为$(A_1,A_2,A_3)$，则共有2种完全括号化的矩阵乘积链：$((A_1{A}_2)A_3)$、$(A_1(A_2{A}_3))$。

对矩阵链加括号的方式会对乘积运算的代价产生巨大影响。我们先来分析两个矩阵相乘的代价：

def matrix_multiply(A:np.array, B:np.array):
    if A.shape[1] != B.shape[0]:
        print('Incompatible Dimensions')
    else:
        ans = np.zeros([A.shape[0], B.shape[1]])
        for row_index in range(len(A)):
            for column_index in range(len(B[0])):
                ans[row_index, column_index] = sum(A[row_index] * B[:, column_index])
    return ans

两个矩阵$A$和$B$只有相容，即$A$的列数等于$B$的行数时，才能相乘。如果$A$是$p\times q$的矩阵，$B$是$q\times r$的矩阵，那么乘积$C$是$p\times r$的矩阵。计算$C$所需时间由$pqr$决定。

我们以矩阵链$(A_1,A_2,A_3)$相乘为例，来说明不同的加括号方式会导致不同的计算代价假设三个矩阵的规模分别为$10\times 100$、$100\times 5$和$5\times 50$。如果按$((A_1{A}_2)A_3)$的顺序计算，为计算$A_1{A}_2$，需要做$10\times 100\times 5=5000$次标量乘法，再与$A_3$相乘又需要做$10\times 5\times 50=2500$ 次标量乘法，共需$7500$次标量乘法。如果按$(A_1(A_2{A}_3))$的顺序，计算$A_2{A}_3$，需$100\times 5\times 50=25000$次标量乘法，$A_1$再与之相乘又需$10\times 100\times 50=50000$次标量乘法共需$75000$次标量乘法。因此，按第一种顺序计算矩阵链乘积要比第二种顺序快10倍。

所以，矩阵链乘法问题可描述如下：给定$n$个矩阵的链）$(A_1.A_2.\cdots ,A_n)$，矩阵$A_i$的规模为$p_{i-1}\times p_i$，求完全括号化方案，使得计算乘积$A_1{A}_2\cdots A_n$所需标量乘法次数最少。

注意，求解矩阵链乘法问题并不是要真正进行矩阵相乘运算，我们的目标只是确定代价最低的计算顺序。确定最优计算顺序所花费的时间通常要比随后真正进行矩阵相乘所节省的时间要少。

在用动态规划方法求解矩阵链乘法问题之前，我们先来证明穷举所有可能的括号化方案不会产生一个高效的算法。对一个$n$个矩阵的链，令$P(n)$表示可供选择的括号化方案的数量。当$n=1$时，由于只有一个矩阵，因此只有一种完全括号化方案。当$n\ge 2$时，完全括号化的矩阵乘积可描述为两个完全括号化的部分积相乘的形式，而两个部分积的划分点在第$k$个矩阵和第$k+1$个矩阵之间，$k$为$1,2,\cdots ,n-1$中的任意一个值。因此，当$n\ge 2$时我们可以得到如下递归公式：
$$p(n)=\sum _{k=1}^{n-1}p(k)P(n-k)$$

因此，括号化方案的数量与$n$呈指数关系，通过暴力搜索穷尽所有可能的括号化方案来寻找最优方案，是个糟糕的策略。

应用动态规划方法

下面用动态规划方法来求解矩阵链的最优括号化方案，我们还是按照《算法设计与分析——动态规划（一）：基础知识》开头提出的4个步骤进行：

1. 刻画一个最优解的结构特征。
2. 递归地定义最优解的值
3. 计算最优解的值，通常采用自底向上的方法。
4. 利用计算出的信息构造一个最优解。

我们按顺序进行这几个步骤，清楚地展示针对本问题每个步骤应如何做。

步骤1：最优括号化方案的结构特征

动态规划方法的第一步是寻找最优子结构，然后就可以利用这种子结构从子问题的最优解构造出原问题的最优解。在矩阵链乘法问题中，此步骤的做法如下所述。为方便起见，我们用符号$A_{i,j}(i\le j)$表示$A_i{A}_{i+1}\cdots A_j$乘积的结果矩阵。可以看出，如果问题是非平凡的，即$i\le j$，那么为了对$A_i{A}_{i+1}\cdots A_j$进行括号化，我们就必须在某个$A_k$和$A_{k+1}$之间将矩阵链划分开。也就是说，对某个整数$k$，我们首先计算矩阵$A_{i,k}$和$A_{k+1,j}$，然后再计算它们的乘积得到最终结果$A_{i,j}$。此方案的计算代价等于矩阵$A_{i,k}$的计算代价，加上矩阵$A_{k+1,j}$，的计算代价，再加上两者相乘的计算代价。

下面我们给出本问题的最优子结构。假设$A_i{A}_{i+1}\cdots A_j$的最优括号化方案的分割点在$A_k$和$A_{k+1}$之间。那么，继续对“前缀”子链$A_{i,k}$进行括号化时，我们应该直接采用独立求解它时所得的最优方案。

我们已经看到，一个非平凡的矩阵链乘法问题实例的任何解都需要划分链，而任何最优解都是由子问题实例的最优解构成的。因此，为了构造一个矩阵链乘法问题实例的最优解，我们可以将问题划分为两个子问题$A_i{A}_{i+1}\cdots A_k$和$A_{k+1}{A}_{k+2}\cdots A_j$，的最优括号化问题，求出子问题实例的最优解，然后将子问题的最优解组合起来。我们必须保证在确定分割点时，已经考察了所有可能的划分点，这样就可以保证不会遗漏最优解。

步骤2：一个递归求解方案

下面用子问题的最优解来递归地定义原问题最优解的代价。令$m[i,j]$表示计算矩阵$A_{i,j}$所需标量乘法次数的最小值，那么，原问题的最优解——计算$A_{i,n}$所需的最低代价就是$m[1,n]$。

我们可以递归定义$m[i,j]$如下。对于$i=j$时的平凡问题，矩阵链只包含唯一的矩阵$A_i$，因此不需要做任何标量乘法运算。所以，对所有$i=1,2,\cdots ,n$，有$m[i,i]=0$，若$i<j$，我们利用步骤1中得到的最优子结构来计算$m[i,j]$。我们假设$A_i{A}_{i+1}\cdots A_j$的最优括号化方案的分割点在矩阵$A_k$和$A_{k+1}$之间。那么$m[i,j]$就等于计算$A_k$和$A_{k+1}$的代价加上两者相乘的代价的最小值。由于矩阵$A_i$的规模为$p_{i-1}\times p_i$，易知$A_{i,k}$和$A_{k+1,j}$相乘的代价为$p_{i-1}{p}_k{p}_j$次标量乘法运算。因此，我们得到

$$m[i,j]=m[i,k]+m[k+1,j]+p_{i-1}{p}_k{p}_j(i<j)$$

此递归公式假定最优分割点$k$是已知的，但实际上我们是不知道的。不过，$k$只有$j-i$种可能的取值，即$k=i,i+1,\cdots ，j-1$，由于最优分割点必在其中，我们只需检查所有可能情况，找到最优者即可。

$m[i,j]$的值给出了子问题最优解的代价，但它并未提供足够的信息来构造最优解。为此，我们用$s[i,i]$保存$A_i{A}_{i+1}\cdots A_j$最优括号化方案的分割点位置$k$，即使得$m[i,j]=m[i,k]+m[k+1,j]+p_{i-1}{p}_k{p}_j(i<j)$成立的$k$值。

步骤3：计算最优代价

现在，我们可以很容易地基于递归公式写出一个递归算法，来计算$A_1{A}_2\cdots A_n$相乘的最小代价$m[1,ni]$。像我们在《算法设计与分析——动态规划（二）：钢条切割》中所看到的那样，此递归算法是指数时间的，并不比检查所有括号化方案的暴力搜索方法更好。

注意到，我们需要求解的不同子问题的数目是相对较少的：每对满足$1\le i\le j\le n$的$i$和$j$对应一个唯一的子问题，共有$O(2^n)$个。递归算法会在递归调用树的不同分支中多次遇到同一个子问题。这种子问题重叠的性质是应用动态规划的另一个$s[i,j]$来标识（$m[i,j]$标识是最优子结构）。

我们采用自底向上表格法代替基于公式递归算法来计算最优代价。下面给出的过程 matrix_chain_order(p:list)实现了自底向上表格法。此过程假定矩阵$A_i$的规模为$p_{i-1}\times p_i$。它的输入是一个序列$p_0,p_1,\cdots ,p_n$，其长度为$n+1$，过程用一个辅助表$m[i,j]$来保存，用另一个辅助表$s[i,j]$记录最优值$m[i,j]$对应的分割点$k$。我们就可以利用表$s$构造最优解。

def matrix_chain_order(p:list):
    m = np.array([[-1] * (len(p) - 1)] * (len(p) - 1))
    s = np.array([[-1] * (len(p) - 1)] * (len(p) - 1))
    
    for matrix_length in range(1, len(p)):
        for i in range(len(p) - matrix_length):
            if matrix_length == 1:
                m[i, i] = 0
                s[i, i] = i
            else:
                j = i + matrix_length - 1
                cost = 1e10
                for k in range(i,j):
                    _ = sum(m[i, k] + m[k + 1, j] + p[i - 1] * p[k] * p[j])
                    if _ < cost:
                        cost = _
                        m[i, j] = cost
                        s[i, j] = k
    return m, s

简单分析 matrix_chain_order(p:list)的嵌套循环结构，可以看到算法的运行时间为$O(n^3)$。此外，算法还需要$O(n^2)$的内存空间来保存表m和s。因此， matrix_chain_order(p:list)比起穷举所有可能的括号化方案来寻找最优解的指数阶算法要高效得多。

步骤4：构造最优解

虽然 matrix_chain_order(p:list)求出了计算矩阵链乘积所需的最少标量乘法运算次数，但它并未直接指出如何进行这种最优代价的矩阵链乘法计算。表$s[i,j]$记录了构造最优解所需的信息。每个表$s[i,j]$记录了一个$k$值，指出$A_i{A}_{i+1}\cdots A_j$的最优括号化方案的分割点应在$A_k$和$A_{k+1}$之间。因此，我们通过表$s[i,j]$可以很轻易的得出整个矩阵链的切分点。