【強化学習】複数の探索ポリシーを実装/解説して比較してみた(多腕バンディット問題)

今回は探索アルゴリズムの話をします。
DQNシリーズでは ε-greedy しか基本使っていないですが、方策を決めるアルゴリズムは他にもたくさんあります。
多腕バンディット問題を元に、まずはQ学習で実装しその後にDQNに応用してみたいと思います。

本シリーズ

• 【強化学習】OpenAI Gym×Keras-rlで強化学習アルゴリズムを実装していくぞ（準備編）
• 【強化学習】OpenAI Gym×Keras-rlで強化学習アルゴリズムを実装していくぞ（Q学習編）
• 【強化学習】OpenAI Gym×Keras-rlで強化学習アルゴリズムを実装していくぞ（DQN編）
• 【強化学習】Keras-rlでRainbowを実装/解説
• 【強化学習】Keras-rlでApe-Xを実装/解説（並列プログラミング）
• MNISTでSaliencyMapとGrad_CAMを試してDQN(Dueling Network)を可視化してみた(実装/解説)
• 【強化学習】2018年度最強と噂のR2D2を実装/解説してみた
• 【強化学習】複数の探索ポリシーを実装/解説して比較してみた(多腕バンディット問題)(ここ)
• 【強化学習】DQNのハイパーパラメータを3つのゲームで比較してみた
• 【強化学習】R2D2を実装/解説してみたリベンジ 解説編（Keras-RL）
• 【強化学習】R2D3を実装/解説してみた（Keras-RL）
• 【強化学習】ついに人間を超えた！？Agent57を解説/実装してみた（Keras-RL）

概要

以下のアルゴリズムを実装した。

• ε-greedy
• 焼きなましε-greedy
• Softmax
• UCB1
• UCB1-Tuned
• UCB-V
• KL-UCB
• Thompson Sampling(ベータ分布)
• Thompson Sampling(正規分布)

コード全体

本記事で作成したコードは以下です。
※1ファイル完結です。
※GoogleColaboratoryは実行結果付き

• Q学習版コード(GitHubGist)
• Q学習版コード(GoogleColaboratory)
• DQN版コード(GitHubGist)
• DQN版コード(GoogleColaboratory)

多腕バンディット問題(Multi-armed bandit problem)について

簡単に言うと探索と活用のジレンマを定式化した問題です。

例として、３台(A,B,C)のスロットを考えます。
各スロットは以下の確率で当たりがでます。
A:10%
B:90%
C:50%
ただし、それぞれのスロットの当たりが出る確率をあなたは知りません。
ここでスロットのアームを100回回せる場合にどのアームを選べば貰える報酬が最大となるでしょうか？
という問題です。

スロットの確率を知るためにはそのスロットを回さないといけないが、その為には低い確率のスロットを回す可能性があり損してしまう事から、探索と活用のジレンマとなっています。

これは強化学習にも言えることで、報酬は(直近の)Q値の高いアクションで、探索は未知の空間にある高い報酬となります。

参考
・Multi-armed bandit(wikipediaが英語しかなかった…)
・多腕バンディット問題の理論とアルゴリズム
・強化学習入門：多腕バンディット問題
・バンディットアルゴリズムで最適な介入を見つける（基本編）
・バンディットアルゴリズム　基本編

強化学習の記号

記号をまとめておきます。

記号	強化学習	(多腕バンディットだと…?)
$t$	現在	試行回数
$k$	-	スロット
$K$	-	全スロットの集合
$s_t$	現在の状態	-
$s_{t+1}$	次の状態	-
$a_{t}$	アクション	どのスロットか($\in K$)
$r_{t}$	報酬	-
$r_{t}^k$	-	kのスロットの報酬
$Q(s_{t},a_{t})$	Q値($(s_{t},a_t)$での評価値)	-
$\mu(k,t)$	-	ある地点でのスロットの評価値

実装の前準備

今回複数の探索ポリシーを実装するので、メモリと同様にコマンドパターン化します。

Q学習編のコードを元に以下の変更を加えます。

QLAgent.py

