0x07算法设计与分析复习（二）：算法设计策略-动态规划法3

参考书籍：算法设计与分析——C++语言描述（第二版）

算法设计策略-动态规划法

矩阵连乘

问题描述

给定n个矩阵{A0,A1,⋯,An−1}，其中Ai(i=0,⋯,n−1)的维数为pi×pi+1，并且Ai与Ai+1是可乘的。考虑这n个矩阵的连乘积A0A1⋯An−1，由于矩阵乘法满足结合律，所以计算矩阵的连乘可以与许多不同的计算次序。矩阵连乘问题是确定计算矩阵连乘积的计算次序，使得按照这一次序计算矩阵连乘积，需要的“数乘”次数最少。

设有矩阵A和B，A是m×n矩阵，B是n×p矩阵，则A和B是可乘的。乘积矩阵D=AB是m×p矩阵，它的元素dij为：

dij=∑k=0n−1aik∗bkj  0≤i<m,0≤j<p

矩阵A和B相乘的数乘（两元素相乘）次数是：n×m×p。通常以矩阵乘法中需执行的数乘次数作为两矩阵相乘的计算量。
矩阵连乘的计算次序可以用加括号的方式来确定，一旦一个矩阵连乘的计算次序完全确定，也称该连乘积已完全加括号（fully parenthesized）。

完全加括号的矩阵连乘积可递归地定义为：

• 单个矩阵是完全加括号的；
• 矩阵连乘积A是完全加括号的，则A可表示为两个完全加括号的矩阵连乘积B和C的乘积并加括号，即A=(BC)。

设p(n)是n个矩阵连乘时可能的加括号的方案数目，它们决定不同的计算次序。假定先将连乘的矩阵序列A0A1⋯An−1分解成两个矩阵子序列A0A1⋯Ak和Ak+1Ak+2⋯An−1,0≤k<n−1，然后分别对两个子序列完全加括号，最后加上最外层括号，得到原矩阵序列的完全加括号形式，由此可以得到关于p(n)的递推关系式：

p(n)=⎧⎩⎨⎪⎪⎪⎪1∑k=1n−1p(k)p(n−k)n=1n≥2

此递推式的解是卡特朗（Catalan）数列：P(n)=C(n−1)。式中：

C(n)=1n+1(2nn)=Ω(4nn3/2)

动态规划法求解

最优子结构

矩阵连乘AiAi+1⋯Aj简记为A[i:j]，i≤j。于是矩阵连乘A0A1⋯An−1可记作A[0:n−1]。将这一计算次序在矩阵Ak和Ak+1，0≤k<n−1之间断开，则其相应的完全加括号形式为((A0A1⋯Ak)(Ak+1Ak+2⋯An−1))。可先分别计算A[0:k]和A[k+1:n−1]，然后将两个连乘积再相乘得到A[0:n−1]。

矩阵连乘A[0:n−1]的最优计算次序的计算量等于A[0:k]和A[k+1:n−1]两者的最有计算次序的计算量之和，再加上A[0:k]和A[k+1:n−1]相乘的计算量。如果两个矩阵子序列的计算次序不是最优的，则原矩阵的计算次序也不可能是最优的。所以，矩阵连乘问题的最优解具有最优子结构特性。

最优解的递推关系

定义一个二维数组m，用来保存矩阵连乘时所需的最少计算量。
设m是二维数组，m[i][j]定义为：计算矩阵连乘A[i:j](0≤i≤j≤n−1)所需的最少乘次数。

当i=j时，A[i:j]=Ai是单一矩阵，无需计算因此m[i][i]=0,i=0,1,⋯,n−1。
当i<j时，假定在Ak和Ak+1之间分解，则通过先分别计算两个子序列连乘后再相乘的方式计算A[i:j]的计算量m[i][j]为：

m[i][j]=m[i][k]+m[k+1][j]+pipk+1pj+1

这里Ai的维数是pi×pi+1，m[i][k]和m[k+1][j]分别为A[i:k]和A[k+1:j]的最少计算量，pipk+1pj+1是两个子序列连乘的乘积矩阵再相乘的计算量。由于k=i,i+1,⋯,j−1，因此可得到如下的递推式：

m[i][j]=⎧⎩⎨0mini≤k<j{m[i][k]+m[k+1][j]+pipk+1pj+1}i=ji<j

在上式中，i≤j所以只需计算m的上三角部分元素。使上式取最小值的k就是A[i:j]的最优计算次序中的断开位置，如果将此断开位置k保存在s[i][j]中，就可以计算得到最优解值之后，由s构造出相应的最优解。

重叠子问题

由于不同的有序对(i,j),0≤i≤j≤n−1对应于不同的子问题，因此，不同子问题的个数最多只有

