ベイズ最適化で最高のコークハイを作る

はじめに

コークハイとか酎ハイをお店で飲むと、割り方とかレモンが効いていたりとかでお店によって結構違いが出ますよね
自分好みの最高のコークハイの作り方を知ることは全人類の夢だと思います。
本記事は一足先にそんな夢に挑戦したという記事です。

手法としてはベイズ最適化を使用します。

実データで実験計画と絡めながらベイズ最適化を実際に行う記事はあまり見かけなかったので今回は、
最適化パラメータ
　1. コーラとウイスキーの比
　2. レモン汁の量
目的変数
　コークハイの美味しさ

という2次元入力、1次元出力で実際に実験とチューニングを並行しながら行ってみたいと思います。

目次

ベイズ最適化とは

ブラックボックス最適化で使われる最適化手法のひとつです。
ハイパーパラメータチューニングに使われるOptunaの中身もその一種です。*1

入力xとそれに対応する出力yは分かっているけど、x, yの関係を説明する関数fが分からないときに
過去のデータからガウス過程回帰等*2でブラックボックス関数fの回帰モデルを作成し
その回帰モデルをもとに獲得関数の最大化を行うことで、最適化を効率的に行うというものです。

パラメータサーチで一番わかりやすい方法は
想定される全てのxについて細かくグリッドサーチ的に調べあげ、ノイズが多ければ反復的にデータを取得すること、ですが
時間や費用など多大なコストが必要な実験では現実的ではありません。

そういった実験系ではベイズ最適化による効率的探索が重要となり、応用ドメインはWEBでのABテストから地質調査、材料力学、生物学、クッキーのレシピ*3...などかなり広いです。

ガウス過程回帰(+ベイズ最適化)について数式込みで詳しく知りたい方は
www.kspub.co.jp

ベイズ最適化については
明治大学講演資料「機械学習と自動ハイパーパラメタ最適化」 佐野正太郎
ベイズ的最適化(Bayesian Optimization)の入門とその応用
シンプルなベイズ最適化について | AI tech studio

などありがたいことに資料が公開されています。
あとはみんなでPRMLを頑張りましょう。

実験系の説明

実験条件

実験で考慮しないこと(パラメータ)

コークハイの割り方以外にも味の評価に影響しそうなパラメータはいくつか思いつきますが、実験系の単純化のために今回は考慮しません

• 飲酒のシチュエーションに関するパラメータ

　周囲の気温とか水温は味に影響しそうですが考慮しません
　　例：夏の暑い日に冷たいものを飲むと美味しく感じる

• 被験者の体調に関するパラメータ

　疲労具合であったり、満腹度合いや、喉が乾いているかどうかだったりは影響しそうですが考慮しません
　　例：疲れているときにレモンをかじると甘く感じるとかの話
　
　(酔っ払い具合も影響しそうですが、これについては実験にゆっくり時間をかけて、少量のコークハイを摂取するという形で対処を行っています。)

実験材料

• ウイスキー(シーバスリーガル12年 700ml*4 )
• コカコーラ(500ml)
• レモン汁(ポッカサッポロ　ポッカレモン70ml)
• おつまみ(ハッピーターンわさびマヨ *5 )
• 緑茶

それぞれ常温のものを使用。

◯コークハイの作り方
コーラの量を0~50ml加えた後に、コーラ+ウイスキーの合計が50mlになるようにウイスキーを加える。
その後レモン汁を0~15ml加えてかき混ぜる。

実験方法

1. コークハイを50ml作成する
2. 味のリセット(おつまみを食べて、茶を飲む)
3. コークハイを飲む
4. スコアをつけ、ベイズ最適化を行う
5. 器具の洗浄
6. 10分前後時間を開けた後に1~5を繰り返す

スコアの付け方

スコアの下限：0
スコアの上限：100

◯スコアの目安
　・50mlなら頑張って飲める：40~
　・中ジョッキ(500ml)飲める：60~
　・(コークハイとして)美味しい：80~

基本的に美味しいか美味しくないかで点数を付けています。
が、それだったらありのままのコーラが100点取って終わりなので、コークハイとしてクオリティで点数を付けます。(1度限りの相対評価なので怪しい部分はあります)

実験をやりました（本題）

実装にはGpyを使用しました。
ガウス過程回帰に使用したカーネルはシンプルなRBFカーネル