class MyQLAgent(Agent):
    def __init__(self,
        action_policy,
        (省略)
        **kwargs):
        super(MyQLAgent, self).__init__(**kwargs)

        self.action_policy = action_policy
        (省略)

    # save/load が必要なポリシーがあるので追加
    def save_weights(self, filepath, overwrite=False):
        if overwrite or not os.path.isfile(filepath):
            d = {
                "q_table": self.q_table
            }
            # 保存したい内容を d に書き込む
            self.action_policy.save_weights(d)

            # 保存
            with open(filepath, 'wb') as f:
                pickle.dump(d, f)

    def load_weights(self, filepath):
        with open(filepath, 'rb') as f:
            d = pickle.load(f)
        self.q_table = d["q_table"]
        # 読み込み
        self.action_policy.load_weights(d)

    def forward(self, observation):
        (省略)

        # 行動を決定
        if self.training:

            (observationは文字列化済み)
            action = self.action_policy.select_action(self, observation)

        else:
            # Q値が最大のもの
            action = np.argmax(self.q_table[observation])

        self.prev_observation  = observation
        self.prev_action = action
        return action
    
    （省略）

class ActionPolicyClass():
    def __init__(self, 別途必要な引数):
        初期化

    # 保存例
    def save_weights(self, d):
        d["data"] = self.data

    # load例
    def load_weights(self, d):
        self.data = d["data"]

    def select_action(self, agent, state):
        アクションをきめて返す
        return action

ε-greedy の実装/解説

まずは一番簡単な ε-greedy です。
Q学習編で解説していますがもう一度解説します。

アルゴリズムは、乱数(0.0～1.0)に対して $\epsilon$ 以下ならランダムに行動、
それより大きければ評価値が最大となる行動(スロット)を選びます。

要するに普段は報酬が最大のスロットを選ぶが、ある一定の確率で探索をしようという方策です。
乱数の影響にすごく左右されるアルゴリズムとなります。

以下、実装です。

EpsilonGreedy.py

class EpsilonGreedy():
    def __init__(self, epsilon):
        self.epsilon = epsilon

    def save_weights(self, d):
        pass
    
    def load_weights(self, d):
        pass

    def select_action(self, agent, state):
        if self.epsilon > np.random.uniform(0, 1):
            # アクションをランダムに選択
            action = np.random.randint(0, agent.nb_actions)
        else:
            # 評価が最大のアクションを選択
            action = np.argmax(agent.q_table[state])
        return action

焼きなまし ε-greedy 法(Annealing ε-greedy)

焼きなましはこの記事で勝手に言っているだけで正式な名称ではありません。
DQNで提案された手法ですね。

考え方は、序盤はどのアームも確率が分からないので探索を優先し、後半は確率がある程度分かってくるので報酬を優先するというものです。
焼きなまし ε-greedy 法では $\epsilon$ の値を最初は高くし、後半になるほど低くすることで実現しています。

以下実装です。

AnnealingEpsilonGreedy.py

class AnnealingEpsilonGreedy():
    def __init__(self,  
        initial_epsilon,   # 初期ε
        final_epsilon,     # 最終状態でのε
        exploration_steps  # 初期～最終状態になるまでのステップ数
        ):
        self.epsilon_step = (initial_epsilon - final_epsilon) / exploration_steps
        self.initial_epsilon = initial_epsilon
        self.final_epsilon = final_epsilon
        self.step = 0

    def save_weights(self, d):
        d["step"] = self.step
    
    def load_weights(self, d):
        self.step = d["step"]

    def select_action(self, agent, state):

        # epsilon の計算
        epsilon = self.initial_epsilon - self.step * self.epsilon_step
        if epsilon < self.final_epsilon:
            epsilon = self.final_epsilon
        self.step += 1

        if epsilon > np.random.uniform(0, 1):
            # アクションをランダムに選択
            action = np.random.randint(0, agent.nb_actions)
        else:
            # 評価が最大のアクションを選択
            action = np.argmax(agent.q_table[state])
        return action

Softmax

画像分類の最後の活性化関数でも使われているSoftmax関数でアクションを決める方法です。

$$ y_i = \frac{e^{x_i}}{\sum_{j=1}^N e^{x_j}} $$

Softmax関数は合計が1になるように正規化する関数です。
評価が高いアーム程、高確率で選ばれるようになるアルゴリズムとなります。

参考
・softmax関数を直感的に理解したい
・Pythonでニューラルネットワークの活性化関数softmax関数を実装

実装は、Softmax の値に対して乱数を出し、その乱数が最大のアームを選びます。

