算法设计与分析：动态规划 - 矩阵链式相乘问题

朴灿烈づ我的快乐病毒、 2023-06-06 12:09 13阅读 0赞

### 文章目录 ###

*  最优化问题
 *   *  分治与动态规划
     *  多步决策过程
 *  矩阵链式相乘问题
 *   *  问题描述与分析
     *  解决思路
     *   *  初始版本
         *  改进版本
         *  算法运行过程剖析
         *  用迭代替代递归
 *  总结

> 本文参考UCAS卜东波老师算法设计与分析课程撰写

# 最优化问题 #

什么是最优化问题？考虑一个问题，包含一个目标函数，其输入是可行解（可行解由问题的约束规定），输出是可行解的质量刻画（可行解是否足够好），一般求最小或最大。

> 如一个简单的高中数学问题，已知 x 1 + x 2 + x 3 = 10 , x 1 , x 2 , x 3 ∈ N + , 求 g ( x 1 , x 2 , x 3 ) = x 1 x 2 x 3 的 m a x x\_1+x\_2+x\_3=10,x\_1,x\_2,x\_3\\in N^+,求g(x\_1,x\_2,x\_3)=x\_1x\_2x\_3的max x1\+x2\+x3=10,x1,x2,x3∈N\+,求g(x1,x2,x3)=x1x2x3的max，在这个问题中目标函数是 g g g,约束条件是为正整数与和为10，所有满足约束条件的 x 1 , x 2 , x 3 x\_1,x\_2,x\_3 x1,x2,x3都是可行解，但使得 g g g最大的可行解为最优解，最优化问题就是找这种最优解。

动态规划则十分适合用于处理这类问题

## 分治与动态规划 ##

基于以前的分治思想，如果将分治应用于最优化问题，与原始相比，有以下两点不同：

*  将原问题分解为相似子问题  
    在不知道如何分解才能获得最优的情况下，我们不能像以前一样简单地按下标、值或者区域划分，必须要**枚举所有可能的划分方案，通过比较获得最优**
 *  用子问题的解组成原问题的解  
    不像以往直接将子问题的可行解组合，这里需要通过**子问题的最优解组合得到原问题的最优解**

假设我们能够做到将子问题的最优解组合成原问题的最优解（最优子结构性质），那么就可以应用分治的思想将一个复杂问题的最优化转换成简单小问题的最优化，实现问题的解决。这种方法就是动态规划。

## 多步决策过程 ##

动态规划方法适用于可以描述为多步决策过程的最优化问题。多步决策要求完整解能够被一步步构建出来。

> 这里以著名的背包问题为例，给你一个容量为C的背包，现在有一些财宝可供你带走，这些财宝的价值 v v v和对应占用空间 c c c由一个集合给出 \{ ( v 1 , c 2 ) , ( v 2 , c 2 ) , . . . , ( v n , c n ) \} \\\{(v\_1,c\_2),(v\_2,c\_2),...,(v\_n,c\_n)\\\} \{ (v1,c2),(v2,c2),...,(vn,cn)\}，问你该如何选择，才能让获得的价值最大？

我们将每一次选择的过程看作一次决策，一共可以进行k次决策（直到背包容量无法再放下新的财宝）。这就是一个典型的多步决策过程。我们用**矩阵链式相乘问题**来引入动态规划的应用。

# 矩阵链式相乘问题 #

## 问题描述与分析 ##

