动态规划之矩阵连乘

【问题描述】

给定n个矩阵｛A1,A2,…,An｝，其中Ai与Ai+1是可乘的，i=1，2…，n-1。如何确定计算矩阵连乘积的计算次序，使得依此次序计算矩阵连乘积需要的数乘次数最少。例如，给定三个连乘矩阵{A1，A2，A3}的维数分别是10*100，100*5和5*50，采用（A1A2）A3，乘法次数为10*100*5+10*5*50=7500次，而采用A1（A2A3），乘法次数为100*5*50+10*100*50=75000次乘法，显然，最好的次序是（A1A2)A3，乘法次数为7500次。

分析：
矩阵链乘法问题描述：
给定由n个矩阵构成的序列｛A1，A2，...，An}，对乘积A1A2...An，找到最小化乘法次数的加括号方法。

1）寻找最优子结构
此问题最难的地方在于找到最优子结构。对乘积A1A2...An的任意加括号方法都会将序列在某个地方分成两部分，也就是最后一次乘法计算的地方，我们将这个位置记为k，也就是说首先计算A1...Ak和Ak+1...An，然后再将这两部分的结果相乘。
最优子结构如下：假设A1A2...An的一个最优加括号把乘积在Ak和Ak+1间分开，则前缀子链A1...Ak的加括号方式必定为A1...Ak的一个最优加括号，后缀子链同理。
一开始并不知道k的确切位置，需要遍历所有位置以保证找到合适的k来分割乘积。

2）构造递归解
设m[i,j]为矩阵链Ai...Aj的最优解的代价，则

3）构建辅助表，解决重叠子问题
从第二步的递归式可以发现解的过程中会有很多重叠子问题，可以用一个nXn维的辅助表m[n][n] s[n][n]分别表示最优乘积代价及其分割位置k 。
辅助表s[n][n]可以由2种方法构造，一种是自底向上填表构建，该方法要求按照递增的方式逐步填写子问题的解，也就是先计算长度为2的所有矩阵链的解，然后计算长度3的矩阵链，直到长度n；另一种是自顶向下填表的备忘录法，该方法将表的每个元素初始化为某特殊值(本问题中可以将最优乘积代价设置为一极大值)，以表示待计算，在递归的过程中逐个填入遇到的子问题的解。

第一种 自底向上 的方式

#include<iostream>
using namespace std;

const int N=7;
//p为矩阵链，p[0],p[1]代表第一个矩阵，p[1],p[2]代表第二个矩阵，length为p的长度
//所以如果有六个矩阵，length=7，m为存储最优结果的二维矩阵，s为存储选择最优结果路线的
//二维矩阵
void MatrixChainOrder(int *p,int m[N][N],int s[N][N],int length)
{
    int n=length-1;
    int l,i,j,k,q=0;
    //m[i][i]只有一个矩阵，所以相乘次数为0，即m[i][i]=0;
    for(i=1;i<length;i++)
    {
        m[i][i]=0;
    }
    //l表示矩阵链的长度
    // l=2时，计算 m[i,i+1],i=1,2,...,n-1 (长度l=2的链的最小代价)
    for(l=2;l<=n;l++)
    {
        for(i=1;i<=n-l+1;i++)
        {
            j=i+l-1; //以i为起始位置，j为长度为l的链的末位，
            m[i][j]=0x7fffffff;
            //k从i到j-1,以k为位置划分
            for(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;
                }
            }
        }
    }
    cout << m[1][N-1] << endl;
}
void PrintAnswer(int s[N][N],int i,int j)
{
    if(i==j)
    {
        cout<<"A"<<i;
    }
    else
    {
        cout<<"(";
        PrintAnswer(s,i,s[i][j]);
        PrintAnswer(s,s[i][j]+1,j);
        cout<<")";
    }
}
int main()
{
    int p[N]={30,35,15,5,10,20,25};
    int m[N][N],s[N][N];
    MatrixChainOrder(p,m,s,N);
    PrintAnswer(s,1,N-1);
    return 0;
}

对于 p={30 35 15 5 10 20 25}:

计算顺序为：

第二种 自顶向下 （备忘录）的方式

#include<iostream>
using namespace std;

const int N=7;

int MatrixChainOrder2(int *p,int m[N][N],int s[N][N],int i, int j)
{
    
    if (i == j)
    {
        return 0;
    }
    if (m[i][j] < 0x7fffffff)
    {
        return m[i][j];
    }

    for (int k=i; k<j; ++k)
    {
        
        int tmp = MatrixChainOrder2(p,m,s,i,k) + MatrixChainOrder2(p,m,s,k+1,j) + p[i-1]*p[k]*p[j];
        if (tmp < m[i][j])
        {
            m[i][j] = tmp;
            s[i][j] = k;
        }
    }
    return m[i][j];
}


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


int main()
{
    int p[N]={30,35,15,5,10,20,25};
    int m[N][N],s[N][N];

    for (int i=0; i<N; ++i)
    {
        for (int j=0; j<N; ++j)
            m[i][j] = 0x7fffffff;
    }
    cout << MatrixChainOrder2(p,m,s,1,6)<<endl;
    PrintAnswer(s,1,N-1);
    return 0;
}

矩阵连乘计算次序问题的最优解包含着其子问题的最优解。这种性质称为最优子结构性质。 

•在分析问题的最优子结构性质时，所用的方法具有普遍性：首先假设由问题的最优解导出的子问题的解不是最优的，然后再设法说明在这个假设下可构造出比原问题最优解更好的解，从而导致矛盾。 

•利用问题的最优子结构性质，以自底向上的方式递归地从子问题的最优解逐步构造出整个问题的最优解。最优子结构是问题能用动态规划算法求解的前提。

•递归算法求解问题时，每次产生的子问题并不总是新问题，有些子问题被反复计算多次。这种性质称为子问题的重叠性质。 

•动态规划算法，对每一个子问题只解一次，而后将其解保存在一个表格中，当再次需要解此子问题时，只是简单地用常数时间查看一下结果。 

•通常不同的子问题个数随问题的大小呈多项式增长。因此用动态规划算法只需要多项式时间，从而获得较高的解题效率。