SoftmaxPolicy.py

class SoftmaxPolicy():
    def __init__(self):
        pass

    def save_weights(self, d):
        pass
    
    def load_weights(self, d):
        pass

    def select_action(self, agent, state):

        exp_x = np.exp(agent.q_table[state])
        total = sum(exp_x)

        vals = []
        for i in range(agent.nb_actions):
            softmax = exp_x[i] / total

            # softmax 値以下の乱数を生成
            vals.append( random.uniform(0, softmax) )

        # 乱数の結果一番大きいアクションを選択
        action = np.argmax(vals)
        return action

UCB(Upper Confidence Bound)1

ε-greedy ではアームを回した回数を考慮できていませんでした。
そこで、UCBでは選択回数が少ないアーム=正しく評価されていないアームとして、
そのバランスを計算してアームを選ぶアルゴリズムとなります。

計算式は以下です。
各アクションに対するUCB値を出したら最大のUCB値を選択します。

UCB_a= \bar{x}_a + \sqrt{ \frac{2\log{N} }{n_a} }

記号	強化学習	(多腕バンディットだと…?)
$\bar{x}_a$	取得した報酬の平均(Q値/選択回数)	アームの報酬平均
$N$	$s_t$における$a_t$の選択回数の合計	全アームの選択回数
$n_a$	$s_t$における$a_t$の選択回数	アームの選択回数

以下、実装です。

UCB1.py

class UCB1():
    def __init__(self):
        # 選択回数を数える為の変数を追加
        self.ucb_count = {}

    def save_weights(self, d):
        d["ucb_count"] = self.ucb_count
    
    def load_weights(self, d):
        self.ucb_count = d["ucb_count"]

    def select_action(self, agent, state):

        # 状態がなければ新規追加
        if state not in self.ucb_count:
            # 0を回避するために最低を1で初期化
            self.ucb_count[state] = [1 for _ in range(agent.nb_actions)]

        # 合計を出す(数式ではN)
        total = sum(self.ucb_count[state])
        total = math.log(total)  # 重そうな計算なので for の外で

        # 各アクションのUCB値を計算
        ucbs = []
        for i in range(agent.nb_actions):
            n = self.ucb_count[state][i]

            # 現在のQ値の平均(数式だとx_a)を出す
            x = agent.q_table[state][i] / n

            # 数式を計算
            tmp = x + math.sqrt(2 * total / n )
            ucbs.append(tmp)

        # ucbが最大値となるアクションを選択
        action = np.argmax(ucbs)

        # 選択したアクションの選択回数を増やす
        self.ucb_count[state][action] += 1

        return action

UCBには派生がたくさんあるので有名どころ?をいくつか解説/実装します。

参考
・Finite-time Analysis of the Multiarmed Bandit
Problem*(論文)

UCB1-Tuned

UCB1-Tuned は、UCB1 に分散も考慮して改良したアルゴリズムです。
UCB1より優れた結果を出しますが理論的な保証はないとの事です。

UCBturned_a= \bar{x_a}+\sqrt{\frac{\log{N}}{n_a} \min( \frac{1}{4}, V_i(n_a) ) }

V_i({n_a}) = \hat{x}_a+\sqrt{\frac{2\log{N}}{n_a} }

$\min(a,b)$ は値の低い方を採用する関数、$\hat{x}_a$は取得した報酬の分散を表します。

以下実装です。

UCB1_Tuned.py

