动态规划之矩阵链相乘

在看这张之前，最好看看我写的动态规划详解，里面都是讲理论基础，我下面的分析都是在此基础上进展的。

问题描述（详见算法导论P197-P198）：

已知：给定n个矩阵构成的一个矩阵链(A1, A2, ..., An)，矩阵Ai的维数为pi-1×pi

求：决定该矩阵链的乘法结合顺序（即加括号），使得矩阵链乘法的运行时间最短

几个前提概念

• 矩阵链乘法的运行时间将以标量乘法（单行×单列）的次数来衡量
• A是p×q矩阵，B是q×r矩阵，则A×B的运行时间为pqr
• 矩阵乘法满足结合律

使用动态规划的四个步骤来分析：
①.描述最优解的结构。

书中直接就说利用Ai...Aj的乘积矩阵来描述，跳跃性有点大，不仔细分析，还不知道里面的水有多深。那现在来探讨一下。

还是按照正常思路，从总数n开始划分，总的矩阵相乘代价最小，两个大矩阵的划分有n-1种可能，因此可以任意选择一个i值来分析（有的问题只能由n到n-1，有一定限制，所以第一步划分时要警惕一下）。

选择一个划分结点i，分成两个矩阵相乘，A1...Ai和Ai+1...An。两个子矩阵都不能被对方代替，如果不好看出递归方程，可以再一次划分，将A1...Ai划分成A1...Aj和Aj+1...Ai。从这里看出，子矩阵的左右范围是变化的，不能通过以为数组来描述，这时就考虑用二维数组（一般备份存储空间都是数组来表示的，一维二维居多，可尝试）。

②.递归定义最优解的值。

既然用二维数组，则m[i,j]表示Ai...Aj的乘积。当选取中间点k来划分时，m[i,j] = m[i,k]+m[k+1,j]+分割代价。Ai...Ak的乘积为行数和列数为pi-1和pk,那么划分后子数组相乘的代价为pi-1 * pk *pj。所以m[i,j] = m[i,k]+m[k+1,j]+pi-1 * pk *pj。

但要对下标变量有限制，表达式中i<=k<j，保证不取到右即可，对k的可能性取值进行筛选，最小值即为m[i,j]的取值。考察临界点，当i=j时，m[i,i]=0。整体表达式为

③.按自底而上的方式计算最优解的值。

现在的问题是找到备份数组，在第一步的分析中已经阐述了这个问题。当递归表达式中数组m换成函数时，就是一个暴力求解的递归算法。这里设置二维数组m[1...n,1...n]来存储数据。

但是如何实现自底向上的递归？由于m[i,j]依赖于m[i,k]和m[k+1,j],而k是一直变化的。依赖于长度小于j-i+1的数据。最小是两个相邻的矩阵相乘。因此，这里以长度为主变量来实现自底向上：L←2-n。

④.由计算出的结果创造一个最优解。

通常重构最优解信息保存的数组与备份数组一样，此处为s[n,n]。s[i,j]保存m[i,j] = m[i,k]+m[k+1,j]+pi-1 * pk *pj表达式中的分割点k的取值。这个递归方程很典型，方程中就有关分割点的信息，当更新m[i,j]的时候，同步更新s[i,j]的值。

伪代码：

MATRIX-CHAIN-ORDER(p)
{
    n = p.length-1; //因为p从0到n
    let m[1...n,1...n]and s[1...n-1,2...n]be new tables //m和s的维数一样，用到的参数范围可能不同！要注意，可能s需要的范围大
    for(i=1;i<=n;i++)   //根据方程初始化一些参数，为了防止错误，最后全部初始化
        m[i,i]=0;
    for(L=2;L<=n;L++){  //以长度变化来遍历
        for(i=1;i<=n-L+1;i++){
            j=i+L-1;
            m[i,j] = 100000;
            for(k=i;k<j;k++){   //k的范围也是根据方程表达式来写出
                q = m[i,k] + m[k+1,j] + pi-1 * pk * pj;
                if(q < m[i,j]){
                    m[i,j] = q;
                    s[i,j] = k;
                }
            }
        }
    }
    
    return m and s;
}

很巧妙的方法：上面的分析说，利用以长度为主变量遍历，但网上有一种方法只是用i,j,k来遍历，不过这里最外层遍历的是j。但不建议这种方法，太巧秒，容易出错。

    for(j=1;j<=N;j++){                  
        for (i=j;i>=1;i--)          //当i=j时,m[i][j]=0，                             
        {                           //当i<j时,m[i][j]=min{m[i][k]+m[k+1][j]+p(i-1)p(k)p(j)} i=<k<j                    
            if (j==i)  
                m[i][j]=0;  
            else  
            {  
                m[i][j]=600000;   
                for (k=i;k<j;k++)  
                {  
                    q=m[i][k]+m[k+1][j]+matrix[i-1]*matrix[k]*matrix[j];  
                    if (q<m[i][j])  
                    {  
                        m[i][j]=q;  
                        s[i][j]=k;  
                    }  
                }  
            }  
        } 
    }

输出分配方案

PRINT-OPTIMAL-PARENS(s,i,j)
{
    if(i==j)
        printf("Ai");
    else{
        printf("(");
        PRINT-OPTIMAL-PARENS(s,i,s[i,j]);
        PRINT-OPTIMAL-PARENS(s,s[i,j]+1,j);
        printf(")");
    }
}