动态规划实例（二）求解矩阵链乘法

问题

矩阵链乘法问题(matrix-chain multiplication problem)可描述如下：给定n个矩阵的链<A₁, A₂, ..., A_n>，矩阵A_i的规模为p_i-1×p_i (1<=i<=n)，求完全括号化方案，使得计算乘积 A₁A₂...A_n所需标量乘法次数最少 —— 完全括号化(fully parenthesized)：它是单一矩阵，或者是两个完全括号化的矩阵乘积链的积，且已外加括号。例如，如果矩阵链为<A₁, A₂, A₃, A₄>，则共有5种完全括号化的矩阵乘积链：

(A1(A2(A3A4)))

(A1((A2A3)A4))

((A1A2)(A3A4))

((A1(A2A3))A4)

((A1(A2A3))A4)

对矩阵链加括号的方式会对乘积运算的代价产生巨大影响。我们先来分析两个矩阵相乘的代价。下面的伪代码的给出了两个矩阵相乘的标准算法，属性rows和columns是矩阵的行数和列数。两个矩阵A和B只有相容(compatible)，即A的列数等于B的行数时，才能相乘。

MATRIX-MULTIPKLY(A,B)
if A.columns≠B.rows 
    error "incompatible dimensions"
else let C be a new A.rows×B.columns matrix
    for i = 1 to A.rows
         for j = 1 to B.columns
          c(ij)=0
           for k = 1 to A.columns
             c(ij)=c(ij)+a(ik)*b(kj)
return C

为了说明不同的加括号方式会导致不同的计算代价，我们以矩阵链<A1,A2,A3>为例 —— 假设三个矩阵的规模分别为10×100、100×5和5×50。

• 如果按照((A1A2)A3)的顺序计算，为计算A1A2（规模10×5），需要做10*100*5=5000次标量乘法，再与A3相乘又需要做10*5*50=2500次标量乘法，共需7500次标量乘法。
• 如果按照(A1(A2A3))的顺序计算，为计算A2A3（规模100×50），需要做100*5*50=25000次标量乘法，再与A1相乘又需要做10*100*50=50000次标量乘法，共需75000次标量乘法。

第一种顺序计算要比第二种顺序计算快10倍。

思路

按照动态规划的四个步骤来进行，清楚地展示针对本问题每个步骤应如何做。

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

假设A_iA_i+1...A_j的最优括号方案的分割点在A_k和A_k+1之间。那么，继续对“前缀”子链A_iA_i+1...A_k进行括号化时，我们应该直接采用独立求解它时所得的最优方案。为了构造一个矩阵链乘法问题实例的最优解，我们可以将问题划分为两个子问题(A_iA_i+1...A_k和A_k+1A_k+2...A_j的最优括号化问题)，求出子问题实例的最优解，然后将子问题的最优解组合起来。我们必须保证在确定分割点时，已经考察了所有可能的划分点，这样就可以保证不会遗漏最优解。

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

A_iA_i+1...A_j的最小代价括号化方案的递归求解公式为：

1. 若i=j，m[i,j]=0
2. 若i<j，m[i,j]=min{m[i,k]+m[k+1,j]+p_i-1p_kp_j}, i<=k<j

m[i,j]的值给出了子问题最优解的代价，但它并未提供足够的信息来构造最优解。为此，我们用s[i,j]保存最优括号化方案的分割点位置k，即使得m[i,j]=m[i,k]+[k+1,j]+p_i-1p_kp_j成立的k值。

步骤3：计算最优代价 

现在，我们可以很容易地基于递归公式写出一个递归算法，但递归算法是指数时间的，因为递归算法会在递归调用树的不同分支中多次遇到同一个子问题，并不必检查若有括号化方案的暴力搜索方法更好！

我们采用自底向上表格法代替递归算法来计算最优代价 —— 此过程假定矩阵Ai的规模为p_i-1×p_i (i=1,2,...,n)。它的输入是一个序列p=<p₀, p₁, ..., p_n>，其长度为p.length=n+1。过程用一个辅助表m[1..n, 1..n]来保存代价m[i, j]，用另一个辅助表s[1..n-1, 2..n](s[1,2]..s[n-1,n]这里i<j)记录最优值m[i,j]对应的分割点k。我们就可以利用表s构造最优解。对于矩阵A_iA_i+1...A_j最优括号化的子问题，我们认为其规模为链的长度j-i+1。因为j-i+1个矩阵链相乘的最优计算代价m[i,j]只依赖于那么少于j-i+1个矩阵链相乘的最优计算代价。因此，算法应该按长度递增的顺序求解矩阵链括号化问题，并按对应的顺序填写表m。

步骤4：构造最优解 

在上一个步骤中可以求出计算矩阵链乘积所需的最少标量乘法运算次数，但它并未直接指出如何进行这种最优代价的矩阵链乘法计算。表s[i, j]记录了一个k值，指出A_iA_i+1...A_j的最优括号化方案的分割点应在A_k和A_k+1之间。

因此，我们A(1..n)的最优计算方案中最后一次矩阵乘法运算应该是以s[1, n]为分界的A(1..s[1, n])*A(s[1, n]+1..n)。我们可以用相同的方法递归地求出更早的矩阵乘法的具体计算过程，因为s[1,s[1,n]]指出了计算A(1..s[1,n])时应进行的最后一次矩阵乘法运行；s[s[1,n]+1,n]指出了计算A(s[1, n]+1..n)时应进行的最后一次矩阵乘法运算。

C++实现

#include<iostream>
using namespace std;
const int INT_MAX=2147483647;
int const M=7;
void MATRIX_CHAIN_ORDER(int *p,int Length,int m[][M],int s[][M])
{
    int q,n=Length-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;         /* 等价于 l=j-i+1 */
            m[i][j]=INT_MAX;
            for(int k=i;k<=j-1;k++)
            {
                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 PRINT_OPTIMAL_PARENS(int s[][M],int i,int j)
{
    if(i == j) cout<<"A"<<i;
    else
    {
        cout<<"(";
        PRINT_OPTIMAL_PARENS(s,i,s[i][j]);
        PRINT_OPTIMAL_PARENS(s,s[i][j]+1,j);
        cout<<")";
    }
}
int main()
{
   int p[M]={30,35,15,5,10,20,25};
   int m[M][M],s[M][M];
   MATRIX_CHAIN_ORDER(p,M,m,s);
   cout<<"当n=6时最优解为: \n"<<m[1][6];
   cout<<"\n括号化方案为：\n"; 
   PRINT_OPTIMAL_PARENS(s,1,6);
   return 0;
}

代码引自：【动态规划】矩阵链乘法_小与米的博客-CSDN博客_矩阵链乘法