class UCB1_Tuned():
    def __init__(self):
        # 選択回数を数える為の変数を追加
        self.ucb_count = {}
        self.ucb_var = {}  # 分散用

    def save_weights(self, d):
        d["ucb_count"] = self.ucb_count
        d["ucb_var"] = self.ucb_var
    
    def load_weights(self, d):
        self.ucb_count = d["ucb_count"]
        self.ucb_var = d["ucb_var"]

    def select_action(self, agent, state):

        # 前回の報酬を計算
        state0 = agent.prev_observation  # 前回の状態
        action = agent.prev_action  # 前回のaction
        reward = agent.prev_reward  # 前回の報酬
        if state0 not in self.ucb_count:
            self.ucb_count[state0] = [1 for _ in range(agent.nb_actions)]
            self.ucb_var[state0] = [0 for _ in range(agent.nb_actions)]

        # 分散を更新
        prev_count = self.ucb_count[state0][action]
        prev_ave = agent.q_table[state0][action] / prev_count
        var = self.ucb_var[state0][action]
        var += ((reward - prev_ave) ** 2) / prev_count
        self.ucb_var[state0][action] = var  # 更新

        # 状態がなければ新規追加
        if state not in self.ucb_count:
            # 0を回避するために最低を1で初期化
            self.ucb_count[state] = [1 for _ in range(agent.nb_actions)]
            self.ucb_var[state] = [0 for _ in range(agent.nb_actions)]

        # 合計を出す(数式ではN)
        total = sum(self.ucb_count[state])
        total = math.log(total)  # 重そうな計算なので for の外で

        # 各アクションのUCB値を計算
        ucbs = []
        for i in range(agent.nb_actions):
            count = self.ucb_count[state][i]

            # 平均
            ave = agent.q_table[state][i] / count
            # 分散
            var = self.ucb_var[state][i]

            # 数式を計算
            v = var + math.sqrt(2 * total / count)
            if 1/4 < v:
                v = 1/4
            tmp = ave + math.sqrt( (total / count) * v )
            ucbs.append(tmp)

        # ucbが最大値となるアクションを選択
        action = np.argmax(ucbs)

        # 選択したアクションの選択回数を増やす
        self.ucb_count[state][action] += 1

        return action

参考
・Algorithms for the multi-armed bandit problem(論文)

UCB-V

UCB1-Tuned で分散を意識する事が性能の向上につながる事が分かったので、更に分散を意識したものが UCB-V となります。

UCBV_a= \bar{x}_a + \sqrt{\frac{2 \hat{x}_a E(N) }{n_a}} + c \frac{3b E(N)}{n_a}

E(t) = \zeta\log(t)

記号	意味
$\bar{x}_a$	取得した報酬の平均、アームの報酬平均
$\hat{x}_a$	取得した報酬の分散、アームの報酬分散
$E(t)$	探索関数
$\zeta$	定数 (>0)
c	定数 (>0)
b	定数(?)

UCBv.py

class UCBv():
    def __init__(self):
        # 選択回数を数える為の変数を追加
        self.ucb_count = {}
        self.ucb_var = {}  # 分散用

    def save_weights(self, d):
        d["ucb_count"] = self.ucb_count
        d["ucb_var"] = self.ucb_var
    
    def load_weights(self, d):
        self.ucb_count = d["ucb_count"]
        self.ucb_var = d["ucb_var"]

    def select_action(self, agent, state):

        # 前回の報酬を計算
        state0 = agent.prev_observation  # 前回の状態
        action = agent.prev_action  # 前回のaction
        reward = agent.prev_reward  # 前回の報酬
        if state0 not in self.ucb_count:
            self.ucb_count[state0] = [1 for _ in range(agent.nb_actions)]
            self.ucb_var[state0] = [0 for _ in range(agent.nb_actions)]

        # 分散を更新
        prev_count = self.ucb_count[state0][action]
        prev_ave = agent.q_table[state0][action] / prev_count
        var = self.ucb_var[state0][action]
        var += ((reward - prev_ave) ** 2) / prev_count
        self.ucb_var[state0][action] = var  # 更新

        # 状態がなければ新規追加
        if state not in self.ucb_count:
            self.ucb_count[state] = [1 for _ in range(agent.nb_actions)]
            self.ucb_var[state] = [0 for _ in range(agent.nb_actions)]

        # 合計を出す(数式ではN)
        total = sum(self.ucb_count[state])

        # 各アクションのUCB値を計算
        zeta = 1.2
        c = 1
        b = 1
        e = zeta*math.log(total)
        ucbs = []
        for i in range(agent.nb_actions):
            count = self.ucb_count[state][i]

            # 平均
            ave = agent.q_table[state][i] / count
            # 分散
            var = self.ucb_var[state][i]

            tmp = ave + math.sqrt( (2*var*e)/count ) + c* (3*b*e)/count
            ucbs.append(tmp)

        # ucbが最大値となるアクションを選択
        action = np.argmax(ucbs)

        # 選択したアクションの選択回数を増やす
        self.ucb_count[state][action] += 1

        return action

