生物物理屋が数理最適化でブラックジャックに勝つ手を探してみた話

## はじめに
この記事は「[ただただアウトプットを癖付けるための Advent Calendar 2024](https://qiita.com/advent-calendar/2024/newbie2024)」に投稿した記事です。

[最初の記事](https://qiita.com/yotabouya/items/e63c53c5e18b3cf6a454#%E3%81%AF%E3%81%98%E3%82%81%E3%81%AB)にも書いた通り、私は生物物理の実験を専門にしている研究者です。
最近はデータ解析のため機械学習のコード開発も行っており、幸いにもその成果が[NeurIPSに採択されました](https://neurips.cc/virtual/2024/poster/93465)。

今回は、数理最適化に手を出してみた話を書きます。
数理最適化は、最適な解を求めるための数学的手法です。
学習のため、今回はブラックジャックでの行動を最適化するという課題に取り組んでみました。

## 関連記事

#### 前の記事 「[生物物理屋がAWSに手を出してみた話](https://qiita.com/yotabouya/private/b185754378cb338da070)」

#### 次の記事「[生物物理屋が数理最適化でブラックジャックに勝つ手を探してみた話 後編](https://qiita.com/yotabouya/private/40563070a5c9d8769ee9)」

## 数理最適化
数理最適化は、最適な解を求めるための数学的手法です。
最適な解を求めるためには、目的関数と制約条件を定義し、それらを満たす解を求めます。
実は今回の題材は当初、ポーカーにしようかと思っていました。
ただポーカーはブラックジャックと異なり、カードの数字だけでなくマークも重要な要素となるため、定式化の難易度が高くなります。
そのため、まずはブラックジャックでの行動を最適化することにしました。

## ブラックジャック
ブラックジャックは、ディーラーとプレイヤーがカードを引き、21に近いほうが勝利するゲームです。
まずそれぞれに2枚のカードが配られ、プレイヤーはカードを引くかスタンドするかを選択します。この選択はスタンドを選択するまで繰り返されます。
プレイヤーがスタンドを選択したら、ディーラーはプレイヤーの手よりも大きな合計値になるまでカードを引き続けます（プレイヤーが１人の場合）。
いずれも21を超えるとバーストとなり、即座に負けとなります。
なお、Aは1または11のどちらか都合のいい方で数えられます。
また、J,Q,Kは10として扱われます。

最適化の対象になるのは、プレイヤーの行動です。
プレイヤーはカードを引くかスタンドするかを選択しますが、この選択を最適化することで勝率を上げることができます。
使える情報としては、プレイヤーの手札とディーラーの1枚目のカードがあります。
というのも、ディーラーの1枚目のカードはオープンであり、プレイヤーが選択する際に参考にできるというルールがあるためです。
念のためですが、今回は表情などの情報は使いません。

## 定式化
まずは数理最適化のための定式化を行います。

### 手札の表し方
条件として、現状の手札の数字と、ディーラーの1枚目のカードの数字が与えられます。
これらの表し方ですが、プレイヤーの手札の枚数は変動しうるので、状態空間を固定するために、1から10までのカードの枚数を並べたベクトル

$$h:=(h_{i})_{i=1,\cdots,10} $$

で表すことにします。
ディーラーのカードについても同様に、1から10までのカードの枚数を並べたベクトル

$$d:=(d_i)_{i=1,\cdots,10}$$

で表します。

### 手札を引く行動
カードを引くという行動は、手札にカードを追加することで表現します。
追加されるカードは、残りのカードの中からランダムに選ばれると仮定します。
各数字のカードの合計枚数は決まっていて、

$$t:=(t_i)_{i=1,\cdots,10}$$

で表します。上記のルールから、

$$t_i=4(i=1,\cdots,9),t_{10}=16$$

です。
よって残るカードの枚数は、

$$r(h,d;t):=(r_i)_{i=1,\cdots,10}$$

で表すことができます。もちろん、

$$r_i(h,d;t)=t_i-h_i-d_i$$

です。
これを用いると、選ばれるカードの確率は、以下の

$$p(h,d;t):=(p_i)_{i=1,\cdots,10}$$

で表すことができます。

$$p_i(h,d;t)=\frac{r_i}{\sum_{j=1}^{10}(r_j)}$$

つまり、プレイヤーとディーラーの手札ベクトルが、ある遷移確率によって変化するという状況です。
ある状態から次の状態への遷移確率は、

$$p((h,d;t)\to(h^i,d;t)) = p_i(h,d;t)$$
$$p((h,d;t)\to(h,d^i;t)) = p_i(h,d;t)$$

によって与えられます。ここで、

$$h^i_j = h_j + \delta_{ij}$$
$$d^i_j = d_j + \delta_{ij}$$

です。ただし、$\delta_{ij}$はクロネッカーのデルタです。

### ディーラーの伏せられたカード
ディーラーの伏せられている手札については情報がありません。
なので、初期状態においてはディーラーは1枚しかカードを持っていないとして処理します。

### バースト
続いてバーストについてです。
これは手札の合計値によって決まります。

このために、まずは手札$u$の合計値を計算する関数

$$s(u) = s'(u) + s''(u)$$
$$s'(u) = \sum_{i=1}^{10}iu_i$$
$$s''(u) = a * \min\left(floor\left(\frac{b - s'(u)}{a}\right), u_1\right)$$
$$a = 10$$
$$b = 21$$

を定義します。第2項で、Aの扱いを考慮しています。
これにより、プレイヤーとディーラーそれぞれのバーストを判定する関数

```math
x(u) = \left\{
\begin{array}{cc}
1 & \left( s(u) > b \right) \\
0 & \left( s(u) \leq b \right)
\end{array}
\right.
```

を定義します。

### 勝敗
さて、次は勝敗についてです。
これは、プレイヤーとディーラーの手札の合計値のみによって決まります。
よって、処理が継続しているかどうかにかかわらず、手札に対して勝敗を決定する関数

$$w(h,d)$$

を定義することが可能です。
処理の順序から、プレイヤーがバーストした手札では即座に負けとなります。

$$w(h,d) = 0 \left(x(h) = 1\right)$$

ディーラーのみのバーストについては、勝利となります。

$$w(h,d) = 1 \left(x(d) = 1 \wedge x(h) = 0\right)$$

それ以外の場合は、プレイヤーの手札の合計値がディーラーの手札の合計値よりも大きい場合に勝利となります。

$$w(h,d) = 1 \left(s(h) > s(d) \wedge x(h) = 0 \wedge x(d) = 0\right)$$
$$w(h,d) = 0 \left(s(h) \leq s(d) \wedge x(h) = 0 \wedge x(d) = 0\right)$$

ここで、手札の合計値が等しい場合にはプレイヤーの勝利としています。理由は、この後の方策の定式化を楽にするためです。

### ターン
続いてターンについてです。
プレイヤーのターンとディーラーのターンを定義する変数

$$c$$

を導入しておきます。プレイヤーのターンには0、ディーラーのターンには1を割り当てます。
初期状態は$c=0$とし、スタンドが選択された場合には$c$に１を加えます。$c>1$になった場合には、勝敗を決定します。

### 方策
続いて方策についてです。
方策は、プレイヤーの手札とディーラーの1枚目のカードに対して、カードを引くかスタンドするかを決定する関数

$$\pi(h,d;c)$$

です。カードを引くを1、スタンドを0とします。
方策のうちのいくつかは最適化の対象外です。
まず、バーストの場合にはスタンドを選択するようにします。

$$\pi(h,d;c) = 0 \left(x(h) = 1\right)$$
$$\pi(h,d;c) = 0 \left(x(d) = 1\right)$$

また、ディーラーのターンにおいては、ディーラーの手札がプレイヤーの手札よりも大きくなるまでカードを引くようにします。

$$\pi(h,d;1) = w(h,d) \left(x(h) = 0 \wedge x(d) = 0 \right)$$

残るはプレイヤーのターンにおいて、カードを引くかスタンドするかを決定する方策です。
この部分が最適化の対象となります。

### 状態価値関数

最後に、状態価値関数

$$v_\pi(h,d;c)$$

を定義します。
これは、プレイヤーとディーラーの手札が与えられたときに、方策$\pi$に従って行動した場合の期待値です。
この関数は、以下の再帰的な関係式によって定義されます。

$$v_\pi(h,d;0) = v_\pi(h,d;1) \left(\pi(h,d;0) = 0 \right)$$
$$v_\pi(h,d;1) = w(h,d) \left(\pi(h,d;1) = 0 \right)$$
$$v_\pi(h,d;0) = \sum_{i=1,\cdots,10} p_i(h,d;t) v_\pi(h^i,d;0) \left(\pi(h,d;0) = 1 \right)$$
$$v_\pi(h,d;1) = \sum_{i=1,\cdots,10} p_i(h,d;t) v_\pi(h,d^i;1) \left(\pi(h,d;1) = 1 \right)$$

最適化は、この状態価値関数を最大化するような方策を求めることになります。
すなわち、方策の評価としては、以下のような関数

$$y(\pi) = \sum_{h,d} v_\pi(h,d;0)$$

を最大化することになります。

### 最適化対象の削減

まず今回、異なる状態での状態価値関数の値はトレードオフの関係にありません。
したがって、それぞれの状態の状態価値関数の値を独立に最適化することで、全体の状態価値関数の値を最大化することができます。
ただし、$h,d$はそれぞれ10次元のベクトルであり、$4^9\times16$通りもの状態が存在するため、全探索は現実的ではありません。
実際、バーストした手札からしか遷移できないような手札については、そもそも考慮する必要がありません。
したがって、考慮すべき状態は、手札の合計値が30以下のものに限定します。すなわち、

$$s(h) < b + 10$$
$$s(d) < b + 10$$

です。また、プレイヤーのターンの間は、ディーラーの手札は１枚しかありません。よって

$$n(d) = 1 \left(c = 0\right)$$
$$n(u) := \sum_{i=1,\cdots,10}u_i$$

として、状態空間をさらに削減します。
さらに、プレイヤーがバーストした場合には、ディーラーはカードを引きません。よってこの場合、ディーラーの手札が２枚以上の状況は考慮外とします。

$$n(d) = 1 \left(x(h) = 1\right)$$

最後に、両者のカード枚数の合計はtを超えることはありません。よって

$$h_i + d_i \leq t_i$$

以上の条件を満たす状態について、状態価値関数を求めます。
ここからさらに、バースト状態とディーラーのターンに対して与えられている方策をはめ込み、残った状態についてのみ最適な方策を求めることになります。

### 線形性について

最後に、この問題が線形であるかどうかについて考えます。

状態価値関数は、再帰的な関係式によって定義されています。
状態価値関数の計算においては、ほかの状態における方策を用いています。
したがって、ある状態$(h_0,d_0)$での方策$\pi$を変更した場合、その状態に遷移できる状態$(h_1,d_1)$における状態価値関数の値にも影響を及ぼし得ます。
"及ぼし得る"と書いているのは、たとえばこの遷移がカードを引くことにより生じるもので、$(h_1,d_1)$での方策がスタンドであった場合には、$(h_0,d_0)$での方策の変更が影響を及ぼさないことになります。
もちろん、$(h_1,d_1)$での方策がカードを引くことであった場合には、$(h_0,d_0)$での方策の変更が影響を及ぼすことになります。
すなわち、ある状態での方策を変更した場合に評価関数に及ぶ影響が、ほかの状態での方策に影響されているわけです。

したがって、この問題は線形ではありません。



## 実装
次に、定式化した数理最適化問題を実装します。
と思ったのですが、想像より遥かに重くなってしまったので、今回はここまでとします。
次回の記事で、実装について詳しく書く予定です。