25 Julia编程示例–动态规划
动态规划(dynamic programming)不完全是数学规划（最优化）问题的算法。 它能解决如下问题： 问题可以分解为子问题， 子问题多次重复出现， 就可以将已经出现的子问题记住， 从而大大减小计算量。
25.1 Fibonacci数列
Fibonacci数列满足\(F_0=F_1=1\)， \(F_n = F_{n-1} + F_{n-2}\)， \(n=2,3,\ldots\)。
显然，可以写成简单的递归函数：
function fib_v1(n)
if n==0 || n==1
return 1
else
return fib_v1(n-1) + fib_v1(n-2)
end
end
fib_v1(10)
## 89但是，这里面有许多额外不必要的计算。 比如，算\(F_{10}\)用到\(F_9\)和\(F_8\)， 算\(F_8\)时用到\(F_7\)和\(F_6\)， 算\(F_9\)时又重复计算了\(F_8\)和\(F_7\)，等等。
改进的思想是每计算一个\(F_k\)就将它记住。 具体的算法有从顶向下和从底向上两种。
从顶向下方法：
function fib_v2(n, dict=Dict())
if !haskey(dict, n)
dict[n] = (n==0 || n==1) ? 1 :
fib_v2(n-1, dict) + fib_v2(n-2, dict)
end
return dict[n]
end从底向上方法：
function fib_v3(n)
local f0, f1, f2
f0 = f1 = 1
for i=2:n
f2 = f0 + f1
f1, f0 = f2, f1
end
return f2
end25.2 最优打包方案
设有\(n\)件物品， 价值分别为\(v_i\), 重量分别为\(w_i\)， \(i=1,2,\ldots,n\)。 假设在某种情况下， 只能在重量不超过\(w_{\max}\)的条件下带走这些物品中价值最大的物品子集， 可以写成如下的一个最优化问题：
\[\begin{aligned} & \min \left(- \sum_{i=1}^n v_i x_i \right), \text{ s.t.} \\ & \sum_{i=1}^n w_i x_i \leq w_{\max}, \\ & x_i \in \{0, 1\}, \ \forall i=1,2,\ldots,n . \end{aligned}\]
其中\(x_i=1\)表示装进来此项， \(x_i=0\)表示不装此项。 如果用穷举法， 需要考虑\(2^n\)种组合， 当\(n\)很大时不可行。 设\(w_i\), \(w_{\max}\)都是整数， 可以用动态规划方法求解。 将\(n\)和\(w_{\max}\)看成是变量， 在求解关于\(n\)项和限制\(w_{\max}\)的问题时， 将其拆分为已解决前\(n-1\)项物品的问题， 和第\(n\)项是否要装进来的问题。 要记住已经解决的子问题。
用从底向上的递归和记忆已解决的子问题的方法。 如果\(n=0\)， 则不管重量限制是多少， 带走物品的价值都等于0。
对于从已解决前\(n-1\)项如何打包， 问第\(n\)项是否装入的问题， 先解决前\(n-1\)项物品， 重量限制从0到\(w_{\max}\)的每一个整数值的子问题并记住解。 然后考虑在重量限制为\(w_{\max}\)的条件下第\(n\)项是否应该装进来的问题， 这时，如果第\(n\)项重量已经超过\(w_{\max}\)， 则当然不能装此项； 否则， 就比较一下不装此项的总价值， 与装此项，使得其它\(n-1\)项的总重量不超过\(w_{\max}-w_n\)的最优总价值加上此项的价值， 取其中价值较高的一个方案。 程序：
# v: 价值
# w: 重量
# w_max: 重量限制
function knapsack(v, w, w_max)
n = length(v) # 物品数
# y初值：物品数为0的子问题的解，
# 每一个解是物品个数与重量限制对应的最大价值
y = Dict((0,j) => 0.0 for j in 0:w_max)
for i in 1 : n # 对每一个新增物品
for j in 0 : w_max # 对每一个重量限制
# y[i,j]是前i个物品，设重量限制为j时的最大价值
# 首先考虑第i个物品是否根本在重量限制j下无法装，
# 这时直接不装第i个物品，最大价值为y[i-1,j];
# 否则，有两种选择：
# 不装第i个物品，最大价值为y[i-1,j];
# 装第i个物品，前i-1个物品的重量限制缩小为j-w[i]，
# 最大价值为第i个物品的价值，加上前i-1个物品在重量限制j-w[i]下的最大价值。
# 取这两种选择中较大一个。
y[i,j] = w[i] > j ? y[i-1,j] : max(
y[i-1,j], y[i-1,j-w[i]] + v[i])
end
end
# 从字典y中找到解，
# 注意字典中(n, w_max)对应的项已经是最大价值，
# 我们只需找出对应的装法x，
# 从第n个物品逐个向前求解
# x: 用false表示不装，用true表示装入
# j：从后向前装入部分物品后还剩余的容量
x, j = falses(n), w_max
for i in n: -1 : 1
if w[i] ≤ j && y[i,j] - y[i-1, j-w[i]] == v[i]
# 如果第i个物品不超过剩余可用重量，
# 且前i个物品的最大价值，
# 恰好等于前i-1个物品
# 在假设第i物品占用重量限额条件下的最大价值
# 加上第i物品价值，则表示第i物品被装入
x[i] = true
# 将第i物品的重量从可用重量额度中扣除
j -= w[i]
end
end
return x
end参考：
- Mykel J. Kochenderfer, Tim A. Wheeler. Algorithms for Optimization, §19.5, MIT press, 2019.
25.3 最小能量路径
25.3.1 问题
动态规划是最优化问题的一种算法。 最优化涉及最优决策问题， 此决策问题经常可以分解为若干步骤， 要解决许多子问题， 而这些子问题往往是重复的， 记住已经解决的子问题， 会使得决策问题大大简化， 这是动态规划问题的主要思想。
下面的例子来自MIT计算思想公开课课件。
考虑一个\(n \times n\)的随机矩阵，如：
using Random
Random.seed!(101)
M = rand(0:9, 5, 6)5×6 Matrix{Int64}:
7 7 9 9 7 1
1 3 7 3 3 2
3 6 5 6 1 8
5 1 1 7 6 6
9 2 5 3 9 6考虑从第一行到最后一行的路径， 每条路径从第一行的某个元素出发， 然后进入下一行的正下方或正下方的左侧元素、或右侧元素， 直到进入最后一行的元素为止， 每条路径有\(n\)个元素， 可以计算路径上元素总和。 将\(M\)的\((i,j)\)元素看作是该位置的“能量”， 试图找到从顶层（第一行）某个位置出发的“最小能量”路径。 可以考虑一般的\(m \times n\)矩阵。
25.3.2 穷举法
穷举法可以找到最小能量路径， 还可以罗列所有的路径。 穷举法的计算量将随着\(n\)增大而指数增加。
穷举法的程序：
function matdp_exaust(M; DEBUG=true)
m, n = size(M)
paths = []
# 从第i行开始的递归方法
# 直接修改外层的paths
# 输入为截止到第i行的部分路径path_in和i
function path_remain(path_in, i)
if(i==m) # 最后一行已解决
push!(paths, copy(path_in))
return
else # i < m, 考虑下一行的三种不同情况
# 正下方
path = copy(path_in)
path[(i+1):end] .= 0
path[i+1] = path[i] # 正下方
path_remain(path, i+1)
# 左下方
if path[i] > 1
path = copy(path_in)
path[(i+1):end] .= 0
path[i+1] = path[i] - 1
path_remain(path, i+1)
end
# 右下方
if path[i] < n
path = copy(path_in)
path[(i+1):end] .= 0
path[i+1] = path[i] + 1
path_remain(path, i+1)
end
end # if i==n
# 函数不需要显式地返回值，通过外层变量paths记录遍历结果
end # function path_remain
# 只要遍历第一行
path = fill(0, m)
for start=1:n
path[1] = start
path_remain(path, 1)
end
#println(paths)
npaths = length(paths)
sums = map(sum, [[M[i, p[i]] for i=1:m] for p in paths])
id = argmin(sums)
minen = sums[id]
if(DEBUG)
println("矩阵大小：", (m, n), " 路径数：", npaths)
end
return paths[id], minen
end测试运行：
function test_matdp(m=8, n=m, solver=matdp_exaust; DEBUG=true)
Random.seed!(101)
M = rand(0:9, m, n)
bestpath, minen = solver(M; DEBUG=DEBUG)
if(DEBUG)
end
if(DEBUG)
println("矩阵大小：", (m, n))
println("最小能量和：", minen)
println("最小能量路径：")
println(bestpath)
M
end
end
test_matdp(5, 6, matdp_exaust, DEBUG=true)矩阵大小：(5, 6) 路径数：340
矩阵大小：(5, 6)
最小能量和：13
最小能量路径：
[6, 6, 5, 5, 4]
5×6 Matrix{Int64}:
7 7 9 9 7 1
1 3 7 3 3 2
3 6 5 6 1 8
5 1 1 7 6 6
9 2 5 3 9 6因为是路径上的值是整数值， 路径和为整数值， 所以最小能量路径不一定唯一， 我们只输出了其中一个， 但是因为是穷举法， 所以需要的话也可以输出达到最小能量的所有路径。
25.3.3 递归解法
如果不需要罗列所有路径， 也可以使用递归算法， 但仅考虑达到最小能量的路径。
## 仅考虑一个最优解的递归方法
function matdp_recurse(M; DEBUG=true)
m, n = size(M)
# 递归函数。
# 返回从(i,j)元素向下， 找到的最小能量， 和下一行应去的列号。
function next_ij(M, i, j)
m, n = size(M)
## 初始化情形
if i == m
return M[m, j], n+1
end
## 递归求解情形
# 考虑第i+1行的j-1, j, j+1对应的最小能量
js = unique(max.(1, min.(n, (j-1):(j+1))))
min_ens = [next_ij(M, i+1, jj)[1] for jj in js]
ind = argmin(min_ens)
jnext = js[ind]
min_en = min_ens[ind] + M[i, j]
return min_en, jnext
end
# 只要从第一行开始找到最优开始位置，然后依次下行就可以
minens = [next_ij(M, 1, j)[1] for j=1:n]
j1 = argmin(minens)
minen = minens[j1] # 全局最小能量
bestpath = [j1;]
for i = 2:m
en, j = next_ij(M, i-1, bestpath[i-1])
push!(bestpath, j)
end
return bestpath, minen
end25.3.4 动态规划方法
上述穷举法的缺点是反复计算了许多不必要的下级总和。 以\(6 \times 5\)的如下矩阵为例：
7 7 9 9 7 1
1 3 7 3 3 2
3 6 5 6 1 8
5 1 1 7 6 6
9 2 5 3 9 6上述矩阵的一个最优路径是从上到下列数为： \[ 6, 5, 4, 3, 2 \] 对应的值为 \[ 1, 3, 6, 1, 2 \] 最优值为13。 考虑第3行第3列的元素5； 设某一条路径经过\(M[3,3]=5\)这个点， 则从这个点出发向下易见有9条路径，经过的列号为： \[ (3,2,1), (3,2,2), (3,2,3), (3,3,2), (3,3,3), (3,3,4), (3,4,3), (3,4,4,), (3,4,5) \] 可以计算出从\(M[3,3]=5\)的这个点向下的所有路径以及路径对应的部分和， 部分最优路径为\((3,2,2)\)，部分最小能量为8。 从第一行经过\(M[3,3]=5\)的路径也有9条， 但是可以一次性地计算出从\(M[3,3]=5\)的这个点向下的结果而不必重复计算， 这就是动态规划的思想。
如果仅要求极值而不需要遍历所有路径， 则可以仅考虑从\(M[3,3]=5\)的这个点向下的一条最优路径而不需要考虑任何其它路径， 即将\(M[3,3]=5\)与从这个点向下的最优路径\((5,1, 3)\)记住就可以了。 我们从最后一行开始， 最后一行的\(n\)个点都可以作为一条子路径； 从最后一行(第\(m\)行)到第\(m-1\)行， 则第\(m-1\)行的\(n\)个元素， 每一个元素都只需要考虑下方的2个或3个最优子路径， 从而找到第\(m-1\)行的每个元素的最优子路径和向下的最小能量， 然后再考虑第\(m-2\)行，如此倒推一直到第一行的\(n\)个元素都找到最优子路径， 这时比较第一行的\(n\)个最优子路径然后找到最优即可。 这样可以将计算复杂度降低到\(O(n)\)级别而不是指数级别。 如果需要找到所有最优子路径也很容易。 用这种思想求最优路径的函数：
## 用动态规划方法求解的函数
function matdp_dp(M; DEBUG=true)
m, n = size(M)
# 对M的每个元素，求出并记住最优路径
MM = Array{Any}(undef, (m, n))
# MM[i,j] 存储一条最优子路径和最优子路径的部分和
# 最后一行，最优路径即本身
for j = 1:n
MM[m, j] = ([j,], M[m,j])
end
# 第(m-1)到第1行：
for i=(m-1):-1:1
for j = 1:n
# 考虑M[i,j]的最优路径
# 只要比较MM[i+1, j-1], MM[i+1, j], MM[i+1, j+1]
js = [max(1, j-1), j, min(n, j+1)]
# js: 下一行的三个可能列号
vs = [MM[i+1, jj][2] for jj in js]
# vs: 下一行的三个可能路径的部分和
jmin = js[argmin(vs)]
# jmin: 下一行的最优列号
vmin = MM[i+1, jmin][2]
# vmin: 下一行的最优值
MM[i, j] = ([j; MM[i+1,jmin][1]], M[i,j] + vmin)
end # for j
end # for i
# 比较第一行的n个最优值，找到全局最优路径
best1 = [MM[1,j][2] for j=1:n]
bestj = argmin(best1)
bestpath = MM[1, bestj][1]
besten = MM[1, bestj][2]
return bestpath, besten
end用BenchmarkTools测试比较：
using BenchmarkTools
@benchmark test_matdp(10, 8, matdp_exaust, DEBUG=false)
@benchmark test_matdp(10, 8, matdp_recurse, DEBUG=false)
@benchmark test_matdp(10, 8, matdp_dp, DEBUG=false)结果显示穷举法运行时间中位数为40ms, 递归法运行时间中位数为176ms， 动态规划法运行时间中位数为109μs， 动态规划方法比递归发速度快一千多倍， 穷举法和递归法相差几倍。 动态规划法因为避免了大量重复计算所以速度快得多， \(n\)越大这个差距越大， \(n\)较大时穷举法和递归发已经不能在合理时间内完成算法。
25.3.5 使用自定义的迭代器
关于穷举法，
MIT计算思维公开课的Julia编程使用了自定义迭代器的思想。
自定义的迭代器必须实现Base.iterate(itr)和Base.iterate(itr, state)，
可选定义Base.IteratorSize()函数说明该迭代器类型的长度类型。
这样就可以遍历所有路径。
关于迭代器参见：
具体遍历方法是， 从取第一列为路径开始， 一直到路径为最右边一列结束。 在中间的某一路径， 求得下一路径的方法是， 从最后一行向上倒序查找， 直到当前路径在当前行的点可以向右移动一个格子而且不会与上一行割裂， 这样就右移这个点， 然后保留上面所有行的路径， 从右移的这个点开始找到最左下的向下路径。
# MIT公开课用迭代器进行穷举的方法
struct Paths # 用来表示所有路径的类型
m::Int # 矩阵行数
n::Int # 矩阵列数
end
# 要自定义迭代器，必须定义Base.iterate(itr)和Base.iterate(itr, state)
# 用来初始化的函数，从矩阵最左边一列开始
Base.iterate(p::Paths) = fill(1, p.m), fill(1, p.m)
# 用来定义迭代器长度的类型
Base.IteratorSize(::Type{Paths}) = SizeUnknown()
# 用来迭代进入下一个状态的函数
function Base.iterate(p::Paths, state)
if state ≠ fill(p.n, p.m) # 当路径为最后一列时结束
newstate = next(state, p.n)
return newstate, newstate
end
end
# 给定一条路径，求下一条路径，n为列数
function next(path, n)
m = length(path) # 行数
k = m # 从最后一行开始向上查找
# 从现有路径的最后一行向上查找，直到找到一个可以向右移动一格的位置
# 下面的while结构是如果第k行还没有到第一行，
# 并且当前路径第k行已经在最后一列，
# 或者当前路径第k行如果向右一列会与第k-1行相比过于靠右无法连接，
# 则第k行不变，考虑上一行即第k-1行
while k≥2 && ( path[k]==n || path[k]+1 > path[k-1]+1 )
k -= 1
end
# 上面的循环找到了适当的行号k，
# 此行的当前路径点可以右移且不会与第k-1行割裂，
# 于是右移第k行的这个点，第1:(k-1)行保持不变
path[k] +=1
# 前面用倒序从最后一行开始找到了可以右移且不会与前面行割裂的一个行号k，
# 将当前路径第k行右移一个格，
# 然后找到从第k行开始的最左边的一条向下路径
for j = (k+1):m
# 每一行依次向左下方
path[j] = max(path[j-1]-1,1)
end
return(path)
end
# 用迭代器获得所有路径
# m: 矩阵行数
# n: 矩阵列数
function allpaths(m,n)
v = Vector{Int}[] # v的每个元素是一个整数向量，为一条路径
paths = Paths(m,n) # 一个迭代器
for p ∈ paths
push!(v, copy(p))
end
return v
end
# 返回最优路径和最优值的函数
function matdp_iter(M; DEBUG=true)
(m, n) = size(M)
paths = allpaths(m, n)
npaths = length(paths)
sums = map(sum, [[M[i, p[i]] for i=1:m] for p in paths])
id = argmin(sums)
minen = sums[id]
if(DEBUG)
println("矩阵大小：", n, " 路径数：", npaths)
end
return paths[id], minen
end # function matdp_iter