参考
・Exploration-exploitation trade-off using variance estimates in multi-armed bandits(論文)
・https://github.com/jkomiyama/banditlib/blob/master/policy/policy_ucbv.hpp

KL-UCB

KL divergence とは２つの確率分布の差異を計る尺度です。
KL divergence を元に UCB の Regret 下限を厳密化したアルゴリズムが KL-UCB となります。
Regret とは多腕バンデットアルゴリズムの評価値で、Regretが下限という事は理論上もっとも性能がいいという事になります。

KL divergence を出す関数は以下です。

$$ d(p,q)= p\log{\frac{p}{q}} + (1-p) \log{\frac{1-p}{1-q}} $$

これを元に以下で値を出します。

$$ \max \Bigr( n_a d(\frac{x_a}{n_a},q) \leq \log(N) + c\log(\log(N)\Bigl) $$

この式が最大となる$q$を見つけ、その値がKL-UCB値にとなります。

記号	強化学習	(多腕バンディットだと…?)
$x_a$	$s_t$における取得した報酬の合計	アームの取得報酬の合計
$n_a$	$s_t$における$a_t$の選択回数	アームの選択回数
c	定数、定理上は 3 だけど 0 の方がパフォーマンスが出るとの事	〃

KL_UCB.py

class KL_UCB():
    def __init__(self, C=0, delta=1e-6, eps=1e-10):
        # 選択回数を数える為の変数を追加
        self.ucb_count = {}
        self.C = C
        self.delta = delta  # 探索幅
        self.eps = eps      # 探索の許容誤差

    def save_weights(self, d):
        d["ucb_count"] = self.ucb_count
    
    def load_weights(self, d):
        self.ucb_count = d["ucb_count"]

    def select_action(self, agent, state):

        # 状態がなければ新規追加
        if state not in self.ucb_count:
            # 0を回避するために最低を1で初期化
            self.ucb_count[state] = [1 for _ in range(agent.nb_actions)]

        # 合計を出す(数式ではN)
        total = sum(self.ucb_count[state])

        # 右辺をだしておく
        logndn = math.log(total) + self.C * math.log(math.log(total))
        
        # 各アクションのUCB値を計算
        ucbs = []
        for i in range(agent.nb_actions):
            count = self.ucb_count[state][i]
            p = agent.q_table[state][i] / count

            # 例外処理：p は 0～1
            if p > 1:
                ucbs.append(1)
                continue
            
            # 最大値を探索する
            q = p + self.delta
            converged = False  # debug
            for _ in range(20):
                # kl-divergence
                try:
                    kl = p * math.log(p/q) + (1-p) * math.log((1-p)/(1-q))
                except ValueError:
                    break

                if count * kl + self.eps >= logndn:
                    converged = True
                    break
                
                q += self.delta

            # debug
            #assert not converged, "WARNING:Newton iteration in KL-UCB algorithm did not converge!! p={} logndn={}".format(p, logndn)
            ucbs.append(q)

        # ucbが最大値となるアクションを選択
        action = np.argmax(ucbs)

        # 選択したアクションの選択回数を増やす
        self.ucb_count[state][action] += 1

        return action

参考
・カルバック・ライブラー情報量
・The KL-UCB Algorithm for Bounded Stochastic Bandits
and Beyond(論文)
・確率的バンディット問題
・https://github.com/jkomiyama/banditlib/blob/master/policy/policy_klucb.hpp

Thompson Sampling

Thompson Samplingは、ベイズ推定を元にしたアルゴリズムです。
Thompson Sampling も Regret が下限である(理論上はもっとも性能がいい)アルゴリズムとなります。

ベイズ推定

ベイズ推定は、ベイズ確率を元にした推定方法です。

一般的に確率というと、確率は変わらずに事象を重ねるとその確率に近づく(無限回繰り返せば確率通りになる)ものを指すと思います。
しかし、ベイズ確率は実際に観測された事象を元に確率を求める考え方です。
(厳密ではない言い回しのなので悪しからず…)

例えばここにコインがあり、2回投げたら2回とも表が出たとします。
そのあと表が出る確率は一般的な確率ですと50%ですが、ベイズ確率では66.6% ($\beta(3,1)$)になります。
表の確率が増えているのは、表が出ているので次も表が出るだろう(コインが歪んでいたりして表が出やすいのだろう)という情報を加味しているためとなります。

ベイズ推定では確率がどういった分布をとるかを仮定する必要があり、その仮定が重要となります。
コインの例ですとベルヌーイ分布(表と裏のどちらかしかとらない分布)を仮定しています。

この分布はどんな分布でもいいのですが、特定の分布にすると計算がすごく楽になります。
計算が楽になる分布を「共役事前分布」といい、いくつかあります。

今回は、その中からベータ分布と正規分布の場合を実装します。

参考
・史上最強図解　これならわかる！ベイズ統計学 単行本 – 2012/2/21
・5分でスッキリ理解するベイズ推定
・【ベイズ統計】共役事前分布とは？わかりやすく解説

β分布版

ベイズ推定にて一番簡単な推定です。
取得するデータがベルヌーイ分布に従う場合、事前分布をベータ分布と置くと事後分布もベータ分布となります。(共役事前分布)
ベルヌーイ分布はコインの表か裏かや、成功か失敗か等の2値しかとらない分布となります。

参考
・ベータ分布をPythonで書く
・ベータ分布の事後分布の平均と分散【ベイズ】
・https://github.com/jkomiyama/banditlib/blob/master/policy/policy_ts.hpp

更新式

ベータ分布はベイズ更新の式すら省略できます。
ベータ分布のαが成功回数、βが失敗回数となります。

Thompson Sampling(ベータ分布)の実装

報酬をベルヌーイ分布で表現しないといけないので、報酬が正なら成功、負なら失敗として回数を数えています。

ThompsonSamplingBeta.py

class ThompsonSamplingBeta():
    def __init__(self):
        # 報酬結果を保存する変数を作成
        self.reward_count = {}

    def save_weights(self, d):
        d["reward_count"] = self.reward_count
    
    def load_weights(self, d):
        self.reward_count = d["reward_count"]

    def select_action(self, agent, state):

        # 前回の報酬を加算
        state0 = agent.prev_observation  # 前回の状態
        action = agent.prev_action  # 前回のaction
        if state0 not in self.reward_count:
            # [報酬ありの回数、報酬なしの回数] で初期化
            self.reward_count[state0] = [ [1,1] for _ in range(agent.nb_actions)]
        if agent.prev_reward > 0:
            self.reward_count[state0][action][0] += 1
        else:
            self.reward_count[state0][action][1] += 1

        # 状態がなければ新規追加
        if state not in self.reward_count:
            # [報酬ありの回数、報酬なしの回数] で初期化
            self.reward_count[state] = [ [1,1] for _ in range(agent.nb_actions)]
        
        # アクションを計算
        vals = []
        for i in range(agent.nb_actions):
            data = self.reward_count[state][i]

            # ベータ分布に従って乱数を生成
            v = np.random.beta(data[0], data[1])
            vals.append(v)

        # 乱数が最大値となるアクションを選択
        action = np.argmax(vals)
        return action

正規分布(ガウス分布)版

正規分布の別名がガウス分布です。
あまりベイズ更新にのみ特化したまとめ(計算過程がどうでもいい人)がなかったので更新式を別途載せています。

参考
・正規分布の事後分布の平均と分散【ベイズ】
・確率分布のベイズ推定(正規分布に従う場合)

更新式

更新式はこの本(史上最強図解　これならわかる！ベイズ統計学 単行本 – 2012/2/21))より抜粋しています。