既知のデータ点からの距離をみて、入力変数の値が近いデータ同士は出力も近い値が出るよねという表現

最適化対象の獲得関数はUCB(Upper Confidence Bound)

この式は"活用"と"探索"を表現していて
活用：平均値が高いところをより調べたい
探索：データ点が少ないところを調べたい(今回はデータ数nが大きくなれば係数が小さくなって影響が小さくなる)
という感じです

入力変数についてはそれぞれ0~1の範囲にMinMax正規化を行いました。
コーラの量0~50　→　コーラの割合0~1
レモン汁の量0~15　→　レモン汁の量0~1

また計量カップの都合上
コーラの割合は0.1刻み(元の量で換算すると5ml刻み)
レモン汁の量は0.2刻み(元の量で換算すると3ml刻み)としました。

初期点として
(コーラ割合, レモン汁量) = (1, 1), (0, 1)
すなわち
コーラ50ml + レモン15ml
コーラ0ml (ウイスキー50ml) + レモン0ml
の2点を既知として実験を行いました。

実装コード

import GPy
import matplotlib.colors as colors
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.cm import ScalarMappable

""" 実験データの逐次記入 """
df.loc[0] = [50, 15, 68] #初期点
df.loc[1] = [0, 0, 44]   #初期点
df.loc[2] = [25, 15, 58] #実験1回目のデータ
df.loc[3] = [50, 0, 85]  #2回目
df.loc[4] = [45, 0, 94]  #3回目...
df.loc[5] = [25, 6, 50] 
df.loc[6] = [50, 6, 90]
df.loc[7] = [0, 15, 46]

"""データサイズn"""
n = len(df)

from sklearn.preprocessing import MinMaxScaler

""" 入力変数をMinMaxスケーリング """
scaler = MinMaxScaler()
scaler.fit(np.array([0,50]).reshape(-1,1)) #0~50を0~1に
df["coke"] = scaler.transform(df["coke"].values.reshape(-1,1))

scaler = MinMaxScaler()
scaler.fit(np.array([0, 15]).reshape(-1,1)) #0~15を0~1に
df["lemon"] = scaler.transform(df["lemon"].values.reshape(-1,1))

"""入力の整形"""
x_train = np.stack([df["coke"], df["lemon"]], axis=1)
y_train = df["score"].values

"""既知のデータをもとにガウス過程回帰"""
kern = GPy.kern.RBF(2)#, ARD=True)
gpy_model = GPy.models.GPRegression(X=x_train.reshape(-1, 2),
                                    Y=y_train.reshape(-1, 1),
                                    kernel=kern,
                                    normalizer=True)

"""モデルの可視化"""
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)

gpy_model.optimize() # カーネル最適化
gpy_model.plot_mean(ax=ax, cmap="jet")  # カーネル最適化後の予測


"""データ点の描画(defaultでは見辛い)"""
[ax.plot(xi, yi, marker=".", color="k", markersize=10) for xi, yi in zip(df["coke"].values, df["lemon"].values)]
ax.set_ylabel("normalized_lemon", fontsize=18)
ax.set_xlabel("normalized_coke", fontsize=20)
ax.tick_params(labelsize=20)


"""color bar追加""" 
axpos = ax.get_position()
cbar_ax = fig.add_axes([1, 0.15, 0.02, 0.8])
norm = colors.Normalize(vmin=y_train.min(), vmax=y_train.max())
mappable = ScalarMappable(cmap='jet',norm=norm)
cbar_ax.tick_params(labelsize=10) 
fig.colorbar(mappable, cax=cbar_ax)

fig.tight_layout()
plt.savefig(f"2dim_gaussian_n={n}.png")

"""
作った回帰モデルをもとに新たな入力xに対する出力yをみて
獲得関数acqを最大化する点(=次に実験を行うべき条件)を探す。
"""
x = np.linspace(0, 1.0, 11)
acq_list = []
for i in range(6):
    x_pred = np.array([x, np.full(11, i*0.2)]).T
    y_mean, y_var = gpy_model.predict(x_pred)
    acq = (y_mean + ((np.log(n) / n) ** 0.5 * y_var)).max()
    acq_list.append(acq)
    
next_lem = acq_list.index(max(acq_list)) * 0.2 #次のレモンの分量
print("max lemon:", next_lem)


