Preface

本人最近在做强化学习的内容，我发现强化学习基础当中马尔科夫决策过程（MDP）的求解（策略迭代，值迭代）与算法基础当中的动态规划密切相关。但由于本人在本科阶段没有认真看算法导论，所以现在在研究生阶段就补补，并为强化学习打好基础。

Overview

动态规划是将求最优问题分解为求最优子问题，而各子问题包含着公共子问题。所以动态规划常常适用于有重叠子问题和最优子结构性质的问题。可能大家比较熟悉的一个例子就是斐波那契数列。它正正是满足了这些条件重叠子问题和最优子结构性质。

• 重叠子问题：一个问题的自子问题可能与其他问题的自问题是一样的。虽然看起来要求解问题很多，但是其他很多问题都是重复的，我们只需要计算一次就可以了。这就显得子问题的个数“很少”。
• 最优子结构性质：问题能够分解成子问题来解决。（递归求解子问题）

矩阵链乘法

我们用动态规划来解决矩阵链相乘的问题。给定由n个要相乘的矩阵构成的序列链

公式-1

可将两个矩阵相乘的标准算法作为一个子程序，根据括号的顺序依次执行矩阵乘法。根据矩阵乘法的结合率，无论我们怎么样给他加括号矩阵乘法得到的结果都不会有影响。

括号乘法-1
括号乘法-2

上面两个乘法得到结果都是一致的，但是他们计算代价却不一定一样。我们在这里先说明一下两个矩阵相乘，如A和B，他们的大小为p×q和q×r。

void arymul(martix A, martix B)  
{  
    int i, j, k;
    martix C;
    for(i = 0; i < rows(A); i++){  
        for(j = 0; j < columns(B); j++){  
            c[i][j] = 0;  
            for(k = 0; k < columns(A); k++){  
                c[i][j] += A[i][k] * B[k][j];  
            }
        }
    }  
}  

可见矩阵乘法的代价为p×q×r。
假设我们有两个矩阵A，B和C，它们分别的大小为10×100，100×5和5×50。如果按照(AB)C的顺序来计算矩阵乘法，它的计算代价就是10 × 100 × 5 + 10 × 5 × 50 = 7500。如果按照A(BC)的顺序来计算矩阵乘法，它的计算代价就是100 × 5 × 50 + 10 × 100 × 50 = 75000。计算量相差10倍，可见优化矩阵乘法的顺序有助于提高算法性能。现在我们就提出一个问题：给定一个矩阵乘法链，求最优的求解顺序。这利用到动态规划了。

Analysis

• Step1：寻找最优的子结构

现在我们要对矩阵乘法链进行优化计算顺序（如图1），对于乘积

公式-2

我们需要寻找其中一个矩阵位于第i个矩阵和第j个矩阵之间，假设为第k个矩阵。

公式-3

其使得公式-2的计算最优，这样我们必须保证：

公式-4
公式-5

公式-4和公式-5的计算顺序为最优子顺序。若存在更优子顺序，那么我们必然可以用更优的子顺序代替公式-4和公式-5，也就是k并非最优的分割位置。这里有点类是骨牌原理，但这里是往后递归推导。所以我们要寻找子问题实例的最优解，然后合并这些子问题的最优解，来构建一个矩阵链乘法问题实例的最优解。必须寻找一个正确的位置分割乘积。

• Step2：递归公式
  我们寻找一个最优（代价最小）的位置将矩阵分割。这就需要设计一个递归公式：

公式-6

其中m[i, j]代表A_{i}到A_{j}矩阵链乘积的最小代价。有了这个递归公式我们就可以很好的写出求矩阵链乘积最优顺序解的程序来了。

代码实现

假设现在有一下6个矩阵：

# -*-coding:utf-8-*-#
import numpy as np

matrix_shapes = np.array([[30, 35], [35, 15], [15, 5], [5, 10], [10, 20], [20, 25]])

length = len(matrix_shapes)
m = np.full((length, length), np.inf)
s = np.zeros((length, length), dtype=np.int32)

def maxtrix_chain_order(shapes, start, end):
    if start == end:
        m[start, end] = 0
    else:
        q = 0
        matrix_start = shapes[start]
        matrix_end = shapes[end]
        q = 0
        for i in range(start, end):
            matrix_k = shapes[i]
            q = maxtrix_chain_order(shapes, start, i) + maxtrix_chain_order(shapes, i + 1, end) + matrix_start[0] * matrix_k[1] * matrix_end[1]
            if m[start, end] > q:
                m[start, end] = q
                s[start, end] = i + 1
    return m[start, end]

maxtrix_chain_order(matrix_shapes, 0, length - 1)
print(m)
print(s)

def print_optimal_parens(s, i, j):
    if i == j:
        print("A_{%d}" % (i + 1), end="")
    else:
        print("(", end="")
        print_optimal_parens(s, i, s[i, j] - 1)
        print_optimal_parens(s, s[i, j], j)
        print(")", end="")

print_optimal_parens(s, 0, length - 1)
print()

上图的上三角矩阵视图：

m矩阵和s矩阵

最后打印出来的结果：

打印结果

Conclusion

动态规划的基本步骤：

1. 描述最优解的结构
2. 递归定义最优解的值（重要）
3. 自底往上的方式计算最优解的值
4. 由计算出的结果构建一个最优解

但是我们必须要确保每个子问题之间是独立的不存在依赖关系（无权最长简单路径的每个子问题就不是独立），但每个问题的子问题可以是重叠的。像上图m矩阵那样子，越往底下（越靠近对角线）的值被调用的次数多，证明它是多个问题的自问，所以我们需要把它的值记录下来，重复运算相同的值是很浪费资源的。

好啦，今天就到这里了。早唞！好梦！

References

• 算法导论
• http://blog.csdn.net/tengweitw/article/details/16922431