*  **给定n个矩阵，求 A 1 A 2 . . . A n A\_1A\_2...A\_n A1A2...An，其中 A i A\_i Ai的规模是 p i − 1 p\_\{i-1\} pi−1行， p i p\_i pi列**  
    考虑最简单的做法，就是矩阵从左到右一个个乘过去。但实际上矩阵的乘法满足结合律，我们可以先乘前面的相邻矩阵，也可以先乘后面的相邻矩阵。下面是一个计算例子：  
     A 1 = \[ 1 2 \] A 2 = \[ 1 2 3 1 2 3 \] A 3 = \[ 1 2 3 4 1 2 3 4 1 2 3 4 \] A\_1 =\\left \[ \\begin\{matrix\} 1 & 2 \\end\{matrix\} \\right \] A\_2 = \\left \[ \\begin\{matrix\} 1 & 2 & 3 \\\\ 1 & 2 & 3 \\\\ \\end\{matrix\} \\right \] A\_3 = \\left \[ \\begin\{matrix\} 1 & 2 & 3 & 4 \\\\ 1 & 2 & 3 & 4\\\\ 1 & 2 & 3 & 4 \\end\{matrix\} \\right \] \\\\ A1=\[12\]A2=\[112233\]A3=⎣⎡111222333444⎦⎤  
    计算 A 1 A 2 A 3 A\_1A\_2A\_3 A1A2A3有两种方式， ( A 1 A 2 ) A 3 (A\_1A\_2)A\_3 (A1A2)A3和 A 1 ( A 2 A 3 ) A\_1(A\_2A\_3) A1(A2A3)，即先算前面两个和先算后面两个。那么计算的次数分别是多少呢？
 *   ( A 1 A 2 ) A 3 (A\_1A\_2)A\_3 (A1A2)A3  
     1 ( A 1 行 数 ) × 2 ( A 1 每 行 元 素 和 A 2 每 列 元 素 对 应 相 乘 次 数 ) × 3 ( A 2 列 数 ) + 1 ( A 1 A 2 行 数 ) × 3 ( A 1 A 2 每 行 元 素 和 A 3 每 列 元 素 对 应 相 乘 次 数 ) × 4 ( A 3 列 数 ) = 18 1(A\_1行数)\\times 2(A\_1每行元素和A\_2每列元素对应相乘次数)\\times3(A\_2列数) + 1(A\_1A\_2行数)\\times 3(A\_1A\_2每行元素和A\_3每列元素对应相乘次数)\\times 4(A\_3列数) = 18\\\\ 1(A1行数)×2(A1每行元素和A2每列元素对应相乘次数)×3(A2列数)\+1(A1A2行数)×3(A1A2每行元素和A3每列元素对应相乘次数)×4(A3列数)=18
 *   A 1 ( A 2 A 3 ) A\_1(A\_2A\_3) A1(A2A3)  
     2 ( A 2 行 数 ) × 3 ( A 2 每 行 元 素 和 A 3 每 列 元 素 对 应 相 乘 次 数 ) × 4 ( A 3 列 数 ) + 1 ( A 1 行 数 ) × 2 ( A 1 每 行 元 素 和 A 2 A 3 每 列 元 素 对 应 相 乘 次 数 ) × 4 ( A 2 A 3 列 数 ) = 32 2(A\_2行数)\\times 3(A\_2每行元素和A\_3每列元素对应相乘次数)\\times4(A\_3列数) + 1(A\_1行数)\\times 2(A\_1每行元素和A\_2A\_3每列元素对应相乘次数)\\times 4(A\_2A\_3列数) = 32\\\\ 2(A2行数)×3(A2每行元素和A3每列元素对应相乘次数)×4(A3列数)\+1(A1行数)×2(A1每行元素和A2A3每列元素对应相乘次数)×4(A2A3列数)=32  
    可以发现先算 A 1 A 2 A\_1A\_2 A1A2先计算所用的乘法次数要更少，因此原问题可以转换为一个最优化问题如下：
 *  **给定n个矩阵, A i A\_i Ai的规模是 p i − 1 p\_\{i-1\} pi−1行， p i p\_i pi列，确定 A 1 A 2 . . . A n A\_1A\_2...A\_n A1A2...An的计算顺序，使得需要的乘积计算次数最少（运算速度最快）**  
    我们说最优化问题是从所有可行解中找出最优，因此我们考虑所有可行解的个数。我们设n个矩阵运算可行解个数为 f ( n ) f(n) f(n)，考虑n和n-1的关系，可以先算前n-1个，再将n-1的结果和n相乘，也可以先算前n-2个，再算n-1和n，最后算前 A 1 A 2 . . . A n − 2 A\_1A\_2...A\_\{n-2\} A1A2...An−2和 ( A n − 1 A n ) (A\_\{n-1\}A\_n) (An−1An)，因此我们有如下递推关系：  
     f ( n ) = f ( n − 1 ) + f ( n − 2 ) f(n) = f(n-1)+f(n-2) f(n)=f(n−1)\+f(n−2)  
    很显然这是一个斐波那契数列，它的时间复杂度是 O ( 2 n ) O(2^n) O(2n)，即指数级的。
    
    > 该时间复杂度的求解并不难，用特征根法求出 f ( n ) f(n) f(n)即可发现，或者一项项展开也可发现规律。因此这里不做赘述
    
    这么高的时间复杂度，如果我们一个个枚举，从中找出最优显然不现实。因此我们考虑将复杂问题分为子问题。

## 解决思路 ##

### 初始版本 ###