""" 入力を1次元固定して予測(レモン)"""
x_pred = np.array([x, np.full(11, next_lem)]).T
y_mean, y_var = gpy_model.predict(x_pred)

gpy_model.plot(fixed_inputs=[(1, next_lem)], plot_data=False, plot_limits=[-.01, 1.01])
plt.xlabel("normalized coke", fontsize=14)
plt.ylabel("taste score [ g ]", fontsize=14)

"""獲得関数acq"""
acq = (y_mean + ((np.log(n) / n) ** 0.5 * y_var)) / 5  

"""獲得関数の可視化"""
plt.plot(np.linspace(-.01, 1.01, 11), acq, color="g")
plt.plot(acq.argmax() * 0.1,
         acq.max(),
         marker=".", color="r", markersize=14)

plt.legend(["Mean", "Acquisition", "Acq Max", "Confidence"])
plt.savefig(f"1dim_gaussian_n={n}_lemon={next_lem}.png")

実験開始

やっとコークハイが飲めます...長い道のりでした。
初期点の条件で前述の通りコークハイを作ってスコアを付けます。飲みながら書いた当時のコメント付きです。

• コークハイ1

　ウイスキー：50[ml]
　コーラ：0[ml]
　レモン：0[ml]
　スコア：44
　厳しい…普通にロックですらない、ありのままのウイスキー
　香りがめっちゃ強い。飲めるけどきつい。本音を言うと5mlで良い。

この時点でコークハイというかコーラが好きなだけでウイスキーを好んで飲めるほど味覚が大人でなかったことを思い出しますが、ベイズに全て委ねて実験を続けます。

• コークハイ2

　ウイスキー：0[ml]
　コーラ：50[ml]
　レモン：15[ml]
　スコア：68
　今度はレモンが強すぎる。レモン汁のコーラ風味

上記の2点の初期点をもとにガウス過程回帰を行うと次のようになりました。

next → lemon=1.0, coke=0.5

見方としては上の図の等高線のようなものが平均値を表し、下図はそれを横からみたようなイメージです。
下の図を見るとまだデータが少なく、データ点同士の距離が大きいのでどの範囲でも分散が大きいですね
獲得関数が最大となる条件を見ると、次はlemon=1.0, coke=0.5の条件を調べろとのことです。

実験1回目

• コークハイ3

　ウイスキー：25[ml]
　コーラ：25[ml]
　レモン：15[ml]
　スコア：58
　感覚としてはウイスキーのレモン汁わり
　まだまだウイスキー成分が多い気がするけどだいぶ飲める。
　この辺りで酔ってくる。

next → lemon=0.0, coke=1.0

lemon=0.0の点は既知なのでその周辺は分散が小さくなっています。
上の図を見ると、右上にいくほどスコアが良いことを表しているので
現状ではコーラとレモンが多いほど美味いんでしょ？？？と言われている感じ

実験2回目

• コークハイ4

　ウイスキー：0[ml]
　コーラ：50[ml]
　レモン：0[ml]
　スコア：85
　コーラは美味しい

next → lemon=0.0, coke=0.9

点差がばらけすぎて描画がイマイチになってしまいました

実験3回目

• コークハイ5

　ウイスキー：5[ml]
　コーラ：45[ml]
　レモン：0[ml]
　スコア：94
　とても飲みやすいけどアルコールも感じてちょうどいい.
　コークハイとして美味しいので高得点。

next → lemon=0.4, coke=0.5

右下辺りがスコア高いんじゃないかという示唆が出てきました。良い感じです。

実験4回目

• コークハイ6

　ウイスキー：25[ml]
　コーラ：25[ml]
　レモン：6[ml]
　スコア：50
　コークハイ3に比べてレモン汁が減ってウイスキー感(というかアルコール感)が強い
　もう少しコーラが欲しい

next → lemon=0.4, coke=1.0

　次の条件はコーラ=1.0, lemon=0.4となり、またコークハイ以外のものを飲まされるようです。そろそろ不安です。

実験5回目

• コークハイ7

　ウイスキー：0[ml]
　コーラ：50[ml]
　レモン：6[ml]
　スコア：90
　コーラにレモンはとても合う
　あくまでコークハイを作りたいので少しスコアは控えめに

next → lemon=1.0, coke=0.0