(n2)+n=Θ(n2)

可以证明，直接采用递推计算式，具有指数事件，这与检查每一种加括号方式相同。在递归计算时，许多子问题被重复计算多次。这也是采用动态规划法解决此问题的原因之一，用动态规划法求解此问题，可以采用自底向上的方式进行计算。在计算过程中，保存已求得的子问题的解。每个子问题只计算一次，以后可以直接使用，从而避免重复计算。

矩阵连乘算法

//矩阵连乘算法
 
class MatrixChain
{
public:
    MatrixChain(int mSize,int *q);//创建二维数组m和s，一维数组p，并初始化
    int MChain();//一般动态规划算法求最优解值，计算最优解值m[0][n-1]
    int LookipChain();//备忘录方法计算最优解值
    void Traceback();//构造最优解的共有函数，从s构造最优解，机构造矩阵连乘序列的完全加括号形式
    ...
 
private:
    void Traceback(int i, int j);//构造最优解的私有递归函数
    int LookupChain(int i,int j);//备忘录方法私有递归
    int *p,**m,**s,n;
};
int MatrixChain::MChain()
{
    //求A[0:n-1]的最优解值
    for(int i=0;i<n;i++)
        m[i][i]=0;
    for(int r=2;r<=n;r++)
        for(int i=0;i<=n-r;i++){
            int j=i+r-1;
            m[i][j]=m[i+1][j]+p[i]*p[i+1]*p[j+1];//m[i][j]的初值
            s[i][j]=i;
            for(int k=i+1;k<j;k++){
                int t=m[i][k]+m[k+1][j]+p[i]*p[k+1]*p[j+1];
                if(t<m[i][j]){
                    m[i][j]=t;
                    s[i][j]=k;
                }
            }
        }
    return m[0][n-1];
}
void MatrixChain::Traceback(int i,int j)
{
    if(i==j){
        cout << 'A' << i;
        return ;
    }
    if(i<s[i][j])   cout << '(';
    Traceback(i,s[i][j]);
    if(i<s[i][j])   cout << ')';
 
    if(s[i][j]+1<j) cout << '(';
    Traceback(s[i][j]+1,j);
    if(s[i][j]+1<j) cout << ')';
}
void MatrixChain::Traceback()
{
    cout << '(';
    Traceback(0,n-1);
    cout << ')';
    cout << endl;
}

函数MChain包含三重循环，循环体内的计算量为0(1)，所以算法的时间复杂度为0(n3)，空间复杂度为0(n2)。相对于穷举法的指数时间复杂度而言，动态规划算法更加有效。

备忘录方法

备忘录方法是动态规划法的一个变种。它采用分治法思想，以自顶向下直接递归的方式求最优解，但与分治法不同的是，备忘录方法为每一个已经计算的子问题建立备忘录，即保存子问题的计算结果以备需要时使用，从而避免了相同子问题的求解。

//矩阵连乘的备忘录方法
int MatrixChain::LookupChain(int i,int j)
{
    if(m[i][j]>0)
        return m[i][j];//子问题已经求解，直接引用
    if(i==j)
        return 0;//单一矩阵无需计算
    int u=LookupChain(i+1,j)+p[i]*p[k+1]*p[j+1];//求最小值
    s[i][j]=i;
    for(int k=i+1;k<j;k++){
        int t=LookupChain(i,k)+LookupChain(k+1,j)+p[i]*p[k+1]*p[j+1];
        if(t<u){
            u=t;
            s[i][j]=k;
        }
    }
    //保存并返回子最优解值
    m[i][j]=u;
    return u;
}
 
int MatrixChain::LookupChain()
{
    return LookupChain(0,n-1);//返回A[0:n-1]的最优解值
}

备忘录方法的时间复杂度也是O(n3)（因为总共有O(n2)个m[i][j]需要计算，这些元素的初始化需要O(n2)时间，计算一个元素的时间为O(n)）。

最长公共子序列

问题描述

定义：若给定序列X=(x1,x2,...,xm)，则另一序列Z=(z1,z2,...,zk)为X的子序列（subsequence）是指存在一个严格递增下标序列(i1,i2,..,ik)使得对于所有j=1,..,k有zj=xij。

定义：给定两个序列X和Y，当另一个序列Z既是X的子序列又是Y的子序列时，称Z是序列X和Y的公共子序列（common subsequence）。

最长公共子序列的穷举法：对于长度为m的序列X，其每个子序列都对应于下标集{1,2,…,m}的一个子集，X的子序列数多达2m，算法的时间是指数级的。

动态规划法求解

最优子结构

