强化学习实战[updating]

我一直希望找到一篇长文,能从实用的角度讲清楚强化学习。一年过去了,我决定自己动手写这样一篇文章。

编程语言

本文将使用Julia作为主要的编程语言来实现Reinforcement Learning中的大多数算法。

Why Julia?事实上,我更乐意使用Clojure,然而其基础工具库实在太匮乏,不得不放弃。Why not Python?我只是单纯觉得Julia写出来的代码更好看!不过阅读本文并不需要读者熟练掌握Julia,只需要有基本的编程思想即可,我会尽量将Julia这门语言独有的内容降到最低(必要的时候会给出解释和参考文献)。

1. Tabular Methods

“昨夜西风凋碧树,独上高楼,望尽天涯路。”

1.1 Day 1: Where's Eve?

Eve_Example

Eve_Example

假设我们的机器人Eve落在了一条数轴的原点处,它的任务是寻找到生命,在-3和5处分别有一盆花和一只小狗,Eve每天只能选择向左或向右移动一个单位,只要发现生命(花或者狗)后,任务便结束,然后Eve返回并获得相应的奖励(假设奖励分别是3和5,其它位置没有奖励)。现在先给出如下定义:

  • \(\mathcal{A}\):Actions, 动作空间,在这里包含两种可能{:left, :right}
  • \(a_t\): 在\(t\)时刻采取的行动
  • \(\mathcal{S}\): States, 状态空间,这里就是Eve所有可能的位置{-3, -2, -1, 0, 1, 2, 3, 4, 5}
  • \(s_t\): 在\(t\)时刻所处的状态
  • \(\mathcal{R}\):所有可能的奖励,这里由三个离散值{0, 3, 5}构成
  • \(R_t\):每一天获得的奖励,特别地,我们将游戏结束时的时间记为\(T\),任务结束时获得的奖励为\(R_T\)

假设Eve很笨,它的记忆只有一天(),每天只会等概率地随机选择向左或向右移动,将其记作策略\(\pi\)(尽管目前还只是随机游走,谈不上什么策略)。我们先观察以下两个指标:

  1. \(\bar{R_T} = \frac{\sum_{i=1}^{N} R_T^i}{N}\),N次试验的平均收益(3.746) 任务结束时获得的奖励分布
  2. \(\bar{T} = \frac{\sum_{i=1}^{N} T^i}{N}\),N次试验平均行动的次数(16.164) 任务结束时,行动次数的分布
using Lazy
using Plots
using StatsBase
using Match
using LaTeXStrings

gr()
theme(:dark)

# ### Day 1 ###
INIT_STATE = 4
REWARD = [3, 0, 0, 0, 0, 0, 0, 0, 5]
ACTIONS = [-1, 1]
STATES = 1:length(REWARD)

function env(init_state, policy, reward_calculator, state_transformer)
    s = init_state
    function play() 
        sₜ  = s
        aₜ = policy(sₜ)
        rₜ = reward_calculator(sₜ, aₜ)
        s = state_transformer(sₜ, aₜ)
        (sₜ, aₜ, rₜ)  # TODO: should return a NamedTuple instead in Julia-0.7
    end
end

# get intermediate State-Action-Reward sequence
get_SAR_seq(play) = @>> begin
    play
    repeatedly 
    takeuntil(res -> ((s, a, r) = res; r > 0))  # TODO: use NamedTuple here
end

random_policy(s) = rand(ACTIONS)
get_reward(s, a) = REWARD[s]
state_transform(s, a) = @match (s, a) begin 
    (1, -1) || (9, 1) => s
    _ => s + a
end

begin
    init_play() = env(INIT_STATE, random_policy, get_reward, state_transform)
    get_steps_count() = @> init_play() get_SAR_seq length
    get_final_reward() = @> init_play() get_SAR_seq last last  # TODO: use NamedTuple here
    N = 10000
    steps_samples = convert(Array{Int}, collect(repeatedly(N, get_steps_count)))
    rewards_samples = convert(Array{Int}, collect(repeatedly(N, get_final_reward)))

    bar(counts(steps_samples, maximum(steps_samples)))
    savefig("eve_steps_count.png")
    bar(counts(rewards_samples, maximum(rewards_samples)))
    savefig("eve_reward_count.png")
    println("""Mean of steps: $(mean(steps_samples));
            Mean of rewards: $(mean(rewards_samples))""")
