焼肉最適化問題

　本稿はSFC-RG Advent Calendar 2016の4日目である．

はじめに

　あなたは研究の中間発表を終えて，今晩何を食べようか考えている．たしかに準備不足ではあったけれど，研究の前提をいまいち解さないファカルティの高飛車な質問にはうんざりしたし，今日くらいはパーッと気分転換したいものだ．そういうわけで，あなたは⊿館を飛び出して焼肉 ざんまい 湘南台店に行くことにした．

組合わせ最適化

　さて，着席し，メニューを開いたあなたはしばし考える．限られた予算，限られた時間，限られた胃袋の容量———いったい何を頼めば最も満足できるだろうか？
　そんなとき，組合わせ最適化が役に立つんです．騙されたと思って，メニューを必死に転記してみよう：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import numpy as np, pandas as pd

menu = ['カルビ', '和牛カルビ', '和牛中落ちカルビ', 'ハラミ', '厚切りロース', 'ネギ牛タン塩', '牛タン塩', 
        'イベリコ豚', 'カシラ', '豚トロ', 'ネギ豚タン塩', '豚タン塩', '厚切りベーコン', 'ウインナー', 'チョリソ', 
        'ホルモン', 'シロコロホルモン', 'レバー', '白レバー', 'ハツ', 'ミノ', 
        'お得!! 三種盛り', '本日の塩三種盛り', '本日の味噌三種盛り']
price = [720, 950, 850, 720, 690, 950, 850, 
        600, 550, 580, 680, 580, 500, 380, 400, 
        550, 600, 550, 450, 550, 650, 
        1280, 780, 780]

n = len(menu)
np.random.seed(0)
df = pd.DataFrame({
    '品目': menu,
    '値段': price,
    '満足度': np.random.randint(10, 20, n),
    '焼き時間': np.random.randint(5, 10, n),
    '量': np.random.randint(10, 20, n),
    })

print(df)
値段 品目 満足度 焼き時間 量
0 720 カルビ 15 9 10
1 950 和牛カルビ 10 8 10
2 850 和牛中落ちカルビ 13 5 14
3 720 ハラミ 13 8 15
4 690 厚切りロース 17 5 15
5 950 ネギ牛タン塩 19 7 16
6 850 牛タン塩 13 8 18
7 600 イベリコ豚 15 5 14
8 550 カシラ 12 6 11
9 580 豚トロ 14 8 14
10 680 ネギ豚タン塩 17 8 19
11 580 豚タン塩 16 8 18
12 500 厚切りベーコン 18 5 11
13 380 ウインナー 18 6 11
14 400 チョリソ 11 6 17
15 550 ホルモン 16 6 19
16 600 シロコロホルモン 17 5 19
17 550 レバー 17 7 13
18 450 白レバー 18 9 16
19 550 ハツ 11 8 17
20 650 ミノ 15 8 12
21 1280 お得!! 三種盛り 19 7 10
22 780 本日の塩三種盛り 18 9 13
23 780 本日の味噌三種盛り 19 7 15

　メニューがpandasのデータフレームになった．取り急ぎ（？）満足度・焼き時間・量は乱数で埋めている．
　あなたはこのメニューの中から，最も満足度の高い組合わせを見つけたい．それも，限られた予算，限られた時間，限られた胃袋の容量という条件を満たすような．ここではそのわがままを割当問題として解くことにする．

変数 $x_i \in \{0, 1\}$ i番目の料理を選ぶかどうか
目的関数 $\sum_i{満足度_i x_i}$ $\rightarrow$ 最大
制約条件 $\sum_i{予算_i x_i} \le 12000$ 予算制限
$\sum_i{焼き時間_i x_i} \le 120$ 時間制限
$150 \le \sum_i{量_i x_i} \le 200$ 分量制限

　こういった問題はpulpという整数計画法ライブラリを使って解くことができる：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
from pulp import *

m = LpProblem(sense = LpMaximize)	# 最大化問題
df['x'] = [LpVariable('x%d' % i, cat = LpBinary) for i in range(n)]	# i番目の品目を選択するかしないか