定理：设X=(x1,x2,⋯,xm)和Y=(y1,y2,⋯,yn)为两个序列，Z=(z1,z2,⋯,zk)是他们的最长公共子序列，则
+ 若xm=yn，则zk=xm=yn，且Zk−1是Xm−1和Yn−1的最长公共子序列；
+ 若xm≠yn，且zk≠xm，则Z是Xm−1和Y的最长公共子序列；
+ 若xm≠yn，且zk≠yn，则Z是X和Yn−1的最长公共子序列；

以上定理表明两个序列的最长公共子序列包含了这两个序列的前缀的最长公共子序列，这意味着最长公共子序列具有最优子结构特性。

最优解的递推关系

设有序列X=(x1,x2,⋯,xm)和Y=(y1,y2,⋯,yn)，根据最优解具有最优子结构特性，由此可以推导出以下递推关系：
+ 若xm=yn，则先求Xm−1和Yn−1的最长公共子序列，并在其尾部加上xm便得到Xm和Ym的最长公共子序列；
+ 若xm≠yn，则必须分别求解两个子问题Xm−1和Yn，以及Xm和Yn−1的最长公共子序列，这两个公共子序列中的较长者就是Xm和Ym的最长公共子序列；

与矩阵连乘类似，需要使用一个二维数组来保存最长公共子序列的长度，设c[i][j]保存Xi=(x1,x2,⋯,xi)和Yj=(y1,y2,⋯,yj)的最长公共子序列的长度。对于c[i][j]有以下递推式：

c[i][j]=⎧⎩⎨⎪⎪⎪⎪0c[i−1][j−1]+1max{c[i][j−1],c[i−1][j]}i=0或j=0i,j>0,xi=yji,j>0,xi≠yj

最长公共子序列算法

如果直接根据上面的递推式写一个计算c[i][j]的递归算法，但会得到一个指数时间的算法。采用动态规划法可以避免重复计算子问题，在本问题中的不同子问题的数目总计为Θ(m×n)，采用动态规划法进行自底向下求解，可在多项式时间内完成计算。由于每一个数组元素的计算时间为O(1)，则程序的时间复杂度为O(m×n)。

//动态规划法求LCS长度
class LCS
{
public:
    //创建二维数组c,s和一维数组a,b，并进行初始化
    LCS(int nx,int ny, char *x,char *y);
    //求最优解值（最长公共子序列长度）
    void LCSLength();
    //构造最优解（最长公共子序列）
    void CLCS();
    ...
private:
    void CLCS(int i,int j);
    int **c,**s,m,n;
    char *s,*b;
};
 
int LCS::LCSLength()
{
    for(int i = 1;i<=m;i++)
        c[i][0]=0;
    for(int i=1;i<n;i++)
        c[0][i]=0;
    for(int i = 1;i<=m;i++){
        for(int j=1;j<n;j++){
            if(x[i]==y[j]){
                c[i][j]=c[i-1][j-1]+1;
                s[i][j]=1;
                //由c[i-1][j-1]计算c[i][j]
            } else if(c[i-1][j]>=c[i][j-1]){
                c[i][j]=c[i-1][j];
                s[i][j]=2;
                //由从c[i-1][j]得到c[i][j]
            } else{
                c[i][j]=c[i][j-1];
                s[i][j]=3;
                //由c[i][j-1]得到c[i][j]
            }
        }
    }
    return c[m][n];//返回最优解值
}

设由序列X=(x1,x2,⋯,xm)和Y=(y1,y2,⋯,yn)，上述程序可以构造他们的最长公共子序列。从s[m][n]开始，如果s[i][j]=1，表示它是由Xi−1和Yj−1的最长公共子序列的尾部加上xi形成的；如果s[i][j]=2，表示它与Xi−1和Yj的最长公共子序列相同；如果s[i][j]=3，表示它与Xi和Yj−1的最长公共子序列相同；如果i=0或j=0，则为空子序列。

//构造最长公共子序列
void LCS::CLCS(int i,int j)
{
    if(i==0||j==0)
        return ;
    if(s[i][j]==1){
        CLCS(i-1,j-1);
        cout<<a[i];
    } else if(s[i][j]==2){
        CLCS(i-1,j);
    } else {
        CLCS(i,j-1);
    }
}

算法的改进

求LCS长度的程序中数组s是可以省去的，因为c[i][j]由c[i][j]=c[i−1][j−1]+1、c[i][j]=c[i−1][j]或c[i][j]=c[i][j−1]计算得来，因此为了确定c[i][j]是从者三者中哪一个计算得来的，可以直接由c的值确定。因此可以写一个类似的CLCS算法在O(m+n)时间内构造出最长公共子序列，该算法使用c而不是用s以减少存储空间。
另外，在只求最长公共子序列的长度、无需构造最优解时，也可以只用两行元素空间，时间是O(min{m,n})。