end
# Mean of steps: 15.9168;
# Mean of rewards: 3.78

平均收益只有3.7左右,离最优解(5)还有点距离,此外平均实验次数居然到了16次。接下来,我们将一步步放松约束条件,改进Eve的行动策略。

这里先将前面Eve与环境交互的过程用下图抽象出来:

MDP的序列化描述

MDP的序列化描述

每一步的reward只由当前时刻的state和action共同决定(这里对于Eve来说,每天的reward其实只与state相关,环境并不受Eve的行动影响,当然,真实情况要比这复杂得多,我们一步步来),而下一时刻的state则只受上一个时刻的state和action共同决定。

1.1 Day 2: The Value of State

某一天,Eve感到很迷茫,它不知道未来这么走下去,收益究竟有多少,我们可以给Eve算一下未来收益的总和\(G\)

\[ \begin{equation} G_t = R_t + R_{t+1} + R_{t+2} ... \end{equation} \]

不过,Eve可能并不这么认为,同样的收益,第5天得到还是第3天得到对于Eve来说有着不同的意义(显然后者的意义更大),所以不妨给每天的收益增加一个基于时间的折扣系数\(\gamma\)(这么做有许多好处,数学上处理起来更方便,而且对于infinite的问题也更容易求解,当然,根据实际问题不同,你完全可以设计不同的\(G_t\)计算方式)。

\[ \begin{equation} G_t = R_t + \gamma R_{t+1} + \gamma^2 R_{t+2} + ... \end{equation} \]

那么,根据Eve在\(t\)时刻所处的状态不同,假设可以算出收益的期望,称作value function:

\[ \begin{equation} V_t(s) = \mathbb{E}(G_t | s) = \mathbb{E}(R_t + \gamma R_{t+1} + \gamma^2 R_{t+2} + ... | s) \label{value_equation} \end{equation} \] 从上帝视角来看,在环境不发生变化的情况下,显然每个状态s都有一个对应的固定的\(V(s)\),这样,每次Eve决定怎么走的时候,只需要看下这张表,从可以到达的所有状态中,找到价值最高的状态\(s'\),跳转过去即可。现在问题变成了如何计算\(V(s)\)。我们不妨先模拟下,将Eve放在不同的位置,然后分别执行原来的随机策略,计算平均收益\(\hat{G_t}\),来估计\(V(s)\)

# ### Day 2 ###
calc_gain(s, γ) = begin
    init_play() = env(s, random_policy, get_reward, state_transform)
    rewards = @>> init_play() get_SAR_seq map(res -> last(res))
    r, n = reduce((temp, r) -> ((g, n) = temp; (g + γ^n * r, n + 1)),  # TODO: use Argument destructuring here
                  (0, 0),
                  rewards)
    r
end

begin
    params = [0.95, 0.98, 1.0]
    plot(layout=@layout([a;b;c]))
    for i in 1:length(params)
        γ = params[i]
        V = zeros(length(REWARD));
        for s in STATES
            V[s] = mean(repeatedly(1000, () -> calc_gain(s, γ)))
        end
        println(V)
        bar!(V, subplot=i, label=latexstring("\\gamma = $(γ)"))
    end
    savefig("eve_v_estimate.png")
end
# [3.0, 2.40616, 2.07831, 2.00239, 2.03833, 2.32305, 2.91446, 3.70394, 5.0]
# [3.0, 2.83945, 2.75356, 2.76483, 2.93804, 3.24011, 3.64397, 4.20483, 5.0]
# [3.0, 3.254, 3.498, 3.736, 3.99, 4.252, 4.558, 4.772, 5.0]
随机模拟估计V(s),gamma=1.0

随机模拟估计V(s),gamma=1.0

看起来上图似乎符合我们的直观感受。\(\gamma\)越大,我们就越是看重长期收益,反之则更看重短期收益。

不过,对我们的机器人Eve来说,手上并没有这样一张表,前面虽然通过模拟得到了这样一张表,但是效率还很低。那是否能找到解析解呢?首先引入状态转移的概念,这里Eve每一时刻所处的状态只与其上一时候所处的状态(和行动)有关,不受更早时候状态(和行为)的影响(即具有Markov性质),于是有:

\[ \begin{equation} P(s'|s) = \sum_{a \in A} P(s' | s, a) \end{equation} \]

对于前面Eve机器人随机策略而言,对应的概率转移矩阵为:

\[ \begin{equation} \begin{bmatrix} 1 &0 &0 &0 &0 &0 &0 &0 &0 \\ 0.5 &0 & 0.5 &0 &0 &0 &0 &0 &0 \\ 0 &0.5 &0 & 0.5 &0 &0 &0 &0 &0 \\ 0 &0 &0.5 &0 & 0.5 &0 &0 &0 &0 \\ 0 &0 &0 &0.5 &0 & 0.5 &0 &0 &0 \\ 0 &0 &0 &0 &0.5 &0 & 0.5 &0 &0 \\ 0 &0 &0 &0 &0 &0.5 &0 & 0.5 &0 \\ 0 &0 &0 &0 &0 &0 &0.5 &0 & 0.5 \\ 0 &0 &0 &0 &0 &0 &0 &0 & 1 \\ \end{bmatrix} \end{equation} \]

于是,\(\eqref{value_equation}\)可以写成:

\[ \begin{equation} V_t(s) = R_t(s) + \gamma \sum_{s' \in S} P(s' | s) V_{t+1}(s') \end{equation} \]

上式构建了不同时刻的价值函数之间的关系,在这里由于状态转移函数是固定的,最终会达到一个稳态(TODO: 这里需要展开说下,MRP),即:

\[ \begin{equation} \begin{bmatrix} V^*(s_1) \\ \dots \\ V^*(s_N) \end{bmatrix} = \begin{bmatrix} R(s_1) \\ \dots \\ R(s_N) \end{bmatrix} + \gamma \begin{bmatrix} P(s_1|s_1) & \dots &P(s_N|s_1)\\ \vdots &\ddots &\vdots \\ P(s_1|s_N) &\dots &P(s_N|s_N) \end{bmatrix} \begin{bmatrix} V^*(s_1) \\ \dots \\ V^*(s_N) \end{bmatrix} \end{equation} \]

采用矩阵的表示方法即是:

\[ \begin{equation} \begin{split} \boldsymbol{V} &= \boldsymbol{R} + \gamma \boldsymbol{P V}\\ \boldsymbol{V} &= (\boldsymbol{I} - \gamma \boldsymbol{P})^{-1} \boldsymbol{R} \end{split} \end{equation} \]

稠密矩阵求逆的复杂度为\(O(n^3)\)(TODO: 这里需要插入相关文献),接下来,我们将分别采用Value Iteration和 Policy Iteration来求解。

Value Iteration

让我们先从上帝视角来看看这个问题。我们的本意是希望得到一张表,让Eve每次行动的时候,能够查看它从当前位置往左和往右走的期望价值,然后从中选出期望价值最大的行动\(a\)作为决策Policy的结果。显然,最糟糕的情况下,我们需要遍历所有的可能\(|A|^{|S|}\),不过通过Value Iteration,我们可以先找到稳态的\(V(s)\),然后直接得出最优的Policy。可以证明在最优策略下(TODO: Reference Needed Here, Bellman Function),价值函数满足以下式子:

\[ \begin{equation} V^*(s) = \underset{a \in A}{max} \left( R(s, a) + \gamma \sum_{s' \in S}P(s' | s,a) V^*(s') \right) \end{equation} \]

于是,我们可以先随机初始化\(V(s)\),然后根据上式迭代计算并更新,直至收敛。下图动态地展示了该过程:

# Value Iteration

function gen_value_iteration(γ)
    function value_iteration(values)
        function calc_V(s)
            maximum(map(a -> get_reward(s, a) + γ * values[state_transform(s, a)],
                        ACTIONS))
        end
        map(calc_V, STATES)
    end
end

function find_stable_value(γ)
    is_terminate(Vₜ, Vₜ₊₁) = maximum(abs.(Vₜ₊₁ - Vₜ)) < 0.1;
    value_iter_fn = gen_value_iteration(γ);
    gen_pair(seq) = zip(seq, tail(seq));
    init_values = zeros(Int, length(STATES))
    
    @>> iterate(value_iter_fn, init_values) begin
        gen_pair
        filter(pair -> is_terminate(pair...))
        first
        first
    end
end

begin
    params = [0.90, 0.95, 0.99]
    plot(layout=@layout([a;b;c]))
    for i in 1:length(params)
        γ = params[i]
        V = find_stable_value(γ)
        bar!(V, subplot=i, label=latexstring("\\gamma = $γ"))
    end
    savefig("eve_v_iteration.png")
end
gamma=0.9时,Value 迭代的过程

gamma=0.9时,Value 迭代的过程

不同gamma对应的稳态Value:

不同gamma对应的稳态Value

不同gamma对应的稳态Value

在得到稳态的\(V^*\)之后,即可根据下式得到对应的最优Policy:

\[ \begin{equation} \pi(s) \leftarrow \underset {a \in A} {arg \ max} \left( R(s, a) + \gamma \sum_{s' \in S}P(s' | s,a) V^*(s') \right) \end{equation} \]

Policy Iteration

Policy Iteration的思想是,先将所有状态对应的\(V\)置为0(也可以根据先验设置相应的值),针对每个状态,先随机初始化一个action(或根据前面设置的先验调整),记为\(\pi_0\),在给定Policy之后,我们可以算出下一步所有状态对应的\(V(s')\),即(Policy Evaluation):

\[ \begin{equation} V^{\pi}_t (s) = R(s, \pi(s)) + \gamma \sum_{s' \in S} P(s'|s, \pi(s)) V^{\pi}_{t+1}(s') \end{equation} \]

然后根据新得到的\(V_{t+1}(s)\),重新调整Policy,使得每个state对应的action都是朝着\(V(s'|s)\)最大的方向在移动,即(Policy update):

\[ \begin{equation} \pi_{k+1}(s) = \underset {a \in A} {arg \ max} \left( R(s,a) + \gamma \sum_{s' \in S} P(s' | s, a) V^{\pi_k}(s') \right) \end{equation} \]

如此迭代,直至Policy不再变化。如何证明最优呢?(TODO: 单调性)

# Policy Iteration

function policy_evaluation(V, π, γ)
    get_value(s, a) = @match (s, a) begin 
        (1, -1)  => V[1]
        (9, 1) => V[9]
        _ => V[s + a]
    end
    
    reward = [get_reward(sa...) for sa in zip(STATES, π)]
    future_reward = γ * map(s -> get_value(s, π[s]), STATES)
    reward + future_reward
end

function policy_update(V)
    π = [] 
    push!(π, V[1] < V[2] ? 1 : -1)
    for i in 2:8
        push!(π, (V[i - 1] < V[i + 1]) ? 1 : -1)
    end
    push!(π, V[end] < V[end - 1] ? -1 : 1)  # ensure go right at position 5
end

function policy_iter(PV, γ)  # TODO: use argument dispatch here
    πₜ, Vₜ = PV
    println(PV)
    Vₜ₊₁ = policy_evaluation(Vₜ, πₜ, γ)
    πₜ₊₁ = policy_update(Vₜ₊₁)
    πₜ₊₁, Vₜ₊₁ 
end

begin
    V₀ = zeros(Float32, length(STATES))
    π₀ = [random_policy(s) for s in STATES]
    γ = 0.9  # try to change this value
    temp_PV = @>> iterate(PV -> policy_iter(PV, γ), (π₀, V₀)) take(20)
    plot()
    @gif for PV in temp_PV
        π, V = PV
        bar(V, label="V")
        for pair in zip(π, 1:9)
            dir, i = pair
            if dir == -1
                plot!([i, i - 0.5], [0,0], arrow=10, linewidth=5, color = :orange, label = "")
            else
                plot!([i, i + 0.5], [0,0], arrow=10, linewidth=5, color = :blue, label = "")
            end
        end
    end
end
Policy Iteration

Policy Iteration

2. Approximate Methods

“衣带渐宽终不悔,为伊消得人憔悴。”

3. Deep Reinforcement Learning

“众里寻他千百度,蓦然回首,那人却在灯火阑珊处。”

参考文献

Books

  1. Reinforcement Learning: An Introduction. Second edition
  2. Decision Making Under Uncertainty Theory and Application
  3. Artificial Intelligence: Foundations of Computational Agents, 2nd Edition
  4. Approximate Dynamic Programming: Solving the Curses of Dimensionality, 2nd Edition
  5. Dynamic Programming and Optimal Control
  6. Reinforcement Learning State of the Art
  7. reinforcement learning and dynamic programming using function approximators