动态规划 矩阵链乘法

矩阵链乘问题描述

　　给定n个矩阵构成的一个链<A₁,A₂,A₃,.......A_n>，其中i=1,2,...n，矩阵A的维数为p_i-1p_i，对乘积 A₁A₂...A_n 以一种最小化标量乘法次数的方式进行加全部括号。

　　注意：在矩阵链乘问题中，实际上并没有把矩阵相乘，目的是确定一个具有最小代价的矩阵相乘顺序。找出这样一个结合顺序使得相乘的代价最低。

3、动态规划分析过程

1）最优加全部括号的结构

　　动态规划第一步是寻找一个最优的子结构。假设现在要计算A_iA_i+1....A_j的值，计算A_i...j过程当中肯定会存在某个k值（i<=k<j）将A_i...j分成两部分，使得A_i...j的计算量最小。分成两个子问题A_i...k和A_k+1...j，需要继续递归寻找这两个子问题的最优解。

　　有分析可以到最优子结构为：假设A_iA_i+1....A_j的一个最优加全括号把乘积在A_k和A_k+1之间分开，则A_i..k和A_k+1..j也都是最优加全括号的。

2）一个递归解

　　设m[i,j]为计算机矩阵A_i...j所需的标量乘法运算次数的最小值，对此计算A_1..n的最小代价就是m[1,n]。现在需要来递归定义m[i,j]，分两种情况进行讨论如下：

当i==j时：m[i,j] = 0，（此时只包含一个矩阵）

当i<j 时：从步骤1中需要寻找一个k（i≤k＜j）值，使得m[i,j] ＝min{m[i,k]+m[k+1,j]+p_i-1p_kp_j} （i≤k＜j）。

3）计算最优代价

　　虽然给出了递归解的过程，但是在实现的时候不采用递归实现，而是借助辅助空间，使用自底向上的表格进行实现。设矩阵Ai的维数为pi-1pi，i=1,2.....n。输入序列为：p=<p0，p1,...pn>，length[p] = n+1。使用m[n][n]保存m[i,j]的代价，s[n][n]保存计算m[i,j]时取得最优代价处k的值，最后可以用s中的记录构造一个最优解。书中给出了计算过程的伪代码。具体代码如下：

#include<iostream>
using namespace std;
 
#define maxvalue 1000000

void matrix_chain_order(int p[],int len,int **m,int **s)
{
    int n=len-1;
    //A[i][i]只有一个矩阵，所以相乘次数为0，即m[i][i]=0;
    for(int i=1;i<=n;i++)
    {
        m[i][i]=0;
    }
    for(int L=2;L<=n;L++)//l为chain length
    {
        for(int i=1;i<=n-L+1;i++) ///从第一矩阵开始算起，计算长度为L的最少代价,不是从0
        {
            int j=i+L-1;
            m[i][j]=maxvalue;
            for(int k=i;k<j;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_optimal_parents(int **s,int i,int j)
{
    if(i==j)
    {
        cout<<"A"<<i;
    }
    else
    {
        cout<<"(";
        print_optimal_parents(s,i,s[i][j]);
        print_optimal_parents(s,s[i][j]+1,j);
        cout<<")";
    }
}


int main()
{
    /*
    现有矩阵A1(30×35)、A2(35×15)、A3(15×5)、A4(5×10)、A5(10×20)、A6(20×25)，得到p=<30,35,15,5,10,20,25>
    */
    int p[]={30,35,15,5,10,20,25};
    int len=sizeof(p)/sizeof(int);
    cout<<len<<endl;
    int **m,**s;
    m=new int*[len];
    s=new int*[len];
    for(int i=1;i<len;i++)
    {
        m[i]=new int[len];
        s[i]=new int[len];
    }
     for(int i=1;i<len;i++)
     {
         for(int j=1;j<len;j++)
         {
             m[i][j]=0;
             s[i][j]=0;
         }
     }
     matrix_chain_order(p,len,m,s);
      cout<<"***********M为***********"<<endl;
     for(int i=1;i<len;i++)
     {
         for(int j=1;j<len;j++)
         {
             cout<<m[i][j]<<ends;
         }
         cout<<endl;
     }
     cout<<"***********S为***********"<<endl;
     for(int i=1;i<len;i++)
     {
         for(int j=1;j<len;j++)
         {
             cout<<s[i][j]<<ends;
         }
         cout<<endl;
     }

     print_optimal_parents(s,1,6);

}

输出：

 

我们假设矩阵为A1,A2,A3A4A5A6

数组p[7]恰好存储了Ai(维度为p(i-1)*p(i))

上面的代码我们可以按照下面的方式理解：

首先看对角线,矩阵元素长度为1时，m[i][j]=0;因为一个元素没有乘法为0.

长度为2时有2个元素，可以是(1,2),(2,3),(3,4),(4,5),(5,6)

以此类推，

长度为6时只有一种选择了(1,6).

可以看到m[i][j]形成的上三角。

如果采用递归进行实现，则需要指数级时间Ω(2ⁿ)，因为中间有些重复计算。递归是完全按照第二步得到的递归公式进行计算，递归实现如下所示：

int  recursive_matrix_chain(int p[],int len,int **m,int **s,int i ,int j)
{
    if(i==j)
    {
        m[i][j]=0;
    }
    else
    {
        int k;
        m[i][j]=maxvalue;
        for(k=i;k<j;k++)
        {
            int tmp=recursive_matrix_chain(p, len,m,s,i,k)+recursive_matrix_chain(p,len,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];
}

只需改下调用方法即可：

// matrix_chain_order(p,len,m,s); 改为

cout<< recursive_matrix_chain(p,len,m,s,1,6)<<endl;

改进递归：

递归重复了很多计算，我们可以做一个备忘录，其思想就是备忘（memorize）原问题的自然但低效的递归算法。

加入备忘的递归算法为每一个子问题的解在表中记录一个表项，开始时，美国表现最初都含有一个特殊的值，以表示该值有待填入。当在递归算法的执行中第一次遇到一个子问题时，就计算它的解并填入表中以后。以后每次遇到子问题时，只要查看返回表中先前填入的值即可。

int  memoized_matrix_chain(int p[],int len,int **m,int **s)
{
    int n=len-1;
    for(int i=1;i<=n;i++)
        for(int j=i;j<=n;j++)
            m[i][j]=maxvalue;
    return lookup_chain(p,len,m,s,1,n);
}
int lookup_chain(int p[],int len,int **m,int **s,int i,int j)
{
    if(m[i][j]<maxvalue)
    {
        return m[i][j];
    }
    if(i==j)
    {
        m[i][j]=0;
    }
    else
    {
        for(int k=i;k<=j-1;k++)
        {
            int q=lookup_chain(p, len,m,s,i,k)+lookup_chain(p,len,m,s,k+1,j)+
                    p[i-1]*p[k]*p[j];
            if(q<m[i][j])
            {
                m[i][j]=q;
                s[i][j]=k;
            }
        }
    }

    return m[i][j];
}

像动态规划的复杂度一样，

memoized_matrix_chain的复杂度也是O(n^3)d .

参考了：http://www.cnblogs.com/Anker/archive/2013/03/10/2952475.html

 http://www.cnblogs.com/weixliu/archive/2012/12/25/2832772.html