m += lpDot(df.満足度, df.x)	# 目的関数：満足度の最大化
m += lpDot(df.値段, df.x) <= 12000	# 制約条件：予算
m += lpDot(df.焼き時間, df.x) <= 120	# 制約条件：焼き時間
m += lpDot(df.量, df.x) >= 150	# 制約条件：量
m += lpDot(df.量, df.x) <= 200	# 制約条件：量

m.solve()
if m.status == 1:
    df['val'] = df.x.apply(lambda v: value(v))	# 結果
    print(df[df.val == 1].品目)
    print('満足度 {}'.format(sum(df[df.val == 1].満足度)))
    print('値段 {}'.format(sum(df[df.val == 1].値段)))

>>>
0           カルビ
4        厚切りロース
5        ネギ牛タン塩
7         イベリコ豚
8           カシラ
9           豚トロ
12      厚切りベーコン
13        ウインナー
16     シロコロホルモン
17          レバー
18         白レバー
20           ミノ
21    お得!! 三種盛り
22     本日の塩三種盛り
23    本日の味噌三種盛り
Name: 品目, dtype: object
満足度 251
値段 10060

　はい．順番はともかくとして，これらを食べれば満足できそうだ．
　ここまで考えたところで，あなたは今月の懐具合がよろしくないことを思い出す．なるべく出費を抑えて，それでいてある程度満足できるような品目の組合わせはあるだろうか？
　これも同様に割当問題として考えられる．

変数 $x_i \in \{0, 1\}$ i番目の料理を選ぶかどうか
目的関数 $\sum_i{予算_i x_i}$ $\rightarrow$ 最小
制約条件 $\sum_i{満足度_i x_i} \le 200$ 満足度制限
$\sum_i{焼き時間_i x_i} \le 120$ 時間制限
$150 \le \sum_i{量_i x_i} \le 200$ 分量制限 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
m = LpProblem(sense = LpMinimize)	# 最小化問題
a['v'] = [LpVariable('v%d' % i, cat = LpBinary) for i in range(n)]	# i番目の品目を選択するかしないか

m += lpDot(df.値段, df.x)	# 目的関数：予算の最小化
m += lpDot(df.満足度, df.x) >= 200	# 制約条件：満足度
m += lpDot(df.焼き時間, df.x) <= 120	# 制約条件：焼き時間
m += lpDot(df.量, df.x) >= 150	# 制約条件：量
m += lpDot(df.量, df.x) <= 200	# 制約条件：量
m.solve()
if m.status == 1:
    df['val'] = df.x.apply(lambda v: value(v))	# 結果
    print(df[df.val == 1].品目)
    print('満足度 {}'.format(sum(df[df.val == 1].満足度)))
    print('値段 {}'.format(sum(df[df.val == 1].値段)))

>>>
7         イベリコ豚
10       ネギ豚タン塩
11         豚タン塩
12      厚切りベーコン
13        ウインナー
14         チョリソ
15         ホルモン
16     シロコロホルモン
17          レバー
18         白レバー
22     本日の塩三種盛り
23    本日の味噌三種盛り
Name: 品目, dtype: object
満足度 200
値段 6850

　200の満足度でいいなら豚ばっか食ってろということらしい．

多腕バンディット問題

　いやいやちょっと待った．乱数でお茶を濁しているけど，あらゆる品目の満足度なんて知らないじゃないか．全品目を食べたことがあるならいざ知らず．それに，毎日同じ店で同じ食事をとるわけでもない．焼肉屋にしたって，湘南台には寅屋にえのもとにとあるわけだ．
　そういうわけで，あなたはいろいろな店に行き，いろいろな注文をして，なるべくどれを頼んでも満足度の高い食事のとれる店を見つけたいと考えた．しかしここで疑念が生まれる———いまの店でそこそこ満足できるなら，別の店を探さなくてもよいのではないか？ しかしいまの店に行き続ける限り，別の店の食事を食べることはできないぞ．しかしそこがハズレだとしたら．さてどうしよう———業界用語でこれを探索と利用のトレードオフ (exploration-exploitation tradeoff) という．