既然要求 A 1 A 2 . . . A n A\_1A\_2...A\_n A1A2...An的运算顺序，我们将其分成两部分，先求左边的运算顺序，再求右边的运算顺序，最后将左右合并，但这里就体现了DP和分治不同的地方，我们并不知道如何划分才能获得最好的效果，因此划分点设置为变量 k k k（如果是分治就会选择中间点使得数据规模指数级下降），再对 k k k进行遍历比较，找出最优的。首先，设 O P T ( 1 , n ) OPT(1,n) OPT(1,n)表示 A 1 A 2 . . . A n A\_1A\_2...A\_n A1A2...An的最优运算次数，当确定分点 k k k时，我们有以下递推式：  
 O P T ( 1 , n ) = O P T ( 1 , k ) + O P T ( k + 1 , n ) + p 0 p k p n OPT(1,n) = OPT(1,k)+OPT(k+1,n) + p\_0p\_kp\_n OPT(1,n)=OPT(1,k)\+OPT(k\+1,n)\+p0pkpn  
其中 p 0 p k p n p\_0p\_kp\_n p0pkpn是将 ( 1 , k ) (1,k) (1,k)和 ( k + 1 , n ) (k+1,n) (k\+1,n)合并起来的开销，因为 A 1 A 2 . . . A k A\_1A\_2...A\_k A1A2...Ak是 p 0 p\_0 p0行， p k p\_k pk列的矩阵， A k + 1 . . . A n A\_\{k+1\}...A\_n Ak\+1...An是 p k p\_k pk行， p n p\_n pn列的矩阵，这两个矩阵乘法次数是 p 0 p k p n p\_0p\_kp\_n p0pkpn（自己可以动手算算）。然后，我们需要对 k k k的位置进行遍历，从中找到使得 O P T ( 1 , n ) OPT(1,n) OPT(1,n)最小的，因此有：  
 O P T ( 1 , n ) = m i n k = 1 n − 1 \{ O P T ( 1 , k ) + O P T ( k + 1 , n ) + p 0 p k p n \} OPT(1,n) = min\_\{k=1\}^\{n-1\}\\\{OPT(1,k)+OPT(k+1,n) + p\_0p\_kp\_n\\\} OPT(1,n)=mink=1n−1\{ OPT(1,k)\+OPT(k\+1,n)\+p0pkpn\}  
