该书由清华大学李升波教授撰写的,主要面向工业控制领域的研究者和工程师,曾获得2024年度Springer中国新发展奖(China New Development Awards)。全书按照原理剖析、主流算法、典型示例的架构,系统地介绍了用于动态系统决策与控制的强化学习方法。全书共分为11章,内容涵盖了强化学习的基本概念、蒙特卡洛法、时序差分法、动态规划法、函数近似法、策略梯度法、近似动态规划、状态约束的处理和深度强化学习等知识点。书籍及源代码下载网站:书籍及代码链接点这里。
另外,由于本部分篇幅较长,因此分成三篇博客发表。其余两篇博客的地址参考我的系列博客的地址汇总。
文章目录
书籍链接:Reinforcement Learning for Sequential Decision and Optimal Control
本篇博客主要介绍Dynamic Programming(DP),即动态规划,一种Model-Based的Indirect RL方法。DP与动态环境中的随机最优控制问题有紧密的联系。DP通常通过求解Bellman方程(离散时间)或Hamilton-Jacobi-Bellman方程(连续时间)来得到最优策略。对于大部分工程任务来说,以上两个方程是最优解的充要条件。DP的优点在于它的强大的适应性。如线性/非线性系统、离散/连续系统、确定性/随机性系统、可微分/不可微分系统等。但是,DP在实际应用中面临着维度灾难和缺少泛化能力的问题。但是,DP仍然在理论分析最优解的存在以及收敛速度等方面有着重要的作用。
5.2 DP的策略迭代(Policy Iteration)算法
5.2.1 Policy Evaluation(PEV)
不同与MC、TD等model-free RL那里的PEV是通过采样得到的样本来更新值函数的,像DP这类model-based RL的PEV是通过迭代求解self-consistency条件来更新值函数的。先来回顾一下第一类self-consistency条件:
v
π
(
s
)
=
∑
a
∈
A
π
(
a
∣
s
)
{
∑
s
′
∈
S
P
s
s
′
a
(
r
s
s
′
a
+
γ
v
π
(
s
′
)
)
}
,
∀
s
∈
S
,
v^{\pi}(s)=\sum_{a\in\mathcal{A}}\pi(a|s)\left\{\sum_{s^{\prime}\in\mathcal{S}}P_{ss^{\prime}}^{a}\left(r_{ss^{\prime}}^{a}+\gamma v^{\pi}(s^{\prime})\right)\right\},\forall s\in\mathcal{S},
vπ(s)=a∈A∑π(a∣s){s′∈S∑Pss′a(rss′a+γvπ(s′))},∀s∈S,
那么下面的PEV算法就是利用上述第一类self-consistency条件来更新值函数的:
这里当下标j趋向于无穷的时候就可以得到在当前策略下值函数的的准确值。注意,在上述迭代过程中策略始终不变。其实在实际计算中,并不用到无穷,而是当更新量小于某个阈值之后即可停机。这里正好再解释一下为什么DP被称为model-based的算法。注意到上述迭代公式中有一个
P
\mathcal{P}
P,这个
P
\mathcal{P}
P就是环境的模型,即环境的状态转移概率,而为了知道这个概率就需要要我们对于环境模型事先有充足的了解才行。再看看之前讲过的model-free RL,它就不需要知道环境的模型,只是根据sample来估计。举我们之前讲过的Q-learning的更新公式为例:
Q
(
s
,
a
)
←
Q
(
s
,
a
)
+
α
(
r
+
γ
max
a
Q
(
s
′
,
a
)
−
Q
(
s
,
a
)
)
Q(s,a)\leftarrow Q(s,a)+\alpha\left(r+\gamma\max_{a}Q(s',a)-Q(s,a)\right)
Q(s,a)←Q(s,a)+α(r+γamaxQ(s′,a)−Q(s,a)) 可以发现这个更新公式中并没有出现环境的模型
P
\mathcal{P}
P。因此,有没有在算法中出现环境的模型就是model-based和model-free的一个重要区别。下面讲到的DP的值迭代的更新公式可以按照这个角度对照参看,看看其是否符合这里讲的model-based的特点。
下面从fixed-point iteration的角度来理解一下PEV算法。注意到第一类self-consistency条件的式子并不满足fixed-point iteration的标准形式,因为它左右两侧的变量一个是值函数
v
π
(
s
)
v^{\pi}(s)
vπ(s),一个是
v
π
(
s
′
)
v^{\pi}(s^{\prime})
vπ(s′)。因此,我们需要对上式进行一些变形。首先,设状态空间有n个可能的取值,则我们将对应的n个方程写出来:
V
π
(
s
(
1
)
)
=
∑
π
∑
P
s
(
1
)
s
(
1
)
′
a
(
r
+
γ
V
π
(
s
(
1
)
′
)
)
,
V
π
(
s
(
2
)
)
=
∑
π
∑
P
s
(
2
)
s
(
2
)
′
a
(
r
+
γ
V
π
(
s
(
2
)
′
)
)
,
⋮
V
π
(
s
(
n
)
)
=
∑
π
∑
P
s
(
n
)
s
(
n
)
′
a
(
r
+
γ
V
π
(
s
(
n
)
′
)
)
,
\begin{gathered} V^{\pi}\big(s_{(1)}\big) =\sum\pi\sum\mathcal{P}_{s_{(1)}s_{(1)}^{\prime}}^a\Big(r+\gamma V^\pi\Big(s_{(1)}^{\prime}\Big)\Big), \\ V^\pi(s_{(2)}) =\sum\pi\sum\mathcal{P}_{s_{(2)}s_{(2)}^{\prime}}^a\Big(r+\gamma V^\pi\Big(s_{(2)}^{\prime}\Big)\Big), \\ \vdots \\ V^{\pi}\big(s_{(n)}\big) =\sum\pi\sum\mathcal{P}_{s_{(n)}s_{(n)}^{\prime}}^{a}\Big(r+\gamma V^{\pi}\Big(s_{(n)}^{\prime}\Big)\Big), \end{gathered}
Vπ(s(1))=∑π∑Ps(1)s(1)′a(r+γVπ(s(1)′)),Vπ(s(2))=∑π∑Ps(2)s(2)′a(r+γVπ(s(2)′)),⋮Vπ(s(n))=∑π∑Ps(n)s(n)′a(r+γVπ(s(n)′)), 然后将这n个方程构成的方程组写成矩阵形式:
X
=
γ
B
X
+
b
,
X
=
[
V
π
(
s
(
1
)
)
,
V
π
(
s
(
2
)
)
,
⋯
,
V
π
(
s
(
n
)
)
]
T
∈
R
n
,
B
=
{
B
i
j
}
n
×
n
∈
R
n
×
n
,
b
=
{
b
i
}
n
×
1
∈
R
n
,
B
i
j
=
∑
a
∈
A
π
(
a
∣
s
(
i
)
)
P
s
(
i
)
s
(
j
)
a
,
b
i
=
∑
a
∈
A
π
(
a
∣
s
(
i
)
)
∑
j
P
s
(
i
)
s
(
j
)
a
r
s
(
i
)
s
(
j
)
a
.
\begin{gathered} X=\gamma BX+b, \\ X=\begin{bmatrix}V^{\pi}\big(s_{(1)}\big),V^{\pi}\big(s_{(2)}\big),\cdots,V^{\pi}\big(s_{(n)}\big)\end{bmatrix}^{\mathrm{T}}\in\mathbb{R}^{n}, \\ B=\begin{Bmatrix}B_{ij}\end{Bmatrix}_{n\times n}\in\mathbb{R}^{n\times n}, \\ b=\{b_i\}_{n\times1}\in\mathbb{R}^n, \\ B_{ij}=\sum_{a\in\mathcal{A}}\pi(a|s_{(i)})\mathcal{P}_{s_{(i)}s_{(j)}}^{a},\\ b_{i}=\sum_{a\in\mathcal{A}}\pi(a|s_{(i)})\sum_j\mathcal{P}_{s_{(i)}s_{(j)}}^{a}r_{s_{(i)}s_{(j)}}^a. \end{gathered}
X=γBX+b,X=[Vπ(s(1)),Vπ(s(2)),⋯,Vπ(s(n))]T∈Rn,B={Bij}n×n∈Rn×n,b={bi}n×1∈Rn,Bij=a∈A∑π(a∣s(i))Ps(i)s(j)a,bi=a∈A∑π(a∣s(i))j∑Ps(i)s(j)ars(i)s(j)a. 这里因为
π
\pi
π和环境的模型
P
\mathcal{P}
P是已知的,所以
B
B
B和
b
b
b都是定值。下面我们来利用fixed-point iteration的相关理论来解释两点:
- PEV算法的收敛性:在每轮迭代是,都会扫过整个状态空间,因此本算法的收敛性可以通过Bananach定理来保证。
- 上述矩阵方程只有一个fixed-point:定义矩阵A:
A
=
d
e
f
I
n
×
n
−
γ
B
.
A\overset{\mathrm{def}}{\operatorname*{=}}I_{n\times n}-\gamma B.
A=defIn×n−γB. 如果A满秩则存在唯一解。考虑到我们选取的
γ
\gamma
γ在(0,1)之间,所以由如下定理可得A满秩:
定理3:If 0 < γ < 1 , r a n k ( A ) = n . \text{定理3:}\text{If }0<\gamma<1,\mathrm{rank}(A)=n. 定理3:If 0<γ<1,rank(A)=n. 证明如下:首先,易知 0 ≤ B i j ≤ 1 0\le B_{ij}\le 1 0≤Bij≤1且矩阵B每行的和都为1(这可通过 ∑ j = 1 n B i j = \sum_{j=1}^nB_{ij}= ∑j=1nBij= ∑ j = 1 n ∑ a ∈ A π ( a ∣ s ( i ) ) P s ( i ) s ( j ) a \sum_{j=1}^n\sum_{a\in\mathcal{A}}\pi(a|s_{(i)})\mathcal{P}_{s_{(i)}s_{(j)}}^{a} ∑j=1n∑a∈Aπ(a∣s(i))Ps(i)s(j)a通过双重求和遍历了状态空间里所有i、j得到)。然后,易得 ∥ B ∥ ∞ = max i ∣ ∑ j = 1 n B i j ∣ = 1 \|B\|_{\infty}=\max_{i}|\sum_{j=1}^{n}B_{ij}|=1 ∥B∥∞=imax∣j=1∑nBij∣=1 这里 ∥ ⋅ ∥ ∞ \|\cdot\|_{\infty} ∥⋅∥∞表示矩阵的无穷范数,即矩阵元素的绝对值按行求和后取最大值 ∥ A ∥ ∞ = max 1 ≤ i ≤ n ∑ j = 1 n ∣ a i j ∣ \left\|A\right\|_{\infty}=\max_{1\leq\mathrm{i}\leq n}\sum_{j=1}^{n}\left|a_{ij}\right| ∥A∥∞=max1≤i≤n∑j=1n∣aij∣)。那么因为 γ \gamma γ为常数,所以有: ∥ γ B ∥ ∞ = γ ∥ B ∥ ∞ = γ < 1. \|\gamma B\|_\infty=\gamma\|B\|_\infty=\gamma<1. ∥γB∥∞=γ∥B∥∞=γ<1. 接下来,需要引入矩阵的谱半径的概念: ρ ( A ) = d e f max { ∣ λ ∣ : λ is an eigenvalue of A } . \rho(A)\overset{\mathrm{def}}{=}\max\left\{|\lambda|:\lambda\text{ is an eigenvalue of }A\right\}. ρ(A)=defmax{∣λ∣:λ is an eigenvalue of A}. 那么,有如下不等式: ρ ( γ B ) ≤ ∥ γ B ∥ ∞ < 1 , \rho(\gamma B)\leq\|\gamma B\|_\infty<1, ρ(γB)≤∥γB∥∞<1, 因此, γ B \gamma B γB的特征值中的最大值都小于1。又因为 I n × n I_{n\times n} In×n的特征值都是1,因此A的特征值均大于0。所以A是满秩的。证毕。
5.2.2 Policy Improvement(PIM)
这里主要是通过greedy的搜索来从PEV得到的值函数里面获得一个更优的策略。即:
a
∗
=
arg
max
a
q
π
(
s
,
a
)
.
a^*=\arg\max_aq^\pi(s,a).
a∗=argamaxqπ(s,a).
π
′
(
a
∣
s
)
=
{
1
,
if
a
=
a
∗
0
,
if
a
≠
a
∗
\left.\pi^\prime(a|s)=\left\{\begin{array}{ll}1,&\text{if }a=a^*\\\\0,&\text{if }a\neq a^*\end{array}\right.\right.
π′(a∣s)=⎩
⎨
⎧1,0,if a=a∗if a=a∗或写成:
π
′
(
s
)
=
arg
max
π
′
q
π
(
s
,
π
′
(
s
)
)
,
∀
s
∈
S
.
\pi^{\prime}(s)=\arg\max_{\pi^{\prime}}q^{\pi}\Big(s,\pi^{\prime}(s)\Big),\forall s\in\mathcal{S}.
π′(s)=argπ′maxqπ(s,π′(s)),∀s∈S. Greedy搜索得到的策略
π
′
\pi^{\prime}
π′是比原策略
π
\pi
π更优的,这里证明已经在第二单元的博客里证明过了。详情可以参考那里的博客。
但是在DP这里并不能直接使用上面的做法。因为我们在PEV阶段获得的是状态值函数而不是动作值函数。为此,利用动作值函数和状态值函数之间的关系 q π ( s , a ) = ∑ s ′ ∈ S P s s ′ a ( r s s ′ a + γ v π ( s ′ ) ) q^\pi(s,a)=\sum_{s'\in\mathcal{S}}\mathcal{P}_{ss'}^a (r_{ss'}^a +\gamma v^\pi(s')) qπ(s,a)=∑s′∈SPss′a(rss′a+γvπ(s′))可以对上述式子进行变形得到: π ′ ( s ) ← arg max π ′ { ∑ s ′ ∈ S P s s ′ a ( r + γ V ∞ π ( s ′ ) ) } . \pi^{\prime}(s)\leftarrow\arg\max_{\pi^{\prime}}\biggl\{\sum_{s^{\prime}\in\mathcal{S}}\mathcal{P}_{ss'}^a\left(r+\gamma V_{\infty}^{\pi}\bigl(s'\bigr)\right)\biggr\}. π′(s)←argπ′max{s′∈S∑Pss′a(r+γV∞π(s′))}.
5.2.3 收敛性证明
策略迭代的收敛性取决于以下两个部分:
- 每次PEV利用的self-consistency条件都是可解的(solvable):换句话说,即每个策略对应的状态值函数的值都是有限的。或者说,每个中间策略都是可行的(feasible),能够有效的和环境交互。而且因为策略迭代的算法使然,每次策略都是建立在之前策略的基础上。因此,就要求初始的策略是可行的。
- PIM每次获得的新策略是更优的:这个在上面已经证明过了。
下面分两部分给出上述两个部分收敛性的详细证明。
5.2.3.1 PEV的收敛性
回顾DP的PEV更新公式的矩阵形式:
X
=
γ
B
X
+
b
.
X=\gamma BX+b.
X=γBX+b. 可以看出,这个算子每次相当于对X做了一个变换后把变换后的值赋值给X。我们定义这个变换:
L
(
X
)
=
γ
B
X
+
b
,
\mathcal{L}(X)=\gamma BX+b,
L(X)=γBX+b, 我们下面要证明这个矩阵的收敛速度为
γ
\gamma
γ-contractive。证明如下:
∥
L
(
X
j
+
1
)
−
L
(
X
j
)
∥
∞
=
γ
∥
B
(
X
j
+
1
−
X
j
)
∥
∞
≤
γ
∥
B
max
i
∈
{
1
,
2
,
⋯
,
n
}
(
∣
X
j
+
1
i
−
X
j
i
∣
)
∥
∞
=
γ
∥
B
∥
X
j
+
1
−
X
j
∥
∞
∥
∞
=
γ
∥
B
∥
∞
∥
X
j
+
1
−
X
j
∥
∞
=
γ
∥
X
j
+
1
−
X
j
∥
∞
.
\begin{aligned}\left\|\mathcal{L}(X_{j+1})-\mathcal{L}(X_j)\right\|_\infty&=\gamma\|B(X_{j+1}-X_j)\|_\infty\\ &\leq\gamma\left\|B\max_{i\in\{1,2,\cdots,n\}}(\left|X_{j+1}^i-X_j^i\right|)\right\|_\infty\\ &=\gamma\left\|B\|X_{j+1}-X_j\|_\infty\right\|_\infty \\ &=\gamma\|B\|_\infty\|X_{j+1}-X_j\|_\infty\\ &=\gamma\|X_{j+1}-X_j\|_\infty. \end{aligned}
∥L(Xj+1)−L(Xj)∥∞=γ∥B(Xj+1−Xj)∥∞≤γ
Bi∈{1,2,⋯,n}max(
Xj+1i−Xji
)
∞=γ∥B∥Xj+1−Xj∥∞∥∞=γ∥B∥∞∥Xj+1−Xj∥∞=γ∥Xj+1−Xj∥∞. 从左边到右一式把
L
(
X
)
\mathcal{L}(X)
L(X)的定义式代入即可得。右一式到右二式根据矩阵乘法的意义结合了矩阵的无穷范数的定义可得,注意X是一个列向量。右二式到右三式根据矩阵的无穷范数的定义即可。右三式到右四式是矩阵的无穷范数的性质。右四式到右五式是因为上面已经证明过
∥
B
∥
∞
=
1
\|B\|_\infty=1
∥B∥∞=1。所以,我们得到最终的不等式:
∥
L
(
X
j
+
1
)
−
L
(
X
j
)
∥
∞
≤
γ
∥
X
j
+
1
−
X
j
∥
∞
.
\left\|\mathcal{L}(X_{j+1})-\mathcal{L}(X_j)\right\|_\infty\leq\gamma\|X_{j+1}-X_j\|_\infty.
∥L(Xj+1)−L(Xj)∥∞≤γ∥Xj+1−Xj∥∞. 接下来把
L
(
X
j
+
1
)
=
X
j
+
2
\mathcal{L}(X_{j+1})=X_{j+2}
L(Xj+1)=Xj+2以及
L
(
X
j
)
=
X
j
+
1
\mathcal{L}(X_j)=X_{j+1}
L(Xj)=Xj+1代入上式,就可以得到:
∥
X
j
+
2
−
X
j
+
1
∥
∞
≤
γ
∥
X
j
+
1
−
X
j
∥
∞
.
\left\|X_{j+2}-X_{j+1}\right\|_\infty\leq\gamma\|X_{j+1}-X_j\|_\infty.
∥Xj+2−Xj+1∥∞≤γ∥Xj+1−Xj∥∞. 这就证明了每次相邻两次迭代的值函数的更新量
Δ
X
\Delta X
ΔX是以指数速度收敛的。因此,
L
(
X
)
\mathcal{L}(X)
L(X)是一个
γ
\gamma
γ-contractive的operator。由上述证明以及PEV那里关于只有一个fixed-point的证明,我们可以很轻松得到每次PEV在迭代次数趋于无穷时,值函数的估计值
V
j
π
(
s
)
V_{j}^{\pi}(s)
Vjπ(s)收敛于真实值函数
v
π
(
s
)
v^{\pi}(s)
vπ(s)(即那个唯一的fixed-point)。
5.2.3.2 DP的策略迭代的收敛性
欲证明DP的策略迭代的收敛性,同样需要证明以下两个部分:
- 每次PEV得到的值函数序列{ V ∞ π 0 , V ∞ π 1 , ⋯ V_{\infty}^{\pi_0 },V_{\infty}^{\pi_1},\cdots V∞π0,V∞π1,⋯}是单调递增的
- 每次PIM得到的策略序列{ π 0 , π 1 , ⋯ \pi_0,\pi_1,\cdots π0,π1,⋯}最终收敛到最优策略 π ∗ \pi^* π∗
首先来证明第一点。即证明:
V
∞
π
0
(
s
)
≤
V
∞
π
1
(
s
)
≤
V
∞
π
2
(
s
)
≤
⋯
≤
V
∞
π
k
(
s
)
≤
V
∞
π
k
+
1
(
s
)
≤
⋯
≤
v
∗
.
V_\infty^{\pi_0}(s)\leq V_\infty^{\pi_1}(s)\leq V_\infty^{\pi_2}(s)\leq\cdots\leq V_\infty^{\pi_k}(s)\leq V_\infty^{\pi_{k+1}}(s)\leq\cdots\leq v^*.
V∞π0(s)≤V∞π1(s)≤V∞π2(s)≤⋯≤V∞πk(s)≤V∞πk+1(s)≤⋯≤v∗. 首先,根据我们在PEV收敛性那里的证明,我们知道每次PEV都会收敛到
v
π
k
(
s
)
v^{\pi_k}(s)
vπk(s),即:
V
∞
π
k
(
s
)
=
v
π
k
(
s
)
.
V_\infty^{\pi_k}(s)=v^{\pi_k}(s).
V∞πk(s)=vπk(s). 因为每次PIM更新后的策略都比原来更好,即:
v
π
k
(
s
)
≤
v
π
k
+
1
(
s
)
for
∀
s
∈
S
.
v^{\pi_k}(s)\leq v^{\pi_{k+1}}(s)\quad\text{for}\ \forall s\in S.
vπk(s)≤vπk+1(s)for ∀s∈S. 那么我们用大V(估计值)来替换小v(真实值)就有:
V
∞
π
k
(
s
)
≤
V
∞
π
k
+
1
(
s
)
.
V_\infty^{\pi_k}(s)\leq V_\infty^{\pi_{k+1}}(s).
V∞πk(s)≤V∞πk+1(s). 证毕。
接下来证明第二点。即证明:
When PIM stops improvement, i.e.,
π
∞
=
π
∗
if
π
∗
is unique
\text{When PIM stops improvement, i.e., }\pi_\infty=\pi^*\text{ if }\pi^*\text{ is unique}
When PIM stops improvement, i.e., π∞=π∗ if π∗ is unique 首先,需要澄清两个2概念。可能一开始看到这个式子的时候会感到奇怪,
π
∗
\pi^*
π∗和
π
∞
\pi_\infty
π∞不就是同一个东西吗?为啥还要证明啊?这样想就大错特错了。这里
π
∗
\pi^*
π∗指的是满足Bellman最优性方程的最优策略,而
π
∞
\pi_\infty
π∞指的是当PIM停止迭代的时候得到的策略。这两个策略是不一样的。那么,我们接下来就证明这两个策略是一样的。针对相邻的两个策略
π
k
\pi_k
πk和
π
k
+
1
\pi_{k+1}
πk+1对应的状态值函数,我们有:
v
π
k
+
1
(
s
)
=
max
a
∑
s
′
∈
S
P
s
s
′
a
(
r
s
s
′
a
+
γ
v
π
k
(
s
′
)
)
,
∀
s
∈
S
v^{\pi_{k+1}}(s)=\max_a\sum_{s^{\prime}\in\mathcal{S}}\mathcal{P}_{ss'}^a\big(r_{ss'}^a+\gamma v^{\pi_k}(s^{\prime})\big),\forall s\in\mathcal{S}
vπk+1(s)=amaxs′∈S∑Pss′a(rss′a+γvπk(s′)),∀s∈S 这里的max下面加了一个a表示这里a是可以自由变动的变量。同时,当PIM终止时,有
π
∞
+
1
=
π
∞
\pi_{\infty+1}=\pi_{\infty}
π∞+1=π∞ 综合上面两式就有:
v
π
∞
(
s
)
=
v
π
∞
+
1
(
s
)
=
max
a
∑
s
′
∈
S
P
(
r
+
γ
v
π
∞
(
s
′
)
)
,
∀
s
∈
S
,
v^{\pi_\infty}(s)=v^{\pi_{\infty+1}}(s)=\max_a\sum_{s^{\prime}\in\mathcal{S}}\mathcal{P}\big(r+\gamma v^{\pi_\infty}(s^{\prime})\big),\forall s\in\mathcal{S},
vπ∞(s)=vπ∞+1(s)=amaxs′∈S∑P(r+γvπ∞(s′)),∀s∈S, 即:
v
π
∞
(
s
)
=
max
a
∑
s
′
∈
S
P
(
r
+
γ
v
π
∞
(
s
′
)
)
,
∀
s
∈
S
,
v^{\pi_\infty}(s)=\max_a\sum_{s^{\prime}\in\mathcal{S}}\mathcal{P}\big(r+\gamma v^{\pi_\infty}(s^{\prime})\big),\forall s\in\mathcal{S},
vπ∞(s)=amaxs′∈S∑P(r+γvπ∞(s′)),∀s∈S, 而我们知道根据第一类Bellman方程解出来的最优策略
π
∗
\pi^*
π∗满足:
v
π
∗
(
s
)
=
max
a
∑
s
′
∈
S
P
(
r
+
γ
v
π
∗
(
s
′
)
)
,
∀
s
∈
S
.
v^{\pi^*}(s)=\max_a\sum_{s^{\prime}\in\mathcal{S}}\mathcal{P}\big(r+\gamma v^{\pi^*}(s^{\prime})\big),\forall s\in\mathcal{S}.
vπ∗(s)=amaxs′∈S∑P(r+γvπ∗(s′)),∀s∈S. 那么,根据上面的推导,我们就可以得到
π
∞
=
π
∗
\pi_\infty=\pi^*
π∞=π∗。证毕。
5.2.4 使用Newton-Raphson法来解释DP的策略迭代
可能大家会对于DP的策略迭代背后的原理感到好奇,例如为什么策略迭代为什么往往比值迭代更快的收敛,为什么通常更具有鲁棒性?这背后有没有什么数学原理呢?这里就提供一种使用最优化理论里的Newton-Raphson法来解释DP的策略迭代的原理
首先介绍一下Newton-Raphson法。它的迭代公式为:
X
k
+
1
=
X
k
−
[
∇
g
(
X
k
)
]
−
1
g
(
X
k
)
,
X_{k+1}=X_k-[\nabla g(X_k)]^{-1}g(X_k),
Xk+1=Xk−[∇g(Xk)]−1g(Xk),
g
(
X
)
g(X)
g(X)是一个关于向量X的向量函数,
∇
g
(
X
)
\nabla g(X)
∇g(X)就是
g
(
X
)
g(X)
g(X)的Jacobian矩阵,即:
∇
g
(
X
)
=
[
∂
g
1
(
X
)
∂
X
1
∂
g
1
(
X
)
∂
X
2
⋯
∂
g
1
(
X
)
∂
X
n
∂
g
2
(
X
)
∂
X
1
∂
g
2
(
X
)
∂
X
2
⋯
∂
g
2
(
X
)
∂
X
n
⋮
⋮
⋱
⋮
∂
g
m
(
X
)
∂
X
1
∂
g
m
(
X
)
∂
X
2
⋯
∂
g
m
(
X
)
∂
X
n
]
.
\nabla g(X)=\begin{bmatrix}\frac{\partial g_1(X)}{\partial X_1}&\frac{\partial g_1(X)}{\partial X_2}&\cdots&\frac{\partial g_1(X)}{\partial X_n}\\\frac{\partial g_2(X)}{\partial X_1}&\frac{\partial g_2(X)}{\partial X_2}&\cdots&\frac{\partial g_2(X)}{\partial X_n}\\\vdots&\vdots&\ddots&\vdots\\\frac{\partial g_m(X)}{\partial X_1}&\frac{\partial g_m(X)}{\partial X_2}&\cdots&\frac{\partial g_m(X)}{\partial X_n}\end{bmatrix}.
∇g(X)=
∂X1∂g1(X)∂X1∂g2(X)⋮∂X1∂gm(X)∂X2∂g1(X)∂X2∂g2(X)⋮∂X2∂gm(X)⋯⋯⋱⋯∂Xn∂g1(X)∂Xn∂g2(X)⋮∂Xn∂gm(X)
. 这里的
∇
g
(
X
)
\nabla g(X)
∇g(X)是一个m行n列的矩阵。
现在我们尝试使用Newton-Raphson法来解释DP的策略迭代。考虑一个简化的情况,其中状态空间
S
\mathcal{S}
S和动作空间
A
\mathcal{A}
A都各自只有两个元素。首先回顾一下第一类Bellman方程:
v
∗
(
s
)
=
max
a
∈
A
∑
s
′
∈
S
P
(
s
′
∣
s
,
a
)
(
r
s
s
′
a
+
γ
v
∗
(
s
′
)
)
,
∀
s
∈
S
.
v^*(s)=\max_{a\in\mathcal{A}}\sum_{s^{\prime}\in\mathcal{S}}\mathcal{P}(s^{\prime}|s,a)\left(r_{ss^{\prime}}^a+\gamma v^*(s^{\prime})\right),\forall s\in\mathcal{S}.
v∗(s)=a∈Amaxs′∈S∑P(s′∣s,a)(rss′a+γv∗(s′)),∀s∈S. 那么我们首先为两个可能的状态取值
s
(
1
)
s_{(1)}
s(1)和
s
(
2
)
s_{(2)}
s(2)写出它们的第一类Bellman方程:
V
π
(
s
(
1
)
)
=
max
a
∑
s
′
P
s
(
1
)
s
′
a
(
r
+
γ
V
π
(
s
′
)
)
V
π
(
s
(
2
)
)
=
max
a
∑
s
′
P
s
(
2
)
s
′
a
(
r
+
γ
V
π
(
s
′
)
)
,
a
∈
{
a
(
1
)
,
a
(
2
)
}
.
\begin{gathered} V^\pi(s_{(1)})=\max_a\sum_{s^{\prime}}\mathcal{P}_{s_{(1)}s^{\prime}}^a\big(r+\gamma V^\pi(s^{\prime})\big)\\ V^{\pi}\big(s_{(2)}\big)=\max_a\sum_{s^{\prime}}\mathcal{P}_{s_{(2)}s^{\prime}}^a\big(r+\gamma V^\pi(s^{\prime})\big),\\a\in\{a_{(1)},a_{(2)}\}. \end{gathered}
Vπ(s(1))=amaxs′∑Ps(1)s′a(r+γVπ(s′))Vπ(s(2))=amaxs′∑Ps(2)s′a(r+γVπ(s′)),a∈{a(1),a(2)}. 这里我们引入Bellman operator来做一个形式上的化简:
X
=
B
(
X
)
X=\mathcal{B}(X)
X=B(X) 这里
X
=
[
V
π
(
s
(
1
)
)
V
π
(
s
(
2
)
)
]
X=\begin{bmatrix}V^\pi(s_{(1)})\\V^\pi(s_{(2)})\end{bmatrix}
X=[Vπ(s(1))Vπ(s(2))],这里注意
B
\mathcal{B}
B是一个非线性的operator,因为在Bellman方程中包含max操作。为了进一步化简,可以引入函数
g
(
X
)
:
R
2
→
R
2
g(X){:}\mathbb{R}^2\to\mathbb{R}^2
g(X):R2→R2:
g
(
X
)
=
B
(
X
)
−
X
,
g(X)=\mathcal{B}(X)-X,
g(X)=B(X)−X,
g
(
X
)
g(X)
g(X)可以写成分量的形式:
g
(
X
)
=
[
g
1
(
X
)
g
2
(
X
)
]
g(X)=\begin{bmatrix}g_1(X)\\g_2(X)\end{bmatrix}
g(X)=[g1(X)g2(X)] 其中每个分量
g
i
(
X
)
g_i(X)
gi(X)的具体表达式如下:
g
i
(
X
)
=
max
j
(
γ
β
i
j
T
X
+
λ
i
j
)
−
X
i
,
i
,
j
∈
{
1
,
2
}
,
β
i
j
=
[
P
s
(
i
)
s
(
1
)
a
(
j
)
,
P
s
(
i
)
s
(
2
)
a
(
j
)
]
T
,
λ
i
j
=
∑
s
′
P
s
(
i
)
s
′
a
(
j
)
r
s
(
i
)
s
′
a
(
j
)
.
g_i(X)=\max_j\left(\gamma\beta_{ij}^\mathrm{T}X+\lambda_{ij}\right)-X_i ,i,j\in\{1,2\},\\\beta_{ij}=\left[\mathcal{P}_{s_{(i)}s_{(1)}}^{a_{(j)}},\mathcal{P}_{s_{(i)}s_{(2)}}^{a_{(j)}}\right]^\mathrm{T},\lambda_{ij}=\sum_{s'}\mathcal{P}_{s_{(i)}s'}^{a_{(j)}}r_{s_{(i)}s'}^{a_{(j)}} .
gi(X)=jmax(γβijTX+λij)−Xi,i,j∈{1,2},βij=[Ps(i)s(1)a(j),Ps(i)s(2)a(j)]T,λij=s′∑Ps(i)s′a(j)rs(i)s′a(j). 这里的
g
i
(
X
)
g_i(X)
gi(X)仍然包含max操作,但是经过我们观察可以知道,
g
i
(
X
)
g_i(X)
gi(X)是一个分成两段的各段均为线性函数的分段函数。
为此,可以把
g
i
(
X
)
g_i(X)
gi(X)写成:
g
i
(
X
)
=
{
γ
β
i
1
T
X
+
λ
i
1
−
X
i
,
f
o
r
X
∈
I
i
1
γ
β
i
2
T
X
+
λ
i
2
−
X
i
,
f
o
r
X
∈
I
i
2
,
g_i(X)=\begin{cases}\gamma\beta_{i1}^\mathrm{T}X+\lambda_{i1}-X_i,&\quad\mathrm{for}X\in\mathbb{I}_{i1}\\\gamma\beta_{i2}^\mathrm{T}X+\lambda_{i2}-X_i,&\quad\mathrm{for}X\in\mathbb{I}_{i2}\end{cases},
gi(X)={γβi1TX+λi1−Xi,γβi2TX+λi2−Xi,forX∈Ii1forX∈Ii2, 这里
I
i
1
\mathbb{I}_{i1}
Ii1和
I
i
2
\mathbb{I}_{i2}
Ii2分别表示使得
g
i
(
X
)
g_i(X)
gi(X)在分别取动作
a
(
1
)
a_{(1)}
a(1)和
a
(
2
)
a_{(2)}
a(2)达到max的
X
X
X的子区间。那么因为
g
(
X
)
g(X)
g(X)有两个分量
g
1
(
X
)
g_1(X)
g1(X)和
g
2
(
X
)
g_2(X)
g2(X),而每个分量都包含两个子区间,因此关于
g
(
X
)
g(X)
g(X)的两个分量分别取什么动作达到最值就有四种情况。将平面划分为四个区域:
但是,在每个子区域内部,
g
(
X
)
g(X)
g(X)的表达式都是唯一的。因此,四个子区域内部的表达式可以统一写为:
g
(
X
)
=
γ
β
σ
X
+
λ
σ
−
X
,
f
o
r
X
∈
U
σ
,
σ
∈
{
1
,
2
,
3
,
4
}
β
σ
=
[
P
s
(
1
)
s
(
1
)
a
∗
P
s
(
1
)
s
(
2
)
a
∗
P
s
(
2
)
s
(
1
)
a
∗
P
s
(
2
)
s
(
2
)
a
∗
]
,
λ
σ
=
[
∑
s
′
P
s
(
1
)
s
′
a
∗
r
s
(
1
)
s
′
a
∗
∑
s
′
P
s
(
2
)
s
′
a
∗
r
s
(
2
)
s
′
a
∗
]
,
\begin{aligned}&g(X)=\gamma\beta_\sigma X+\lambda_\sigma-X,\mathrm{for~}X\in\mathbb{U}_\sigma,\sigma\in\{1,2,3,4\}\\&\beta_\sigma=\begin{bmatrix}\mathcal{P}_{s_{(1)}s_{(1)}}^{a^*}&\mathcal{P}_{s_{(1)}s_{(2)}}^{a^*}\\\mathcal{P}_{s_{(2)}s_{(1)}}^{a^*}&\mathcal{P}_{s_{(2)}s_{(2)}}^{a^*}\end{bmatrix},\lambda_\sigma=\begin{bmatrix}\sum_{s^{\prime}}\mathcal{P}_{s_{(1)}s^{\prime}}^{a^*}r_{s_{(1)}s^{\prime}}^{a^*}\\\sum_{s^{\prime}}\mathcal{P}_{s_{(2)}s^{\prime}}^{a^*}r_{s_{(2)}s^{\prime}}^{a^*}\end{bmatrix},\end{aligned}
g(X)=γβσX+λσ−X,for X∈Uσ,σ∈{1,2,3,4}βσ=[Ps(1)s(1)a∗Ps(2)s(1)a∗Ps(1)s(2)a∗Ps(2)s(2)a∗],λσ=[∑s′Ps(1)s′a∗rs(1)s′a∗∑s′Ps(2)s′a∗rs(2)s′a∗], 其中,
a
∗
a^*
a∗表示在此时的区域内应该采取的最优动作(注意,第一行的
a
∗
a^*
a∗和第二行的
a
∗
a^*
a∗不一定一样)。这样写成标准的形式后,就方便后续的求导操作了。不过,还要补充说明的一点是单看上面的分段函数图像,在拐点处并不连续,不过这里随便取左边或右边段的导数值即可。有论文证明过这并不影响Newton-Raphson法的收敛性。那么,我们可以写出这四个区域内Jacobian矩阵的统一表达式:
∇
g
(
X
)
=
γ
β
σ
−
I
,
f
o
r
X
∈
U
σ
,
σ
=
1
,
2
,
3
,
4
\nabla g(X)=\gamma\beta_\sigma-I,\mathrm{for~}X\in\mathbb{U}_\sigma,\sigma=1,2,3,4
∇g(X)=γβσ−I,for X∈Uσ,σ=1,2,3,4 那么Newton-Raphson法的迭代公式就可以写成:
X
k
+
1
=
X
k
−
[
∇
g
(
X
)
]
−
1
g
(
X
k
)
=
X
k
−
(
γ
β
σ
−
I
)
−
1
(
γ
β
σ
X
k
+
λ
σ
−
X
k
)
=
−
(
γ
β
σ
−
I
)
−
1
λ
σ
.
\begin{aligned} X_{k+1}& =X_k-[\nabla g(X)]^{-1}g(X_k) \\ &=X_k-(\gamma\beta_\sigma-I)^{-1}(\gamma\beta_\sigma X_k+\lambda_\sigma-X_k) \\ &=-(\gamma\beta_\sigma-I)^{-1}\lambda_\sigma . \end{aligned}
Xk+1=Xk−[∇g(X)]−1g(Xk)=Xk−(γβσ−I)−1(γβσXk+λσ−Xk)=−(γβσ−I)−1λσ. 至此,我们就获得了从Newton-Raphson法的角度来得到的每轮迭代更新公式:
X
k
+
1
=
−
(
γ
β
σ
−
I
)
−
1
λ
σ
.
X_{k+1}=-(\gamma\beta_\sigma-I)^{-1}\lambda_\sigma .
Xk+1=−(γβσ−I)−1λσ.
下面我们从DP的策略迭代的角度来看看获得的迭代公式能不能于与上面刚得到的更新公式等价。在第k轮迭代时,经历过PEV后,可以得到
X
k
+
1
X_{k+1}
Xk+1的表达式:
X
k
+
1
=
γ
B
X
k
+
b
.
X_{k+1}=\gamma BX_k+b.
Xk+1=γBXk+b. 这里的
B
B
B和
b
b
b的具体表达式如下:
B
π
(
i
,
j
)
=
∑
a
∈
A
π
(
a
∣
s
(
i
)
)
P
(
s
(
j
)
∣
s
(
i
)
,
a
)
,
b
π
(
i
)
=
∑
a
∈
A
π
(
a
∣
s
(
i
)
)
∑
j
P
(
s
(
j
)
∣
s
(
i
)
,
a
)
r
s
(
i
)
s
(
j
)
a
.
B_{\pi_{(i,j)}}=\sum_{a\in\mathcal{A}}\pi\big(a\big|s_{(i)}\big)\mathcal{P}\big(s_{(j)}|s_{(i)},a\big),\\b_{\pi_{(i)}}=\sum_{a\in\mathcal{A}}\pi\big(a\big|s_{(i)}\big)\sum_j\mathcal{P}\big(s_{(j)}|s_{(i)},a\big)r_{s(i)s(j)}^a.
Bπ(i,j)=a∈A∑π(a
s(i))P(s(j)∣s(i),a),bπ(i)=a∈A∑π(a
s(i))j∑P(s(j)∣s(i),a)rs(i)s(j)a. 那么在第k轮的PIM中,更新策略其实就是寻找一个策略
π
k
+
1
\pi_{k+1}
πk+1使得:
π
k
+
1
=
arg
max
π
{
γ
B
π
X
k
+
b
π
}
.
\pi_{k+1}=\arg\max_\pi\{\gamma B_\pi X_k+b_\pi\}.
πk+1=argπmax{γBπXk+bπ}. 在进行第k+1轮的PEV时,
X
k
+
1
X_{k+1}
Xk+1应该还满足在策略
π
k
+
1
\pi_{k+1}
πk+1下的self-consistency条件:
X
k
+
1
=
γ
B
π
k
+
1
X
k
+
1
+
b
π
k
+
1
,
X_{k+1}=\gamma B_{\pi_{k+1}}X_{k+1}+b_{\pi_{k+1}},
Xk+1=γBπk+1Xk+1+bπk+1, 这个方程的解为(从策略迭代角度得到的更新公式):
X
k
+
1
=
(
I
−
γ
B
π
k
+
1
)
−
1
b
π
k
+
1
.
X_{k+1}=\left(I-\gamma B_{\pi_{k+1}}\right)^{-1}b_{\pi_{k+1}}.
Xk+1=(I−γBπk+1)−1bπk+1. 如果策略采取确定性的greedy策略,那么
B
B
B和
b
b
b可以简化为:
B
i
j
∗
=
∑
a
∈
A
π
∗
(
a
∣
s
(
i
)
)
P
(
s
(
j
)
∣
s
(
i
)
,
a
)
=
P
s
(
i
)
s
(
j
)
a
∗
,
b
i
∗
=
∑
a
∈
A
π
∗
(
a
∣
s
(
i
)
)
∑
j
P
(
s
(
j
)
∣
s
(
i
)
,
a
)
r
=
∑
s
′
∈
S
P
s
(
i
)
s
′
a
∗
r
s
(
i
)
s
′
a
∗
.
B_{ij}^{*}=\sum_{a\in\mathcal{A}}\pi^{*}\big(a\big|s_{(i)}\big)\mathcal{P}\big(s_{(j)}|s_{(i)},a\big)=\mathcal{P}_{s_{(i)}s_{(j)}}^{a^{*}} ,\\b_{i}^{*}=\sum_{a\in\mathcal{A}}\pi^{*}\big(a\big|s_{(i)}\big)\sum_{j}\mathcal{P}\big(s_{(j)}\big|s_{(i)},a\big)r=\sum_{s^{\prime}\in\mathcal{S}}\mathcal{P}_{s_{(i)}s^{\prime}}^{a^{*}}r_{s_{(i)}s^{\prime}}^{a^{*}} .
Bij∗=a∈A∑π∗(a
s(i))P(s(j)∣s(i),a)=Ps(i)s(j)a∗,bi∗=a∈A∑π∗(a
s(i))j∑P(s(j)
s(i),a)r=s′∈S∑Ps(i)s′a∗rs(i)s′a∗. 其中,
a
∗
a^*
a∗是最优策略(值函数最大)。比较
B
π
k
+
1
B_{\pi_{k+1}}
Bπk+1和
β
σ
\beta_\sigma
βσ,
b
π
k
+
1
b_{\pi_{k+1}}
bπk+1和
λ
σ
\lambda_\sigma
λσ,我们可以发现这两个矩阵和向量是一样的。因此,从策略迭代角度得到的更新公式也可以写成:
X
k
+
1
=
(
I
−
γ
B
π
k
+
1
)
−
1
b
π
k
+
1
=
−
(
γ
β
σ
−
I
)
−
1
λ
σ
.
X_{k+1}=\left(I-\gamma B_{\pi_{k+1}}\right)^{-1}b_{\pi_{k+1}}=-(\gamma\beta_\sigma-I)^{-1}\lambda_\sigma .
Xk+1=(I−γBπk+1)−1bπk+1=−(γβσ−I)−1λσ. 这恰与Newton-Raphson法的更新公式一致!因此,我们就证明了策略迭代的两步cycle与Newton-Raphson法是等价的。尽管上面的分析是在一个简化的情况下进行的,但是这个结论是可以推广到一般的情况的。同时,由于Newton-Raphson法具有二阶收敛性,因此策略迭代也具有二阶收敛性。这也解释了为什么策略迭代往往比值迭代更快的收敛,更具有鲁棒性。