これまでの学習結果を利用 (exploitation) しようとすると，探索 (exploration) が減ってしまい，機会損失が増えてしまう．一方，探索を増やせば，学習した最良の行動とは異なる行動をとることが増えるため，得られる報酬が減ってしまう．学習した最良の行動との差が，探索のコストということになる．— 牧野，et al. 『これからの強化学習』

　このトレードオフを解消する試みを多腕バンディット問題という．多腕バンディット問題は，複数のスロットマシンのレバー（腕）を次々と引いていき，最終的に得られる報酬を最大化するというもので，強化学習の一種といえる．各スロットマシンの報酬はそれぞれ一定の確率分布に従っている．言い換えれば，いろいろな店のメニューにある品目を試していき，最終的に得られる満足度を最大化していく，ということになる．もちろん，品目によって得られる満足度に違いはあるが，なるべく何を食べても満足度の高い店を絞り込んでいきたい．
　そのためのアルゴリズムのうち，ここでは$\epsilon$-GreedyとUCB1を紹介したい．

$\epsilon$-Greedy

　デフォルトで現在最良な選択肢を選ぶが，一定の確率でいま最良と思っていないような選択肢を選びにいく手法．

  • $\epsilon - 1$の確率で最適だと思われる選択肢を利用する
  • $\epsilon / 2$の確率で最適だと思われる選択肢を探索する
  • $\epsilon / 2$の確率で最悪だと思われる選択肢を探索する

　$\epsilon$ Greedyはつまり行きつけの店に入り浸るタイプのアルゴリズムだ．

UCB1

　それもいいが，実はいまの行きつけよりもっといい店なのに，一度行って微妙だったからといって行かないままになっている店がないだろうか？ UCB1は，それぞれの店についてどれくらい知っているかを考慮に入れ，なるべく多くの店のことを知ろうとするアルゴリズムだ．具体的には，各店（腕）についてUCB (Upper Confidence Bound) という値を計算する．

　$\overline {x}_{j}+c\sqrt {\dfrac {2\log n} {x_{j}}}$

　ただし$\overline {x}_{j}$は$_{j}$番目の店の満足度（報酬）の期待値，$n$はそれまでに店を回った回数の合計，$n_{j}$は$_{j}$番目の店に行った回数，$c$は定数．UCB1は，この値が最大になる店に飛び込んでいく．あなたが好奇心旺盛なら，こちらのアルゴリズムを使って考えたほうがいいだろう．
　この手法のメリットとして，ベストでない店に行く回数を確率$1$で$O(\log n)$以内に抑えられる．長々とした証明は論文を参照していただくとして，これは店に行く回数が十分大きい場合の理論限界に到達している．デメリットとしては，探索のためによくない店にあたってしまうことが多いという点が挙げられる．

実験

　では，$\epsilon$-GreedyとUCB1が最良の店を選ぶ過程はどうなっているだろうか？ 『バンディットアルゴリズムによる最適化手法』のサンプルコードをPython3 + numpy向けに改変して実験．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import random
import math
import seaborn as sns

class BernoulliArm():
    def __init__(self, p):
        self.p = p

    def draw(self):
        if random.random() > self.p:
            return 0.0
        else:
            return 1.0

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

    def initialize(self, n_arms):
        self.counts = np.zeros(n_arms)
        self.values = np.zeros(n_arms)
        return

    def select_arm(self):
        if random.random() > self.epsilon:
            return np.argmax(self.values)
        else:
            return np.random.randint(0, len(self.values))

    def update(self, chosen_arm, reward):
        self.counts[chosen_arm] += 1
        n = self.counts[chosen_arm]

        value = self.values[chosen_arm]
        new_value = ((n-1) / float(n)) * value + (1 / float(n)) * reward
        self.values[chosen_arm] = new_value
        return

