矩阵链乘法的动态规划算法_1

矩阵链乘法的动态规划算法(Matrix Chain Multiplication)

上次我们学习了rod-cutting 问题的动态规划算法，初步了解求解动态规划过程的CRCC步骤，此步骤对于可以运用动态优化的问题非常有用，类似给大家提供了一套思维模板，让我们能更系统的思考和解决问题。本次我们将讨论矩阵链乘法的动态规划算法。

矩阵乘法

在讨论矩阵链乘法之前，我们先了解一下矩阵乘法的基本要素，对于矩阵乘法C(i,k)=A(i,j) * B(j,k), 如果两个矩阵可以相乘，那么A矩阵的列数(j)必须等于B矩阵的行数(j)，否则我们成为两个矩阵乘法不相容，相乘的结果为矩阵C, 它的行等于A矩阵的行数，列数等于B矩阵的列数，对于两个矩阵的乘法，我们一般采用标量相乘然后求和即可。对于A*B相乘，给出一般的算法如下，我们注意到，标量相乘的m3->data[i] [j]+=(m1->data[i] [k]) * (m2->data[k] [j])的执行次数对程序的效率影响非常大，对于两个矩阵相乘，其循环的次数已经固定，无法优化。

typedef struct matrix
{
    int row;//row 矩阵行数
    int col;//col 矩阵列数
    int **data; //矩阵实际储存数据，本次探讨的为二维矩阵
}matrix, *matrix_ptr;

/**
	m3=m1*m3, 相乘后返回 指向矩阵 的指针，使矩阵运算统一化
*/
matrix_ptr matrix_multiply(matrix_ptr m1, matrix_ptr m2)
{
    matrix_ptr m3;
    int i;
    int j;
    int k;

    m3=(matrix_ptr)malloc(sizeof(matrix)*1);
    m3->row=m1->row;
    m3->col=m2->col;

    m3->data=(int **)malloc(sizeof(int *)*m3->row);

    for(i=0;i<m3->row;i++)
    {
        *(m3->data+i)=(int *)malloc(sizeof(int)*m3->col);
    }

    for(i=0;i<m1->row;i++)
    {
        for(j=0;j<m2->col;j++)
        {
            m3->data[i][j]=0;
            for(k=0;k<m1->col;k++)
            {
                m3->data[i][j]+=(m1->data[i][k])*(m2->data[k][j]);
            }
        }
    }

    return m3;
}

矩阵链乘法

那么矩阵链的乘法如何定义呢，我们假定存在矩阵<A1, A2, A3…An>，其之间相乘互容，那么此时标量相乘的次数就可以随着不同矩阵之间的结合而发生变化，列举简单示例来阐述这个观点，为了方便，我们假定有3个矩阵相乘(A1,A2,A3)，假定A1的行数=6，列数=7； A2的行数=7，列数=3; A3的行数=3，列数=1。

由于矩阵之间相互相容，此时就有两种途径来解决问题：

• C=(A1 * A2) * A3, 标量相乘的次数=6 * 7 * 3 + 6 * 3 * 1=144
• C=(A1 * (A2 * A3))， 标量相乘的次数=6 * 7 * 1 + 7 * 3 * 1 = 63

我们很容易发现，第二种方法的标量相乘次数仅有第一种算法的一半不到，示例中仅有三个矩阵相乘，如果矩阵数量增加，找到合适的矩阵相结合可以明显降低计算的次数，这将非常显著地提升矩阵链乘法的效率。

矩阵链乘法的动态规划

本文主要探讨如何按照合适的次序对矩阵进行相乘，取得最小的标量相乘次数，需要注意的是，我们并不需要实际对矩阵做乘法运算，仅仅求解矩阵链相乘某个最优次序，使得标量相乘的次数最少。首先确认此问题属于最优化问题的范畴，这有再这个前提之下，我们才能着手开展动态规划算法。

我们接下来按照动态规划的CRCC步骤依次对此问题进行剖析，最终求得动态规划的各个条件。

a) Characterize the structure of the optimal solution, 表征最优问题的基本结构

我们要找到基本结构，最重要的事情是能找到不同层次的子结构，并且其求解具有相似性，这样才能满足后面的递归或者迭代的计算条件。让我们假定矩阵链的下标可以用数组p来表示p[]=<p0,p1,p2,p3…pn>，此数组中包含了A=<A1,A2,A3…An>所有的下标的集合。为了方便阐述问题，本文的研究对象就转化为求解p数组内部标量运算的次数的最小值的某一序列（可以用括号划分和嵌套，表示矩阵之间乘法运算的优先级，从而使标量运算的总次数最小）。

假定我们把矩阵链<A1,A2,A3…An>划分为两个子问题，<A1,A2,A3…Ak>和<Ak+1,Ak+2,Ak+3…An>，我们就会发现第一个子问题和第二个子问题的表示形式存在差异，因为第一个子问题的前段是固定的，永远都是A1；与此同时，第二个子问题的后段是固定的，永远都为An，显然这两个子结构不相同。