事前分布 $N(\mu_0,\sigma_0^2 )$ に対し、正規分布$N(\mu,\sigma^2 )$に従うデータn個を与えた後の事後分布 $N(\mu_1,\sigma_1^2 )$ は以下となります。
($\mu$ は平均、$\sigma^2$は分散、$\bar x$はn個のデータの平均値)

$$ \mu_1 = \frac{ \frac{n\bar x}{\sigma^2} + \frac{ \mu_0 }{\sigma_0^2} }{ \frac{n}{\sigma^2} + \frac{1}{\sigma_0^2} }$$

$$ \sigma_1^2 = \frac{1}{ \frac{n}{\sigma^2} + \frac{1}{\sigma_0^2} }$$

強化学習ではデータは $x$ の1個で更新する事になるので、以下に変形しておきます。

$$ \mu_1 = \frac{ \frac{x}{\sigma^2} + \frac{ \mu_0 }{\sigma_0^2} }{ \frac{1}{\sigma^2} + \frac{1}{\sigma_0^2} }$$

$$ \sigma_1^2 = \frac{1}{ \frac{1}{\sigma^2} + \frac{1}{\sigma_0^2} }$$

Thompson Sampling(正規分布)の実装

正規分布を用いた場合、報酬が正規分布に従うと仮定されます。
また、報酬の分散($\sigma^2$)が事前にわかっている事も前提となります。

