动态规划之矩阵链乘法

上一篇介绍了动态规划的基本概念，又介绍了一个典型的例子链条分割，这一篇来介绍另一个例子-----矩阵链乘法

给定一个n个矩阵的序列，我们希望计算它的乘积 A1A2A3...An，我们可以用括号标定计算次序，然后再利用标准的矩阵乘法算法来计算，由于矩阵乘法满足结合律，所以任何一种计算顺序最后的计算结果都是相同的。

完全括号化的矩阵乘积链：它是单一矩阵，或者两个完全括号化的矩阵乘积链的积，并且已经外加括号；例如矩阵链（A1，A2，A3，A4），则有下边5种完全括号化的矩阵链：

对矩阵链加括号虽然对计算结果没影响，但是计算过程中所花费的代价却是影响巨大。我们先来做简单的分析，首先A，B矩阵能相乘的条件是两个矩阵相容（即A矩阵的列数等于B矩阵的行数），设A是n*m阶，B是m*r阶，那么要进行n*m*r次标量乘法运算；举个例子：设A1，A2，A3三个矩阵分别为10x100,100x5,5x50阶，若按（A1（A2A3））顺序计算，那么总共进行（100x5x50）+（10x100x50）=75000次标量乘法运算；若按（（A1A2）A3）顺序，那么总共进行（10x100x5）+（10x5x50）=7500次标量乘法运算，第二种顺序比第一种运算快10倍。

引入矩阵乘法链问题：给定n个矩阵的链（A1，A2，...An），矩阵Ai的规模pi-1xpi，求完全括号化方案，使得A1A2...An的乘积所需的标量乘法次数最小。

在此之前先来考虑一下完全括号化方案的总数，用P(n)表示方案总数，当n=1时，只有一个矩阵，只有一种方案；n>=2时，这样考虑，完全括号化的矩阵乘积可以看成两个完全括号化的部分积的形式，总结成下边的式子：

下边根据动态规划的4个步骤来设计算法：

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

用Ai...j表示AiAi+1...Aj的乘积，为了进行完全括号化，必须在Ak,Ak+1之间将链断开，即对一个整数k，我们先计算Ai...k，Ai+1...j，再计算两者的乘积，此方案的代价是Ai...k，Ak+1...j的计算代价，加上两者相乘的计算代价。有一个重要的结论：在原问题Ai...Aj的最优括号化方案中，对子链Ai...Ak，Ak+1...Aj括号化的方案就是原问题的括号化方案。

2、一个递归求解方案：

用m[i,j]表示计算Ai...Aj所需的标量乘法次数的最小值，那么m[1,n]就表示整个问题的解。这样来定义递归式，当i=j时，输入只有一个矩阵，即m[i,i]=0，当i<j时，可以在Ak和Ak+1之间划分，那么m[i,j]就是Ai...Ak的次数和Ak+1...Aj的次数，再加上两者相乘次数的最小值。

最后再加一个量s[i,j]保存划分位置。

3、计算代价，给出实现：

自底向上表格法

//自底向上表格法
void Min_ChainOrder_down2up(DataType p[],int n,DataType m[][N],DataType s[][N])
{
    for(int l=2;l<=n;l++)
    {
        for(int i=1;i<=n-l+1;i++)
        {
            int j=l+i-1;
            m[i][j]=MAX_LIMITS;
            for(int k=i;k<=j-1;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;
                }
            }
        }
    }
}

算法时间复杂度O(n^3)，并且需要O(n^2)的内存来维护m和s这两个表

4、构造最优解：

void print_result(DataType s[][N],int i,int j)
{
    if(i==j)
        cout<<"A"<<i;
    else
    {
        cout<<"(";
        print_result(s,i,s[i][j]);
        print_result(s,s[i][j]+1,j);
        cout<<")";
    }
}

下边给出规模为6的测试代码，关于此算法的详细分析和理解另开一帖详细解释

#include <iostream>
#define N 7
#define MAX_LIMITS 100000
typedef int DataType;
using namespace std;
//自底向上表格法
void Min_ChainOrder_down2up(DataType p[],int n,DataType m[][N],DataType s[][N])
{
    for(int l=2;l<=n;l++)
    {
        for(int i=1;i<=n-l+1;i++)
        {
            int j=l+i-1;
            m[i][j]=MAX_LIMITS;
            for(int k=i;k<=j-1;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 print_result(DataType s[][N],int i,int j)
{
    if(i==j)
        cout<<"A"<<i;
    else
    {
        cout<<"(";
        print_result(s,i,s[i][j]);
        print_result(s,s[i][j]+1,j);
        cout<<")";
    }
}
int main()
{
    //6个矩阵，规模为6，N=7
    DataType p[N]={30,35,15,5,10,20,25};
    int n=N-1;
    //
    DataType m[N][N]={0};
    DataType s[N][N]={0};
    Min_ChainOrder_down2up(p,n,m,s);
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            cout<<"m["<<i<<"]["<<j<<"]="<<m[i][j]<<endl;
            cout<<"s["<<i<<"]["<<j<<"]="<<s[i][j]<<endl;
        }
    }
    cout<<"result:"<<endl;
    print_result(s,1,n);
    return 0;
}