コーラが多いほど飲みやすいという感覚に近い結果が現れてきました。
とはいえコークハイを飲まずにコーラばかり飲んだ弊害か、スコアの高い部分が右下に振り切れてしまいました。

実験6回目

• コークハイ8

　ウイスキー：50[ml]
　コーラ：0[ml]
　レモン：15[ml]
　スコア：46
　助けてコーラ…コーラ無しウイスキーは厳しい...

next → lemon=0.0, coke=1.0

次の実験条件に既に済んだ条件が提示されてしまいました。
これは雑に0.1刻みに丸め込んだことや、探索する範囲があまり多くないことが原因として考えられそうです。

上の図を見るとcoke < 0.5の範囲についてはスコアが低いので熱心に探さなくても良いことが感覚ではなく、モデリングの結果から判断出来るは嬉しいポイントです。
全部で50通り考えられる入力変数の組み合わせのうち、8回の実験でこの結果が得られるのは予備実験としては中々良いと思います。

とはいえ点数の付け方的にも右下が一番良いというのは望ましい結果でなく、coke=0.7辺りの直感的に美味しく作れそうな範囲をもう少し調べて欲しいという思いは残ります(この辺りは後述の初期点や、変数の範囲設定の問題)

ARDありver.

ここでガウス過程回帰のカーネルにARD(関連度自動決定カーネル)を使用してみます。
これまでは2つの変数コーラの割合、レモン汁の量を同列に扱っていましたが、ARDを使用するとそれぞれに重み付けをして最適化を行います。
先ほどまでのデータを使って、ARDを使用した結果の変化もみてみたいと思います。

(ドメイン知識(=僕自身のコークハイの好み)的にもレモン汁の量よりもウイスキーがどれくらい薄まるか( ≒ コーラの割合)の方が変数として重要なのでこの2つが同列であることは考えにくいです。)

実験1回目(ARD)

実験2回目(ARD)

実験3回目(ARD)

実験4回目(ARD)

実験5回目(ARD)

実験6回目(ARD)

実験7回目(ARD)

一番下の図を見ると、ウイスキーの割合が多いうちは変数cokeが支配的でレモンの量に関わらずスコアが低い一方で、coke > 0.7辺りからレモン汁の量もスコアに影響を及ぼすという結果が得られました。こっちの方が納得感ありますね。

ARDでの事後分布(lemon=0.0)

反省点

今回はコーラの割合とレモン汁の量を最適化し、ARDを使うことでそれっぽい回帰モデルを作りパラメータサーチをすることができました。ただしいくつかの改善点や反省点があり

・初期点の荒さ
今回は初期点にウイスキーonlyの点と、コーラ+レモン汁の点のみを使用しましたが、最初からウイスキーとコーラの比が1:3くらいのコークハイとして一般的な条件を初期点として設定しておくことでより効率的な探索が可能になったと思います。(結果的にコーラばっかり飲んでたのはあまりよろしくない)

・説明変数の設定
コーラの割合とレモン汁の量というシンプルな変数を使用しましたが、実際の実験を通してアルコール濃度( = ウイスキーの薄さ)がピンポイントで効くように感じたので変数を変えることで結果が変わるやもしれません。

また想定外の取りこぼしを防ぐために探索範囲を広めに取りましたが、今回の結果をもとに、今後は探索範囲を狭めても良さそうです。(これもデータを取得したから判断できることですが)

さいごに

今回の実験ではある程度雑な条件設定でも個人的に納得できる結果が得られ、実験計画と絡めたベイズ最適化の面白さを確認しました。
このまま続けていけば少ない試行回数で絞り込むことができそうです。

とはいえ使える分野、場面が多い分ドメイン知識を持った上で運用しないと最適化する必要のない変数をチューニングするような結果になるかもしれないので注意が必要であると思います。

次は体力があれば今回の結果をベースに角ハイを使ったマルチタスクベイズ最適化を行って更にコークハイを突き詰めたいと思います。

めちゃくちゃ長くなりましたが読んでくださりありがとうございました。

*1:ハイパーパラメータ自動最適化ツール「Optuna」公開 | Preferred Networks Research & Development

*2:最近はNeural Processesなんかもあるみたいです

*3:www.blog.google

*4:買ってから2年経っていたので14年モノかも

*5:全お菓子の中で2番目くらいに美味しい