ThompsonSamplingGaussian.py

class ThompsonSamplingGaussian():
    def __init__(self, dispersion=1):
        # 既知の分散
        self.dispersion = dispersion

        # 現在の分散と平均
        self.recent_sigma = {}
        self.recent_mu = {}

    def save_weights(self, d):
        d["recent_sigma"] = self.recent_sigma
        d["recent_mu"] = self.recent_mu
    
    def load_weights(self, d):
        self.recent_sigma = d["recent_sigma"]
        self.recent_mu = d["recent_mu"]

    def select_action(self, agent, state):
        # 前回の報酬を加算
        reward = agent.prev_reward
        state0 = agent.prev_observation  # 前回の状態
        action = agent.prev_action       # 前回のaction

        # なければ追加
        if state0 not in self.recent_sigma:
            self.recent_sigma[state0] = [1 for _ in range(agent.nb_actions)]
            self.recent_mu[state0] = [0 for _ in range(agent.nb_actions)]
        sigma = self.recent_sigma[state0][action]
        mu = self.recent_mu[state0][action]

        # 更新(平均)
        tmp1 = reward/self.dispersion + mu/sigma
        tmp2 = 1/self.dispersion + 1/sigma
        self.recent_mu[action] = tmp1/tmp2

        # 更新(分散)
        self.recent_sigma[action] = 1/( (1/self.dispersion) + (1/sigma) )

        # なければ追加
        if state not in self.recent_sigma:
            self.recent_sigma[state] = [1 for _ in range(agent.nb_actions)]
            self.recent_mu[state] = [0 for _ in range(agent.nb_actions)]
        sigma = self.recent_sigma[state]
        mu = self.recent_mu[state]

        # アクションを計算
        vals = []
        for i in range(agent.nb_actions):
            # 正規分布に従い乱数を生成
            v = np.random.normal(mu[i], sigma[i])
            vals.append(v)

        # 乱数が最大値となるアクションを選択
        action = np.argmax(vals)
        return action

DQNへ応用

論文を探せばあるかもしれませんが特に探していません。
本記事オリジナルの実装です。

UCB や TS アルゴリズムを DQN に応用する場合、ある状態におけるアクション回数だったり報酬の分散だったりを数える必要があります。
DQN である状態とはNNモデルで近似化されてた状態なので、状態を数えるだけのNNモデルを別途作成して数えるというのがアイデアです。

実装例

ある状態(state)を元に数えている変数をNNモデル化します。
UCB1_Tunedを例に挙げると以下です。

class UCB1_Tuned():
    def __init__(self):
        pass

    def compile(self, agent, optimizer, metrics):
        # 選択回数を数える為の変数を追加
        self.ucb_count = agent.build_network()
        self.ucb_count.compile(loss=clipped_error_loss, optimizer=optimizer, metrics=metrics)
        self.ucb_var = agent.build_network()
        self.ucb_var.compile(loss=clipped_error_loss, optimizer=optimizer, metrics=metrics)

    def save_weights(self, d):
        d["ucb_count"] = self.ucb_count.get_weights()
        d["ucb_var"] = self.ucb_var.get_weights()
    
    def load_weights(self, d):
        self.ucb_count.set_weights(d["ucb_count"])
        self.ucb_var.set_weights(d["ucb_var"])

    def select_action(self, agent, state):
        # 前回の報酬を計算
        state0 = agent.get_prev_observation()  # 前回の状態
        action = agent.get_prev_action()  # 前回のaction
        reward = agent.get_prev_reward()  # 前回の報酬
        
        counts = self.ucb_count.predict(state0, 1)[0]
        qvals = agent.model.predict(state0, 1)[0]
        ucb_vars = self.ucb_var.predict(state0, 1)[0]

        # 分散を更新
        prev_count = counts[action]
        prev_ave = qvals[action] / prev_count
        var = ucb_vars[action]
        var += ((reward - prev_ave) ** 2) / prev_count
        ucb_vars[action] = var
        # 更新
        self.ucb_var.train_on_batch(state0, np.asarray(ucb_vars))

        counts = self.ucb_count.predict(state, 1)[0]
        qvals = agent.model.predict(state, 1)[0]
        ucb_vars = self.ucb_var.predict(state, 1)[0]

        # 合計を出す(数式ではN)
        total = sum(counts)
        total = math.log(total)  # 重そうな計算なので for の外で

        # 各アクションのUCB値を計算
        ucbs = []
        for i in range(agent.nb_actions):
            count = counts[i]

            # 平均
            ave = qvals[i] / count
            # 分散
            var = ucb_vars[i]

            # 数式を計算
            v = var + math.sqrt(2 * total / count)
            if 1/4 < v:
                v = 1/4
            tmp = ave + math.sqrt( (total / count) * v )
            ucbs.append(tmp)

        # ucbが最大値となるアクションを選択
        action = np.argmax(ucbs)

        # 選択したアクションの選択回数を増やす
        counts[action] += 1
        self.ucb_count.train_on_batch(state, np.asarray(counts))

        return action