class UCB1():
  def __init__(self, counts, values):
    self.counts = counts
    self.values = values
    return

  def initialize(self, n_arms):
    self.counts = np.zeros(n_arms)
    self.values = np.zeros(n_arms)
    return

  def select_arm(self):
    n_arms = len(self.counts)
    for arm in range(n_arms):
      if self.counts[arm] == 0:
        return arm

    ucb_values = [0.0 for arm in range(n_arms)]
    total_counts = sum(self.counts)
    for arm in range(n_arms):
      bonus = math.sqrt((2 * math.log(total_counts)) / float(self.counts[arm]))
      ucb_values[arm] = self.values[arm] + bonus
    return np.argmax(ucb_values)

  def update(self, chosen_arm, reward):
    self.counts[chosen_arm] = self.counts[chosen_arm] + 1
    n = self.counts[chosen_arm]

    value = self.values[chosen_arm]
    new_value = ((n - 1) / float(n)) * value + (1 / float(n)) * reward
    self.values[chosen_arm] = new_value
    return

def test_algorithm(algo, arms, num_sims, horizon):
    chosen_arms = np.zeros(num_sims * horizon)
    rewards = np.zeros(num_sims * horizon)
    cumulative_rewards = np.zeros(num_sims * horizon)
    sim_nums = np.zeros(num_sims * horizon)
    times = np.zeros(num_sims * horizon)

    for sim in range(num_sims):
        sim += 1
        algo.initialize(len(arms))

        for t in range(horizon):
            t += 1
            index = (sim - 1) * horizon + t - 1
            sim_nums[index] = sim
            times[index] = t

            chosen_arm = algo.select_arm()
            chosen_arms[index] = chosen_arm

            reward = arms[chosen_arm].draw()
            rewards[index] = reward

            if t == 1:
                cumulative_rewards[index] = reward
            else:
                cumulative_rewards[index] = cumulative_rewards[index - 1] + reward

            algo.update(chosen_arm, reward)

    return [sim_nums, times, chosen_arms, rewards, cumulative_rewards]

means = np.array([0.1, 0.1, 0.1, 0.1, 0.9])
n_arms = len(means) # 腕は5本
random.shuffle(means)

arms = [BernoulliArm(x) for x in means]

for epsilon in [0.1, 0.2, 0.3, 0.4, 0.5]: # epsilonを変えたらどうなるか？
    algo = EpsilonGreedy(epsilon, [], [])
    algo.initialize(n_arms)
    results = test_algorithm(algo, arms, 5000, 500)

    df = pd.DataFrame({"times": results[1], "rewards": results[3]})
    grouped = df["rewards"].groupby(df["times"])

    plt.plot(grouped.mean(), label="epsilon="+str(epsilon))

algo = UCB1([], [])
algo.initialize(n_arms)
results = test_algorithm(algo, arms, 5000, 500)

df = pd.DataFrame({"times": results[1], "rewards": results[3]})
grouped = df["rewards"].groupby(df["times"])

plt.plot(grouped.mean(), label="UCB1")

plt.legend(loc="best")
plt.show()

　好奇心旺盛な人は序盤それなりに苦労することがわかる．

おわりに

　こうしたナイーブな実装で焼肉の頼み方を最適化したり，店の巡り方を最適化できるかどうかというとまあ実際微妙（たとえば焼肉の各品目から得られる満足度は，限界効用逓減の法則に従ってすり減っていく）だが，日々の意思決定をアルゴリズム的に考えてみる遊びはそれなりにおもしろい．ソーシャルゲームにどう課金するか，というのもこの俎上に載せられる．
　『Algorithms to Live By: The Computer Science of Human Decisions』という本はそういう，情報系の人間がよくやる与太話をまじめに考察したものだ———書類はどのキャッシュアルゴリズムに従って並べるべきかとか．先に挙げた探索と利用のトレードオフについても述べられている．YouTubeで著者らの講演を観てほしい．

　はい．

参考文献

追記 (2016.12.06)

　今回のような設定では多腕バンディット問題というより最適腕識別として考えたほうがよさそう．多腕バンディット問題は累積報酬の最大化が目的だけれど，最適腕識別はどの腕が最良か発見するのが目的．将来の報酬が最大の腕を見つける，ということ．『バンディット問題の理論とアルゴリズム』を読めばいいんだとか．