算法分析与设计：矩阵连乘问题（动态规划）

问题描述

给定$n$个矩阵$A1，A2，\dots \dots ，An$,其中$A_i$与$A_{i+1}$是可乘的，$i=1，2，\dots \dots ，n-1$。

例如：
计算三个矩阵连乘$A_1,A_2,A_3$；维数分别为$10\times 100$ , $100\times 5$ , $5\times 50$

按此顺序计算需要的次数$((A_1\ast A_2)\ast A_3)$：$10\times 100\times 5+10\times 5\times 50=7500$次
按此顺序计算需要的次数$(A_1\ast (A_2\ast A_3))$： $10\times 5\times 50+10\times 100\times 50=75000$次

所以要解决的问题是：如何确定矩阵连乘积$A_1,A_2,\dots \dots A_n$的计算次序，使得按此计算次序计算矩阵连乘积需要的数乘次数达到最小化。

---

动态规划

将待求解问题分解为若干子问题，通过子问题的解得到原问题的解，这是问题求解的有效途径。但是如何实施分解呢？分治策略的基本思想是将规模为$n$的问题分解$k$规模较小的子问题，各子问题相互独立但与原问题求解策略相同。但并不是所有问题都可以这样处理。

问题分解的另一个途径是将求解过程分解为若干阶段（级），一次求解每个阶段即得到原问题的解。通过分解得到的各子阶段不要求相互独立，但希望它们具有相同类型，而且前一阶段的输出可以作为下一阶段的输入。这种策略特别适合求解具有某种最优性质的问题。贪心法属于这类策略：对问题$P(n)$，其求解过程中各贪心选择步骤构成决策序列$D=D_1,D_2,\dots ,D_k$， 的最优性仅仅依赖于$D_1,D_2,\dots ,D_{(i-1)}$。贪心法不保证最后求出解的最优性。

动态规划法也是一个分阶段判定决策过程，其问题求解策略的基础是决策过程的最优原理：为达到某问题的最优目标$T$，需要一次作出决策序列$D=D_1,D_2,\dots ,D_k$。如果$T$是最优的，则对任意$i(1\le i\le k)$，决策子序列$D_{(i+1)},D_{(i+2)},\dots ,D_k$也是最优的，即当前决策的最优性取决于其后续决策序列是否最优。由此追溯至目标，再由最终目标决策向上回溯，导出决策序列$D=D_1,D_2,\dots ,D_k$。因此动态规划方法可以保证问题求解是全局最优的。

---

矩阵连乘问题的伪码算法

---

矩阵连乘问题C++源代码

#include <iostream>
#include <iomanip>
#include <malloc.h>
#include <string.h>
using namespace std;

/*****************************************************************
* 函数描述： 构造矩阵（n阶方阵）
* 函数参数： n——矩阵行数
			m——矩阵列数
* 函数返回： 构造的矩阵指针
*****************************************************************/
int **GenerateMatrix(int n, int m)
{
    int **s = (int **)malloc(sizeof(int *) * n);
    for (int i = 0; i < n; ++i)
    {
        s[i] = (int *)malloc(sizeof(int) * m);
        memset(s[i], 0, sizeof(int) * m);
    }
    return s;
}

/*****************************************************************
* 函数描述： 释放内存（n阶方阵）
* 函数参数： p——需要被释放内存的矩阵
			n——阶数
* 函数返回： void
*****************************************************************/
void FreeMatrix(int **p, int n)
{
    for (int i = 0; i < n; ++i)
    {
        free(p[i]);
    }
    free(p);
}

/*****************************************************************
* 函数描述： 打印矩阵（不打印0索引）
* 函数参数： p——被打印矩阵
			n——行数
			m——列数
* 函数返回： void
*****************************************************************/
void PrintMatrix(int **p, int n, int m)
{
    for (int i = 1; i < n; ++i)
    {
        for (int j = 1; j < m; ++j)
        {
            cout << setw(7) << p[i][j];
        }
        cout << endl;
    }
    cout << endl;
}