推广到一般形式(i,j)，得到整体递推公式如下：  
 O P T ( i , j ) = \{ 0 if  i = = j ( 如 果 是 一 个 矩 阵 不 用 计 算 ) m i n k = 1 n − 1 \{ O P T ( 1 , k ) + O P T ( k + 1 , n ) + p 0 p k p n \} else  OPT(i,j) = \\begin\{cases\} 0 &\\text\{if \} i==j (如果是一个矩阵不用计算) \\\\ min\_\{k=1\}^\{n-1\}\\\{OPT(1,k)+OPT(k+1,n) + p\_0p\_kp\_n\\\}&\\text\{else \} \\end\{cases\} OPT(i,j)=\{ 0mink=1n−1\{ OPT(1,k)\+OPT(k\+1,n)\+p0pkpn\}if i==j(如果是一个矩阵不用计算)else   
由此我们就可以设计出伪代码，如下：

RECURSIVE-MATRIX-CHAIN(P,i,j)
    	if i==j then 
    		return 0;
    	end if
    	OPT(i,j) = +INF  // 初始化无穷，保证会被更新
    	for k=i to j-1 do
    		v = RECURSIVE-MATRIX-CHAIN(P,i,k) + RECURSIVE-MATRIX-CHAIN(P,k+1,j) + p[i-1]p[k]p[j]
    		if v < OPT(i,j) then
    			OPT(i,j) = v;
    		end if
    	end for
    	return OPT(i,j)

程序实现起来很简单，但这个程序十分费时。原因在于，有很多重复计算的地方。这里以n=4为例，对 A 1 A 2 A 3 A 4 A\_1A\_2A\_3A\_4 A1A2A3A4划分，如下图a（这个图很有用，请记住）：  
![在这里插入图片描述][20191018113055460.png_pic_center]

> 其中相同颜色区域表示一种划分方式，红色节点表示重复计算（在前面的划分中都能找到计算结果）

可以发现实际的计算只需要9次就够了，很多结果被重复计算多次。这个导致了时间耗费严重。

### 改进版本 ###

既然是因为计算过的东西被重复计算造成耗时，那么我们将计算结果保存下来，就能避免重复计算了。可以发现上述树的计算次数和它的节点个数有关，我们用一个二维数组保存状态 （ i , j ） （i,j） （i,j）的计算结果，每一次树要往下分支的时候，先查询以下数组，如果已经存在，则直接返回结果，避免产生新节点。因此，我们有改善后的伪代码如下：

RECURSIVE-MATRIX-CHAIN(p,i,j)
    	if DP[i][j] exists then // 说明已经计算过
    		return DP[i][j]
    	if i==j then 
    		return 0;
    	end if
    	OPT(i,j) = +INF  // 初始化无穷，保证会被更新
    	for k=i to j-1 do
    		v = RECURSIVE-MATRIX-CHAIN(P,i,k) + RECURSIVE-MATRIX-CHAIN(P,k+1,j) + p[i-1]p[k]p[j]
    		if v < OPT(i,j) then
    			OPT(i,j) = v;
    		end if
    	end for
    	DP[i][j] = OPT(i,j) // 更新数组
    	return OPT(i,j)

> 这种用空间换时间的策略在动态规划中十分常见，它能减少的重复计算量是指数级别的！！！你不妨写一个斐波那契数列的程序，看看用DP数组存储和不用DP数组存储结果程序运行耗费的时间比较。

### 算法运行过程剖析 ###

前面已经完成了程序的设计，并对算法进行了优化，但实际上，算法内部运行过程才是我们关心的，因为这涉及运算顺序。为了理解简便，我们还是考虑n=4的情况。

*  求 A 1 A 2 A 3 A 4 A\_1A\_2A\_3A\_4 A1A2A3A4  
     O P T ( 1 , 4 ) = m i n \{ O P T ( 1 , 1 ) + O P T ( 2 , 4 ) + p 0 p 1 p 4 O P T ( 1 , 2 ) + O P T ( 3 , 4 ) + p 0 p 2 p 4 O P T ( 1 , 3 ) + O P T ( 4 , 4 ) + p 0 p 3 p 4 OPT(1,4) = min\\begin\{cases\} OPT(1,1) + OPT(2,4) + p\_0p\_1p\_4 \\\\ OPT(1,2) + OPT(3,4) + p\_0p\_2p\_4 \\\\ OPT(1,3) + OPT(4,4) + p\_0p\_3p\_4 \\end\{cases\} OPT(1,4)=min⎩⎪⎨⎪⎧OPT(1,1)\+OPT(2,4)\+p0p1p4OPT(1,2)\+OPT(3,4)\+p0p2p4OPT(1,3)\+OPT(4,4)\+p0p3p4
    
    可以发现要求 O P T ( 1 , 4 ) OPT(1,4) OPT(1,4)会进而转化为求 O P T ( 1 , 1 ) , O P T ( 2 , 4 ) . . . OPT(1,1),OPT(2,4)... OPT(1,1),OPT(2,4)...， O P T ( 2 , 4 ) OPT(2,4) OPT(2,4)会转换为求 O P T ( 2 , 2 ) , O P T ( 3 , 4 ) . . . OPT(2,2),OPT(3,4)... OPT(2,2),OPT(3,4)...这个递归会不断进行，直到每个OPT能够被直接求出，就递归返回结果。因此我们可以构造如下表格：  
    ![在这里插入图片描述][20191018205649727.png_pic_center]
    
    > 整个计算过程有三轮，其中橙色背景区域表示已经获得计算结果（也可以理解为DP数组已经被更新），在第一轮中，需要计算间隔为1的值(12,23,34)，在第二轮中需要计算间隔为2的值(13,24),最后一轮，计算间隔为3的值(14)。左侧的表格记录了 A i A j A\_iA\_j AiAj的运算次数，右侧的表格记录了选择的分割点。如求 A 1 A 2 A\_1A\_2 A1A2时，分割点选择了1，由公式可知 O P T ( 1 , 2 ) = O P T ( 1 , 1 ) + O P T ( 2 , 2 ) + p 0 × p 1 × p 2 = 0 + 0 + 1 × 2 × 3 = 6 OPT(1,2) = OPT(1,1)+OPT(2,2) + p\_0\\times p\_1\\times p\_2 = 0+0+1\\times 2\\times 3 = 6 OPT(1,2)=OPT(1,1)\+OPT(2,2)\+p0×p1×p2=0\+0\+1×2×3=6,同理可得 O P T ( 2 , 3 ) , O P T ( 3 , 4 ) OPT(2,3),OPT(3,4) OPT(2,3),OPT(3,4)，有了这三个就使得 O P T ( 1 , 3 ) , O P T ( 2 , 4 ) OPT(1,3),OPT(2,4) OPT(1,3),OPT(2,4)得以计算得到结果。计算的顺序过程建议参考上文的图a，父节点的计算结果依赖于子节点的计算结果。

### 用迭代替代递归 ###

我们再思考上述过程，递归从根节点实际上是先往下走，再往上走，当我想知道 A 1 . . . A 4 A\_1...A\_4 A1...A4时，我会向子节点询问，子节点如果没有直接结果，则再继续向孙子节点询问。在询问的过程中，为了避免重复计算，则用一个数组用于存储询问的中间结果。既然如此，我们完全可以从底层开始，子节点主动向自己的父节点汇报自己的计算结果，父节点再向祖先节点汇报，最终根节点汇报它的最终结果就是我们想要的答案，就可以避免自上而下的过程。由此，我们可以把递归的伪代码转换为三重for循环(也就是用迭代替换递归)：

ITERATIVE-MATRIX-CHAIN(p):
    	set OPT(i,i) = 0 for all i(1<=i<=n) // 对角线运算次数为0
    	for l=2 to n do    
    		for  i=1 to n-l+1 do
    			j = i+l-1
    			OPT(i,j) = +INF
    			for k=i to j-1 do 
    				v = OPT(i,k) + OPT(k+1,j) + p[i-1]p[k]p[j]
    				if v < OPT(i,j) then
    					OPT(i,j) = v
    					SPLITTER(i,j) = k  // 记录分割点
    				end if
    			end for
    		end for
    	end for
    return OPT(1,n)  // 返回根节点的计算结果

这里实际上就是添加了两个外层for循环用于更新OPT这个数组，第一层l（2-n）指的是 A i . . . A j A\_i...A\_j Ai...Aj中j的变化过程，我们需要考虑的有 A 1 A 2 , A 1 A 3 , A 2 A 3 , . . . A\_1A\_2,A\_1A\_3,A\_2A\_3,... A1A2,A1A3,A2A3,...j从2开始一直到n，第二层循环的i，j实际上保证了每次取的都是l对应列的对角线元素。如果不好理解的可以看下图：  
![在这里插入图片描述][20191018213202418.png_pic_center]

注意到过程中我们有额外的一个数组SPLITTER用于存储分割点，例如SPLITTER(1,4) = 3意味着当对 A 1 A 2 A 3 A 4 A\_1A\_2A\_3A\_4 A1A2A3A4计算时，以3为分割点，先计算 A 1 A 2 A 3 A\_1A\_2A\_3 A1A2A3

> 算法的时间复杂度是 O ( n 3 ) O(n^3) O(n3)，可以看出由于避免了递归调用的额外开销，迭代的方式实现要比递归调用更快。但这不意味着迭代就一定比递归快，在迭代中，所有的子问题都会被要求计算一遍，从而得到原问题的解，但递归的时候可能某些子问题我们并不需要计算（只计算实际要用到的），另一方面，递归的方式在程序编写上要更加简洁易懂，这也是其优势。

# 总结 #

本文通过矩阵链式相乘问题来引入动态规划的实际应用，可以发现在应用动态规划的时候我们能够将一个原本指数级的计算问题转换成多项式时间的计算问题。为什么能做到这一点？我们对比递归算法的初始版本和改进版本，发现我们将很多冗余的计算用数组存储结果，从而节省了很多时间（这个节省量是很大的，还记得前文我说过是指数级的）。如果仅仅只用递归，那仅仅只是简便了程序的实现而已。

实际上动态规划中很有名的就是背包问题（我在文首也提到过），一个物品装与不装是2个选择，n个问题使得问题规模可以达到指数级，而我们可以通过动态规划将这个计算规模大幅下降。总结以下几点：

*  1、碰到最优化问题，首先考虑是否是多步决策
    
     *  （1）解能否被逐步构建（这与分治思想同源）
     *  （2）目标函数可否分（优化问题独有），在矩阵相乘的例子中目标函数就是 A 1 A 2 . . . A n A\_1A\_2...A\_n A1A2...An
 *  2、 动态规划为什么快？避免了许多冗余计算！！！

> 如果你觉得文章对你有用，不妨顺手点个赞哦～

*  上一篇：[算法设计与分析：分治思想（4）- 快速傅立叶变换（对数组的归约）][4_-]
 *  下一篇：[算法设计与分析：动态规划（2） - 序列联配问题(最优对齐方案)][2_ -]
 *  [目录][Link 1]

[20191018113055460.png_pic_center]: /images/20230601/a76bc5917d0a4029bf136b15af08aecb.png
[20191018205649727.png_pic_center]: /images/20230601/d76a65bf48f94490b8770792623a2221.png
[20191018213202418.png_pic_center]: /images/20230601/cd56b05f6bf84c7d9780dd9d1d973a3d.png
[4_-]: https://blog.csdn.net/GentleCP/article/details/102246489
[2_ -]: https://blog.csdn.net/GentleCP/article/details/102779619
[Link 1]: https://blog.csdn.net/GentleCP/article/details/102609055