Algorithm Experiment：动态规划算法——矩阵连乘问题【Dynamic programming algorithm -- matrix multiplication problem】

一、算法的基本思想和应用要点

动态规划与分治法相似，都是采用将大问题分成小问题，然后组合小问题的解德奥大问题的解的方法。不同的是，分治法中的小问题之间是相互独立的，而动态规划中的小问题之间是重叠的。即：在子问题重叠的情况下，如果使用分治法就会在递归的过程中重复的执行某项工作（公共的子问题）。而动态规划算法只会对这项工作求解一次并保存下来，在下一次需要的时候直接使用这个结果，从而得到性能优于分治法的策略。
在应用该算法之前，首先应明确适合应用动态规划算法求解的最优化问题应具备的两个要素：

（1）最优子结构，即问题的最优解应包含其子问题的最优解。
（2）重叠子问题，即“问题的递归算法会反复地求解相同的子问题，而不是一直生成新的子问题”【1】。

在应用动态规划算法解最优化问题时，应按照以下步骤设计算法【1】：

刻画一个最优解的结构特征。
递归地定义最优解的值。
计算最优解的值，通常采用自底向上的方法。
利用计算出的信息构造一个最优解。

二、 问题描述

给定n个矩阵{A¬1,A2,…,An}，其中Ai与Ai+1是可乘的，i=1,2,…,n-1。确定计算矩阵连乘积的计算次序，使得到正确结果需要付出的代价最少。

输入：矩阵个数，每个矩阵的规模
输出：计算次序

例如：以{A1，A2，A3}为例，假设矩阵规模分别为23，34,4*5共有两种不同的计算次序：

((A1A2)A3)：需要进行2*3*4+2*4*5=64次
(A1(A2A3））：需要进行3*4*5+2*3*5=90次

第一种计算顺序的速度明显优于第二种。

三、理论分析

由于矩阵乘法满足结合律，所以无论哪种加括号的方式得出的结果都是一样的。所以计算次序可以通过不同的加括号方式来确定。
如果使用穷举法，对于n个矩阵相乘，设P(n)是计算次序的个数，则：

解此递归方程可得这个序列的增长速度是Ω（4^n/n(3/2) ）,是呈几何增长的，所以穷举法对这个问题不是个有效算法。
如果采用动态规划方法：
分析最优解的结构：
用A[i,j]表示AiAi+1…Aj乘积的结果矩阵，在对AiAi+1…Aj进行加括号时就相当于在某个位置k（i≤k≤j）处将矩阵链划分成两个部分。首先计算AiAi+1…Ak，在计算Ak…Aj，再将二者的结果在进行乘积运算得到最终结果。计算代价有三部分构成：A[i,k]的计算代价、A[k,j]的计算代价、二者相乘的计算代价。其中，在第一次得出k的划分后，在计算A[i,k]时，可继续进行寻找最优解k‘的操作，后半部分亦是如此。最终用二者的最优解方案得出的结果相乘得到总的最优解。
建立递归方案：用m[i][j]表示最少乘次数
①i=j时，m[i][j]=0；
②i<j时，m[i][j]=m[I,k]+m[k+1,j]+pi-1pkpj，由于在计算是并不知道断开点k的位置，所以k还未定。不过k的位置只有j-i个可能。因此，k是这j-i个位置使计算量达到最小的那个位置。(此处的k是假定k已知，实际上我们是不知道的，需要通过穷举得到)
综上，有如下递推关系：

这个式子给出了最优解的代价，但没有给存储最优位置分割点k的方案，所以在编程的时候可以使用s[i,j]来保存k的位置。在计算出最优值m[i][j]后可递归地由s[i][j]构造出相应的最优解。

四、算法实现

4.1代码实现

#include<stdio.h>
#include <iostream>
using namespace std;

void MatrixChain(int n, int p[], int m[][100], int s[][100])
{
    for (int i = 1; i <= n; i++) 
        m[i][i] = 0;
  
    for (int r = 2; r <= n; r++) {
        for (int i = 1; i <= n - r + 1; i++) {//第i行
            int j = i + r - 1;
            m[i][j] = m[i + 1][j] + p[i - 1] * p[i] * p[j];//在i处取得最优解
            s[i][j] = i;//断开位置
            for (int k = i + 1; k < j; k++) {//如果有更小的则被替换
                int t = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];//从k处断开
                if (t < m[i][j])
                {
                    m[i][j] = t;
                    s[i][j] = k;//断点
                }
            }
        }
    }
}
void printConclusions(int s[100][100], int i, int j)//输出结果
{
    if (i == j)
        cout << "A" << i;
    else
    {
        cout << "(";
        printConclusions(s, i, s[i][j]);
        printConclusions(s, s[i][j] + 1, j);
        cout << ")";
    }

}
int main()
{
    int p[100];//存储第一个矩阵的行数以及所有矩阵的列数
    int m[100][100];//存储最优值
    int s[100][100];//存储断开位置
    memset(p, 0, sizeof(p));//初始化p[]
    memset(m, 0, sizeof(m));//初始化m[][]
    memset(s, 0, sizeof(s));//初始化s[][]
    cout << "请输入矩阵的个数" << endl;
    int n;
    cin >> n;
    cout << "请依次输入第一个矩阵的行数以及所有矩阵的列数" << endl;
    for (int i = 0; i <= n; i++) {
        cin >> p[i];
    }
    MatrixChain(n, p, m, s);
    cout << "最少的相乘次数是:" << m[1][n] << endl;

    cout << "结果为:" << endl;
    printConclusions(s, 1, n);
    return 0;
}

运行结果：

4.2关键代码说明

void MatrixChain(int n, int p[], int m[][100], int s[][100])
//n存储矩阵个数
//p[]存储矩阵的列数，p[0]存储的是第一个数组的行数
//m[][]存储的是最优值
//s[][]存储的是最优值的断开点的位置

{
    for (int i = 1; i <= n; i++) 
        m[i][i] = 0;//从0开始
  
    for (int r = 2; r <= n; r++) {
        for (int i = 1; i <= n - r + 1; i++) {//第i行
            int j = i + r - 1;//r是指Ai和Aj之间的长度，即A的个数
            m[i][j] = m[i + 1][j] + p[i - 1] * p[i] * p[j];//在i处取得最优解
            s[i][j] = i;//断开位置
            for (int k = i + 1; k < j; k++) {//如果有更小的则被替换
                int t = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];//从k处断开
                if (t < m[i][j])
                {
                    m[i][j] = t;
                    s[i][j] = k;//断点位置
                }
            }
        }
    }

实验总结

balabalabala…

参考书籍

[1] [美] Thomas H.Cormen，[美] Charles E.Leiserson，[美] Ronald L.Rivest，[美] Clifford Stein 著，殷建平，徐云，王刚 等 译. 机械工业出版社, 2012. ISBN: 97871114070