接下来我们继续尝试，把矩阵链<A1,A2,A3…An>划分为<A1,A2,A3…Ak>和<Ak+1,Ak+2,Ak+3…Aj>两个子问题，我们就会发现两个子问题仍然不同，第一个子问题的左端仍然固定，第二个子问题两端可以进行变化，如果我们总能保证k=j-1,那么问题也可以划分成这类形式，因为第二个问题的标量运算总是0,因为第二个子问题会退化为< Aj>的单个矩阵的形式，从而可以直接递归返回，再往下看，实际上就退化为一个子问题，显然这与实际情况不相符合。

最后，通过前述的尝试，我们如果把矩阵链划分为<Ai,Ai+1…Ak>和<Ak+1,Ak+2,…Aj>两个子问题，那么问题就迎刃而解，此问题可以描述为，如果要求解<Ai,Ai+1,…Ak-1,Ak>的标量运算的总数的最小值，我们需要把此问题分解为两个子问题<Ai,Ai+1…Ak>和<Ak+1,Ak+2,…Aj>。如果<Ai,Ai+1,…Ak-1,Ak>取得最优解，那么其两个子问题一定是最优解，也即<Ai,Ai+1…Ak>和<Ak+1,Ak+2,…Aj>此时也是最优解。此时就具备了递归的雏形，可以再接下来一步总来定义递归表达式的值得求解。

可能我们会有疑惑，为什么在<Ai,Ai+1,…Ak-1,Ak>获得标量运算次数最小的时候，它的两个子问题也同时具备运算次数最小呢？ 其实这个问题就回归到 原问题和子问题的逻辑关系，可以用反正法来证明，如果子问题不具备最优的解，那么一定在这个子问题中找到最优解，那么此时原问题的最优解就不是最优解，这与事实相矛盾。

b) Recursively define the value of the optimal solution-递归定义最优解的值

这一步是后续程序实现的关键，如果能顺利递归定义表达式，那么后面的递归程序就是水到渠成的事情了。让我们继续回到问题的本身，为了方便问题描述，我们引入 m[i,j]数组，此数组表示矩阵链A i:j 标量运算次数的最小值，引入m[i,j]数组以后，我们就会发现当i=j的时候，实际上是递归的m[i,j]=0， 这个结论非常重要，它其实肩负着递归的终止点和迭代的起始点的角色，没有i=j的判断，我们的递归将无限循环，同时如果采用迭代也无法定义起点。我们接着分析i<j的情况，这时候我们可以使用递归公式：

m[i,j]=m[i,k]+m[k+1,j]+p[i-1] * p[k] * p[j];

如果我们要找到某个k值，使其m[i,j]最小，综合上述的分析，此表达式可以进一步表述为：

这时候你可以选择使用普通递归（接近暴力算法），top-down形式的memo递归 或者 bottom-up形式的迭代，都可以对此问题作出合理的解答。

c) Compute the value of the optimal solution - 计算最优值

我们上述完成递归定义的最优解，也有人称之为状态转移方程，意味着问题从某个子状态（底层）转移到父状态（上一层），这个定义是求解所有动态规划问题核心中的核心，如果无法完成上述定义，那么后续的计算最优解和构造解的详细信息就无法完成。

为了获得矩阵链的分割位置，我们可以定义式split数组s[i,j]来记录过程中的k值，k值便是Ai:Aj范围内的最优化的分割点。

c-1) 迭代(bottom-up)计算方法的伪代码

在bottom-up问题计算过程中，我们按照长度从小到大进行排列，也就是先从“小”的子问题开始计算，然后再计算较大的子问题，确保每次迭代之前，计算某个问题的前置子问题的解已经求出，从而可以计算当前的某个问题。所以l从小到大求值非常关键。另外就规模和长度而言，子问题都要小于父问题，比如m[i,k]和m[k+1,j]两个子问题长度都小于m[i,j]，这就保证了问题求解的前置性可以得到很好满足。

c-2) 递归top-down（带Memo)计算方法的伪代码

在这个过程当中，我们特别需要注意的是递归的终点，递归的终点为i=j, 同时可以看到，除了递归本身内部的循环，找到合适的k值外，里面还有两个递归函数，每个递归函数的复杂度为O(n)，两个递归函数嵌套调用，其复杂度为O(n2)，如果再考虑k=i to j-1的复杂度为O(n)，那么整个递归运算的复杂度就为O(n3)，与前述的迭代运算复杂度相同。

d) Construct the the optimal solution from computed information 构建最优解的具体信息

很多情况下，动态规划要解决的问题时候都涉及到相应的选择，本例中就涉及到矩阵的划分位置，我们在递归或迭代过程中采用了s[i,j]数组进行了记录，根据s[i,j]的基本信息，我们可以通过循环或者递归构建出最优解的具体方案。就本问题而言，我们构建最优解的伪代码表述为。

总结

我们采用了《算法导论》介绍的CRCC步骤对矩阵链的乘法运算的最小标量运算次数进行了分析，这个过程涉及到两个子问题的求解过程。

通过学习，进一步理解了更复杂动态规划的过程以及求解的具体步骤。

参考资料

1. 《Introduction to Algorithm, 4ed, edition》