/*****************************************************************
* 函数描述： 矩阵乘法 a*b=c
* 函数参数： a——矩阵a
			b——矩阵b
			c——结果矩阵c
			ra——矩阵a的行数
			ca——矩阵a的列数
			rb——矩阵b的行数
            ca——矩阵b的列数
* 函数返回： void
*****************************************************************/
void MatrixMultiply(int **a, int **b, int **c, int ra, int ca, int rb, int cb)
{
    if (ca != rb) // 如果不满足矩阵乘法规则（左列数=右行数），则返回
    {
        cout << "Matrix multiplication error." << endl;
        return;
    }
    // 矩阵乘法
    for (int i = 0; i < ra; ++i)
    {
        for (int j = 0; j < cb; ++j)
        {
            int tmp = 0;
            for (int k = 0; k < ca; ++k)
            {
                tmp += a[i][k] * b[k][j];
            }
            c[i][j] = tmp;
        }
    }
}

/*****************************************************************
* 函数描述： 矩阵连乘核心算法
* 函数参数： p——连乘矩阵行列值数组
			n——连乘矩阵的个数
			m——储存分步最优值数组
			s——储存分步最优解，最优值断开的位置
* 函数返回： 最优值
*****************************************************************/
int MatrixChain(int *p, int n, int **m, int **s)
{
    int i, j;
    //对角线上的元素都设置成0
    for (i = 1; i <= n; ++i)
        m[i][i] = 0;
    for (int len = 2; len <= n; ++len)
    { // 遍历连乘长度，从2开始，子问题规模
        for (i = 1; i <= n - len + 1; ++i)
        {
            j = i + len - 1;
            // 将链ij划分为A(i) * ( A[i+1:j] )
            m[i][j] = 0 + m[i + 1][j] + p[i - 1] * p[i] * p[j];
            s[i][j] = i;
            for (int k = i + 1; k < j; ++k)
            { // 将链ij划分为( A[i:k] ) * ( A[k+1:j] )
                int t = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];
                if (m[i][j] > t)
                {
                    m[i][j] = t;
                    s[i][j] = k;
                }
            }
        }
        PrintMatrix(m, n + 1, n + 1);
    }
    return m[1][n]; // 返回最优值
}

/*****************************************************************
* 函数描述： 递归求最优解
* 函数参数： s——最优解断开位置矩阵
			i——连乘矩阵链起始位置
			j——连乘矩阵链终止位置
* 函数返回： void
*****************************************************************/
void Traceback(int **s, int i, int j)
{
    if (i == j)
        return;
    Traceback(s, i, s[i][j]);
    Traceback(s, s[i][j] + 1, j);
    cout << "Matrix A" << i << "," << s[i][j] << " and A" << s[i][j] + 1 << "," << j << endl;
}

int main()
{
    // 连乘的矩阵个数
    const int n = 6;
    // n个矩阵连乘，行列顺序。30*35 35*15 15*5 5*10 10*20 20*25
    int p[n + 1] = {30, 35, 15, 5, 10, 20, 25};
    int **m = GenerateMatrix(n + 1, n + 1);
    int **s = GenerateMatrix(n + 1, n + 1);
    int min = MatrixChain(p, n, m, s);
    cout << "Min: " << min << endl;
    PrintMatrix(m, n + 1, n + 1); // 打印储存最优值的矩阵
    PrintMatrix(s, n + 1, n + 1); // 打印储存最优解的矩阵
    Traceback(s, 1, n);           // 打印最优解
    return 0;
}

---

测试数据和执行结果

下图展示的了行列顺序为 $(30\times 35),(35\times 15),(15\times 5),(5\times 10),(10\times 20),(20\times 25)$ 的六个矩阵的连乘结果。由图可知，最小的运算次数为$15125$次。由最优解的记录矩阵可以看出首先在③号矩阵处断开，分成 $①-③$ 和 $④-⑥$ 两部分。$①-③$部分在①处断开，$④-⑥$ 部分在 ⑤ 处断开。

---

总结

矩阵连乘算法的关键在于理解最优值矩阵和最优解矩阵的产生和分析。产生主要是通过连乘长度（$len$）、子问题起始位置（$i$）、子问题分割点（$k$）这三重循环的遍历。每次进行比较后将局部最优的结果储存，同时将子问题的分割方法储存。在分析时，逆向倒退，自顶向下，通过分割点的记录得出最优解。

---

动态规划 VS 分治法

本文中用到的动态规划的方法与分治法十分类似，二者都要求原问题具有最优子结构性质，都是将原问题分而治之，分解成若干个规模较小很容易解决的子问题，之后将子问题的解合并，形成原问题的解。但是动态规划是将将分解后的子问题理解为相互间有联系，有重叠的部分，需要记忆储存，通常用迭代来做。而分治法是将分解后的子问题看成相互独立的子问题，通常使用递归来解决。