ひとつ前の情報は分かりやすく関数化しています。

def get_prev_observation(self):
    return np.asarray([self.recent_observations[-self.input_sequence:]])
def get_prev_action(self):
    return self.recent_action[-1]
def get_prev_reward(self):
    return self.recent_reward[-1]

他の policy の書き換えは全体コードを見てください。

比較

ゲームは今まで通り Pendulum-v0 で比較します。
Processor も以前作成したものを使います。

比較内容として、学習経過と学習までにかかった時間、学習後に100回テストして得られた報酬の平均を出力しています。

また、事前に言っておきますが方策による学習はゲームの相性にものすごく左右されます。
ゲームが違うと全然違う結果になったりするのでこの結果は参考程度に見てください。

Q学習の比較結果

学習回数は1_000_000回です。

まずは、学習過程で取得している報酬です。
これは量が多くて見づらかったので1000個の移動平均でグラフ化しています。

Greedy、焼きなましGreedy、Softmaxが少し立ち上がりが遅い感じですが、あまり大差はありませんね。

次に学習終了までかかった時間です。(単位は秒)

KL-UCB が遅いですが他はあまり変わりません。
KL-UCBが遅いのは値を出すためにforによる探索(最急降下法?)を実施しているからです。

最後に、100回テストした結果の平均報酬です。
(-1600ぐらいが最悪、0が最高)

TS(正規分布)が良い値になっていますね。
UCB-V も学習効率が高そうです。

ε-Greedy系列は乱数の影響がかなりでます。
最初に良い乱数を引けないとそのまま局所解に入ってしまい学習できなくなったりします。
(実際、GoogleColaboratoryの結果ではε-Greedyが良い結果を残しています)

DQNの比較結果

学習回数は30_000回です。

まずが学習過程の取得している報酬です。
こちらも同様に50個の移動平均をグラフ化しています。

Greedy系列の学習が高いですね。

次に学習時間です。

Greedy系列とSoftmaxが早いのはNNモデルがないからです。(他はNNモデルの学習があります)
だいたい1.5～3倍ぐらい時間が違いますね。
また、ここはQ学習と同じ単位なので比較ができますが、Q学習に比べ学習回数こそ少ないですが、学習時間はこちらの方がかかっています。

最後に、100回テストした結果の平均報酬です。

ここはなぜか学習過程では報酬が低いにもかかわらず学習できているアルゴリズムがいくつかあります。
(テストはQ値最高のみを見ているのでその差分でもあるのでしょうか)
TS(ベータ分布)は報酬を0と1しか見れないので学習できていないのはまあそうでしょう。
他はUCB-V以外は学習できていますね。(UCB-Vが出来ていないのは…乱数?)

あとがき

実装してみてわかりましたが速さを重視して ε-Greedy で実装はありだと思いました。
ただ、乱数の影響とゲームの種類の影響が大きいので使い分けするのもありだと思います。

R2D2やApe-xでは各Actorが複数いるので、Actor毎に違う探索ポリシーを使うのはありかもしれませんね。
今後の強化学習のアルゴリズムに期待です。

さて、今まで実装ばっかりでしたが、次回は各アルゴリズムについて比較してみたいと思います。