text:
# Generating Wick's theorem contractions

## Preliminaries

Let's start by importing wick&d and defining a Slater determinant reference

```python
import wicked as w
from IPython.display import display, Math, Latex

def latex(expr):
    """Function to render any object that has a member latex() function"""
    display(Math(expr.latex()))

w.reset_space()
w.add_space("o", "fermion", "occupied", ['i','j','k','l','m'])
w.add_space("v", "fermion", "unoccupied", ['a','b','c','d','e','f'])
```

## Wick's theorem contractions

Consider the product of operators
\begin{equation}
\hat{F} \hat{T}_1 =
\left( \sum_{pq}^\mathbb{O \cup V} f_{p}^{q} \{ \hat{a}^{p}_{q} \} \right)
\left( \sum_{i}^\mathbb{O} \sum_{a}^\mathbb{V} t_{a}^{i} \{ \hat{a}^{a}_{i} \} \right)
\end{equation}
which we can create with the following code:

```python
F = w.utils.gen_op('f',1,'ov','ov')
T1 = w.op('t',['v+ o'])
FT1 = F @ T1
print(FT1)
```

    + f { v+ v } t { v+ o }
    + f { v+ o } t { v+ o }
    + f { o+ v } t { v+ o }
    + f { o+ o } t { v+ o }

To apply Wick's theorem to this product, you need to create a WickTheorem object. You can compute the contraction with the `contraction` function.
This function takes as arguments the operator expression to contract and the minimum and maximum number of second quantized operators in the final expression (`minrank/maxrank`). Here we request all terms with zero to four second quantized operators:

```python
wt = w.WickTheorem()
expr = wt.contract(FT1,minrank=0,maxrank=4)
expr
```

    f^{o1}_{o0} t^{o2}_{v0} a+(o0) a+(v0) a-(o2) a-(o1)
    +f^{v0}_{o0} t^{o0}_{v0}
    +f^{v0}_{o0} t^{o1}_{v0} a+(o0) a-(o1)
    -f^{v0}_{o0} t^{o0}_{v1} a+(v1) a-(v0)
    -f^{v0}_{o0} t^{o1}_{v1} a+(o0) a+(v1) a-(v0) a-(o1)
    -f^{o0}_{o1} t^{o1}_{v0} a+(v0) a-(o0)
    +f^{o0}_{v0} t^{o1}_{v1} a+(v0) a+(v1) a-(o1) a-(o0)
    +f^{v1}_{v0} t^{o0}_{v1} a+(v0) a-(o0)
    -f^{v1}_{v0} t^{o0}_{v2} a+(v0) a+(v2) a-(v1) a-(o0)

You can get a LaTeX expression with the `latex()` function

```python
print(expr.latex())
```

    + {f}^{j}_{i} {t}^{k}_{a} \hat{a}^{i} \hat{a}^{a} \hat{a}_{k} \hat{a}_{j} \\
    + {f}^{a}_{i} {t}^{i}_{a} \\
    + {f}^{a}_{i} {t}^{j}_{a} \hat{a}^{i} \hat{a}_{j} \\
    - {f}^{a}_{i} {t}^{i}_{b} \hat{a}^{b} \hat{a}_{a} \\
    - {f}^{a}_{i} {t}^{j}_{b} \hat{a}^{i} \hat{a}^{b} \hat{a}_{a} \hat{a}_{j} \\
    - {f}^{i}_{j} {t}^{j}_{a} \hat{a}^{a} \hat{a}_{i} \\
    + {f}^{i}_{a} {t}^{j}_{b} \hat{a}^{a} \hat{a}^{b} \hat{a}_{j} \hat{a}_{i} \\
    + {f}^{b}_{a} {t}^{i}_{b} \hat{a}^{a} \hat{a}_{i} \\
    - {f}^{b}_{a} {t}^{i}_{c} \hat{a}^{a} \hat{a}^{c} \hat{a}_{b} \hat{a}_{i}

Or render the expression directly in a notebook

```python
latex(expr)
```

$\displaystyle + {f}^{j}_{i} {t}^{k}_{a} \hat{a}^{i} \hat{a}^{a} \hat{a}_{k} \hat{a}_{j} \\
+ {f}^{a}_{i} {t}^{i}_{a} \\
+ {f}^{a}_{i} {t}^{j}_{a} \hat{a}^{i} \hat{a}_{j} \\
- {f}^{a}_{i} {t}^{i}_{b} \hat{a}^{b} \hat{a}_{a} \\
- {f}^{a}_{i} {t}^{j}_{b} \hat{a}^{i} \hat{a}^{b} \hat{a}_{a} \hat{a}_{j} \\
- {f}^{i}_{j} {t}^{j}_{a} \hat{a}^{a} \hat{a}_{i} \\
+ {f}^{i}_{a} {t}^{j}_{b} \hat{a}^{a} \hat{a}^{b} \hat{a}_{j} \hat{a}_{i} \\
+ {f}^{b}_{a} {t}^{i}_{b} \hat{a}^{a} \hat{a}_{i} \\
- {f}^{b}_{a} {t}^{i}_{c} \hat{a}^{a} \hat{a}^{c} \hat{a}_{b} \hat{a}_{i}$

Let's evaluate a more complex expression
\begin{equation}
\langle \Phi | \hat{\Lambda}_1 [\hat{F}, \hat{T}_1] | \Phi \rangle
\end{equation}
where
\begin{equation}
\hat{\Lambda}_1 =
\sum_{i}^\mathbb{O} \sum_{a}^\mathbb{V} \lambda_{i}^{a} \{ \hat{a}^{i}_{a} \}
\end{equation}

```python
L1 = w.op('lambda',['o+ v'])
opexpr = L1 @ w.commutator(F,T1)
expr = wt.contract(opexpr,0,0)
latex(expr)
```

$\displaystyle - {f}^{j}_{i} {\lambda}^{a}_{j} {t}^{i}_{a} \\
+ {f}^{b}_{a} {\lambda}^{a}_{i} {t}^{i}_{b}$

## Example: Coupled cluster energy expression

Here we compute the coupled cluster energy expression
\begin{equation}
E = \langle \Phi | e^{-\hat{T}} \hat{H} e^{\hat{T}} | \Phi \rangle
= E_0 + \sum_{i}^\mathbb{O} \sum_{a}^\mathbb{V} f^{a}_{i} t^{i}_{a} +
\frac{1}{4} \sum_{ij}^\mathbb{O} \sum_{ab}^\mathbb{V}
 (t^{i j}_{a b} + 2 t^{i}_{a} t^{j}_{b}) v^{a b}_{i j}
\end{equation}
We first define the full Hamiltonian operator and the cluster operator

```python
E0 = w.op("E_0",[""])
F = w.utils.gen_op('f',1,'ov','ov')
V = w.utils.gen_op('v',2,'ov','ov')
H = E0 + F + V

T = w.op("t",["v+ o", "v+ v+ o o", "v+ v+ v+ o o o"])
```

Next, we compute the operator $e^{-\hat{T}} \hat{H} e^{\hat{T}}$ up to the term quadratic in $\hat{T}$ and compute the fully contracted term

```python
Hbar = w.bch_series(H,T,2)
expr = wt.contract(Hbar,0,0)
latex(expr)
```

$\displaystyle + {E_0}^{}_{} \\
+ {f}^{a}_{i} {t}^{i}_{a} \\
+\frac{1}{2} {t}^{i}_{a} {t}^{j}_{b} {v}^{a b}_{i j} \\
+\frac{1}{4} {t}^{i j}_{a b} {v}^{a b}_{i j}$

max_stars_repo_name: fevangelista/pyWicked

max_issues_repo_name: fevangelista/pyWicked

max_forks_repo_name: fevangelista/pyWicked


text:
## 前準備（万有引力の働く系）

万有引力が働く系
\begin{equation}
    m\ddot{\vec{r}}=-G\frac{Mm}{r^2}\vec{e}_r
    \Leftrightarrow
    \ddot{\vec{r}}=-\frac{GM}{r^2}\vec{e}_r
\end{equation}
を考える。
簡便のため $GM=1$ となるように単位系を取る。
このとき系を直交座標系に書き換えると
\begin{equation}
    \ddot{x}=-\frac{x}{(x^2+y^2)^{5/2}},\quad
    \ddot{y}=-\frac{y}{(x^2+y^2)^{5/2}}
\end{equation}
となる。
更にハミルトン系へと書き換えると
\begin{align}
    &\dot{q}_x=p_x,\quad
    \dot{p}_x=-\frac{q_x}{(q_x^2+q_y^2)^{3/2}},\\
    &\dot{q}_y=p_y,\quad
    \dot{p}_y=-\frac{q_y}{(q_x^2+q_y^2)^{3/2}}
\end{align}
となる。
エネルギーは運動の積分で
\begin{equation}
    E=\frac{1}{2}(p_x^2+p_y^2)-\frac{1}{(q_x^2+q_y^2)^{1/2}}
\end{equation}
である。
初期条件に $(q_x,q_y,p_x,p_y)=(1,0,0,1)$ を選ぶと
$E = -1/2$ である。

```
from typing import List, Callable
from numpy import ndarray

def dqxdt(x: ndarray) -> float:
    return x[2]

def dqydt(x: ndarray) -> float:
    return x[3]

def dpxdt(x: ndarray) -> float:
    return -x[0]/(x[0]**2 + x[1]**2)**1.5

def dpydt(x: ndarray) -> float:
    return -x[1]/(x[0]**2 + x[1]**2)**1.5

# 運動方程式
dxdt: List[Callable[[ndarray], float]] = [dqxdt,dqydt,dpxdt,dpydt]
```

## 一段階法（陽的解法）

一段階法は実装が容易。ただし非シンプレクティックで硬い方程式系に不適。

### オイラー法

オイラー法は単純な差分近似
\begin{equation}
    \vec{x}(t+dt)
    \simeq\vec{x}(t)+\frac{d\vec{x}}{dt}(t,\vec{x}(t))dt
\end{equation}

```
import math
import numpy as np

def euler(dxdt: List[Callable[[ndarray], float]],
          init: ndarray, t_end: float) -> ndarray:
    # 次元チェック
    if len(dxdt) != len(init):
        raise Exception("Dimension error!")

    # ソルバーのパラメータ
    dt: float = 1e-2

    # 初期状態
    x_list: ndarray = np.array([init])
    e_list: ndarray = np.array([(init[2]**2+init[3]**2)/2
                                -(init[0]**2+init[1]**2)**(-1/2)])
    # ソルバー実施
    steps: List[int] = range(1, math.floor(t_end/dt), 1)
    for t in steps:
        # オイラー法
        x_next: ndarray = np.zeros(len(init))
        for k in range(len(init)):
            x_next[k] = x_list[-1][k] + dxdt[k](x_list[-1]) * dt
        # 結果を追加
        x_list = np.vstack((x_list,x_next))
        e_list = np.append(e_list,
            (x_next[2]**2+x_next[3]**2)/2-(x_next[0]**2+x_next[1]**2)**(-1/2))

    # 位置・運動量・エネルギーを戻り値として戻す
    return np.vstack((x_list.T, e_list))
```

軌道を計算して表示。誤差によって軌道が広がる。

```
import matplotlib.pyplot as plt

ans: ndarray = euler(dxdt, np.array([1,0,0,1]), 100)
plt.plot(ans[0], ans[1])
```

誤差でエネルギーが上昇。

```
plt.plot(ans[4])
```

### ホイン法

ホイン法は次数 2 の 1 段階法。
\begin{equation}
\vec{k}_1\equiv\frac{d\vec{x}}{dt}(t,\vec{x}(t)),\quad
\vec{k}_2\equiv\frac{d\vec{x}}{dt}\left(t+dt,\vec{x}(t)+\vec{k}_1dt\right)
\end{equation}
とすれば
\begin{equation}
    \vec{x}(t+dt)
    \simeq\vec{x}(t)
    +(\vec{k}_1+\vec{k}_2)\frac{dt}{2}
\end{equation}

```
def heun(dxdt: List[Callable[[ndarray], float]],
          init: ndarray, t_end: float) -> ndarray:
    # 次元チェック
    if len(dxdt) != len(init):
        raise Exception("Dimension error!")

    # ソルバーのパラメータ
    dt: float = 1e-2
    steps: List[int] = range(1, math.floor(t_end/dt), 1)

    # 初期状態
    x_list: ndarray = np.array([init])
    e_list: ndarray = np.array([(init[2]**2+init[3]**2)/2
                                -(init[0]**2+init[1]**2)**(-1/2)])
    # ソルバー実施
    for t in steps:
        # ホイン法
        k: ndarray = np.zeros(len(init))
        for d in range(len(init)):
            k[d] = x_list[-1][d] + dxdt[d](x_list[-1]) * dt

        x_next: ndarray = np.zeros(len(init))
        for d in range(len(init)):
            x_next[d] = x_list[-1][d] + dt/2 * (dxdt[d](x_list[-1])+dxdt[d](k))

        # 結果を追加
        x_list = np.vstack((x_list,x_next))
        e_list = np.append(e_list,
            (x_next[2]**2+x_next[3]**2)/2-(x_next[0]**2+x_next[1]**2)**(-1/2))

    # 位置・運動量・エネルギーを戻り値として戻す
    return np.vstack((x_list.T, e_list))
```

```
ans: ndarray = heun(dxdt, [1,0,0,1], 100)
plt.plot(ans[0], ans[1])
```

オイラー法と比較して誤差が低減する

```
plt.ticklabel_format(useOffset=False)
plt.plot(ans[4])
```

### ルンゲ・クッタ法

ルンゲ・クッタ法は次数 4 の 1 段階法。
\begin{align}
    &\vec{k}_1\equiv\frac{d\vec{x}}{dt}\left(t,\vec{x}\right),\\
    &\vec{k}_2\equiv\frac{d\vec{x}}{dt}\left(t+\frac{dt}{2},\vec{x}+\vec{k}_1\frac{dt}{2}\right),\\
    &\vec{k}_3\equiv\frac{d\vec{x}}{dt}\left(t+\frac{dt}{2},\vec{x}+\vec{k}_2\frac{dt}{2}\right),\\
    &\vec{k}_4\equiv\frac{d\vec{x}}{dt}\left(t+dt,\vec{x}+\vec{k}_3dt\right),
\end{align}
として
\begin{equation}
    \vec{x}(t+dt)
    \simeq\vec{x}(t)+(\vec{k}_1+2\vec{k}_2+2\vec{k}_3+\vec{k}_4)\frac{dt}{6}
\end{equation}

```
def rk4(dxdt: List[Callable[[ndarray], float]],
        init: ndarray, t_end: float) -> ndarray:
    # 次元チェック
    if len(dxdt) != len(init):
        raise Exception("Dimension error!")

    # ソルバーのパラメータ
    dt: float = 1e-2
    steps: List[int] = range(1, math.floor(t_end/dt), 1)

    # 初期状態
    x_list: ndarray = np.array([init])
    e_list: ndarray = np.array([(init[2]**2+init[3]**2)/2
                                -(init[0]**2+init[1]**2)**(-1/2)])
    # ソルバー実施
    for t in steps:
        # ルンゲ・クッタ法
        k: ndarray = np.zeros((4, len(init)))
        for d in range(len(init)):
            k[0,d] = dxdt[d](x_list[-1])
        for d in range(len(init)):
            k[1,d] = dxdt[d](x_list[-1]+dt/2*k[0])
        for d in range(len(init)):
            k[2,d] = dxdt[d](x_list[-1]+dt/2*k[1])
        for d in range(len(init)):
            k[3,d] = dxdt[d](x_list[-1]+dt*k[2])

        x_next: ndarray = np.zeros(len(init))
        for d in range(len(init)):
            x_next[d] = x_list[-1][d] + dt * (k[0,d]+2*k[1,d]+2*k[2,d]+k[3,d])/6

        # 結果を追加
        x_list = np.vstack((x_list,x_next))
        e_list = np.append(e_list,
            (x_next[2]**2+x_next[3]**2)/2-(x_next[0]**2+x_next[1]**2)**(-1/2))

    # 位置・運動量・エネルギーを戻り値として戻す
    return np.vstack((x_list.T, e_list))
```

```
ans: ndarray = rk4(dxdt, [1,0,0,1], 100)
plt.plot(ans[0], ans[1])
```

誤差が更に低減するが、誤差の増大方向が逆

```
plt.ticklabel_format(useOffset=False)
plt.plot(ans[4])
```

## 線形多段法

線形多段法は一段階法と比べて導出が容易で計算量が少ない。
一方で計算開始前に一段階法などで数回計算を行う必要がある。
また一段階法と同様に非シンプレクティックで硬い方程式に不適。

### アダムス・バッシュフォース法（陽的解法）

アダムス・バッシュフォース法は補外による陽的線形多段法。4次の公式は
\begin{equation}
    \vec{x}(t+dt)
    \simeq \vec{x}(t)+\left(-\frac{9}{24}\vec{x}'(t-3dt)+\frac{37}{24}\vec{x}'(t-2dt)-\frac{59}{24}\vec{x}'(t-dt)+\frac{55}{24}\vec{x}'(t)\right)dt
\end{equation}

```
def adams_bashforth(dxdt: List[Callable[[ndarray], float]],
                    init: ndarray, t_end: float) -> ndarray:
    # 次元チェック
    if len(dxdt) != len(init):
        raise Exception("Dimension error!")

    # ソルバーのパラメータ
    dt: float = 1e-2

    # 初期状態
    x_list: ndarray = np.array([init])
    e_list: ndarray = np.array([(init[2]**2+init[3]**2)/2
                                -(init[0]**2+init[1]**2)**(-1/2)])
    # 多段法の準備
    for t in range(1,4):
        # オイラー法
        x_next: ndarray = np.zeros(len(init))
        for k in range(len(init)):
            x_next[k] = x_list[-1][k] + dxdt[k](x_list[-1]) * dt
        # 結果を追加
        x_list = np.vstack((x_list,x_next))
        e_list = np.append(e_list,
            (x_next[2]**2+x_next[3]**2)/2-(x_next[0]**2+x_next[1]**2)**(-1/2))

    # ソルバー実施
    steps: List[int] = range(5, math.floor(t_end/dt), 1)
    for t in steps:
        # アダムス・バッシュフォース法
        x_next: ndarray = np.zeros(len(init))
        for k in range(len(init)):
            x_next[k] = x_list[-1][k] + (
                - 9 * dxdt[k](x_list[-4])
                +37 * dxdt[k](x_list[-3])
                -59 * dxdt[k](x_list[-2])
                +55 * dxdt[k](x_list[-1])) * dt/24

        # 結果を追加
        x_list = np.vstack((x_list,x_next))
        e_list = np.append(e_list,
            (x_next[2]**2+x_next[3]**2)/2-(x_next[0]**2+x_next[1]**2)**(-1/2))

    # 位置・運動量・エネルギーを戻り値として戻す
    return np.vstack((x_list.T, e_list))
```

```
ans: ndarray = adams_bashforth(dxdt, np.array([1,0,0,1]), 100)
plt.plot(ans[0], ans[1])
```

前準備を除いたエネルギーをプロット

```
plt.ticklabel_format(useOffset=False)
plt.plot(ans[4,10:])
```

### アダムス・ムルトン法（陰的解法）

アダムス・ムルトン法は補間による陰的線形多段法。4次の公式は
\begin{equation}
    \vec{x}(t+dt)
    \simeq \vec{x}(t)+\left(\frac{1}{24}\vec{x}'(t-2dt)-\frac{5}{24}\vec{x}'(t-dt)-\frac{19}{24}\vec{x}'(t)+\frac{9}{24}\vec{x}'(t+dt)\right)dt.
\end{equation}

### アダムス法（予測子修正子法）

アダムス法（予測子修正子法）は以下の通り：
1. アダムス・バッシュフォース法（予測子）で $\vec{x}(t)$ から $\tilde{\vec{x}}(t+dt)$ を計算
2. $\tilde{\vec{x}}(t+dt)$ から $\vec{x}'(t+dt)$ を計算
3. アダムス・ムルトン法（修正子）で $\tilde{\vec{x}}(t+dt)$ を $\vec{x}(t+dt)$に修正
4. 修正の差が一定値以下なら終了.そうでないなら 2. に戻って修正結果を再修正.

```
def adams(dxdt: List[Callable[[ndarray], float]],
                    init: ndarray, t_end: float) -> ndarray:
    # 次元チェック
    if len(dxdt) != len(init):
        raise Exception("Dimension error!")

    # ソルバーのパラメータ
    dt: float = 1e-2
    iter_max: int = 10

    # 初期状態
    x_list: ndarray = np.array([init])
    e_list: ndarray = np.array([(init[2]**2+init[3]**2)/2
                                -(init[0]**2+init[1]**2)**(-1/2)])
    # 多段法の準備
    for t in range(1,4):
        # オイラー法
        x_next: ndarray = np.zeros(len(init))
        for k in range(len(init)):
            x_next[k] = x_list[-1][k] + dxdt[k](x_list[-1]) * dt
        # 結果を追加
        x_list = np.vstack((x_list,x_next))
        e_list = np.append(e_list,
            (x_next[2]**2+x_next[3]**2)/2-(x_next[0]**2+x_next[1]**2)**(-1/2))

    # ソルバー実施
    steps: List[int] = range(5, math.floor(t_end/dt), 1)
    for t in steps:
        # アダムス・バッシュフォース法
        x_next: ndarray = np.zeros(len(init))
        for k in range(len(init)):
            x_next[k] = x_list[-1][k] + (
                - 9 * dxdt[k](x_list[-4])
                +37 * dxdt[k](x_list[-3])
                -59 * dxdt[k](x_list[-2])
                +55 * dxdt[k](x_list[-1])) * dt/24

        for i in range(0,iter_max):
            # 微分値を計算
            x_prime: ndarray = np.zeros(len(init))
            for k in range(len(init)):
                x_prime[k] = dxdt[k](x_next)

            # アダムス・ムルトン法
            x_nnext: ndarray = np.zeros(len(init))
            for k in range(len(init)):
                x_nnext[k] = x_list[-1][k] + (
                          dxdt[k](x_list[-3])
                    - 5 * dxdt[k](x_list[-2])
                    +19 * dxdt[k](x_list[-1])
                    + 9 * x_prime[k]) * dt/24

            diff: float = np.linalg.norm(x_next-x_nnext, ord=2)
            x_next = x_nnext

            if diff < 1e-10:
                break

        # 結果を追加
        x_list = np.vstack((x_list,x_next))
        e_list = np.append(e_list,
            (x_next[2]**2+x_next[3]**2)/2-(x_next[0]**2+x_next[1]**2)**(-1/2))

    # 位置・運動量・エネルギーを戻り値として戻す
    return np.vstack((x_list.T, e_list))
```

```
ans: ndarray = adams(dxdt, np.array([1,0,0,1]), 100)
plt.plot(ans[0], ans[1])
```

```
plt.ticklabel_format(useOffset=False)
plt.plot(ans[4,10:])
```

## シンプレクティック数値積分

ハミルトン系の保存量を近似的に維持する。
陽公式はリープフロッグ法をベースとした合成公式。
陰公式はガウス・ルジャンドル積分を利用。

### リープ・フロッグ法（陽的解法）

リープ・フロッグ法は中点差分による2段階法で2次のシンプレクティック数値積分。
\begin{align}
    &\vec{q}(t+2dt)\simeq\vec{q}(t)+\vec{p}(t+dt)\cdot2dt,\\
    &\vec{p}(t+2dt)\simeq\vec{p}(t)+\frac{d\vec{p}}{dt}(t+dt)\cdot2dt
\end{align}
分離可能なハミルトン系では
\begin{align}
    &\vec{q}(t+dt/2)\simeq\vec{q}(t)+\frac{d\vec{q}}{dt}(p(t))dt/2,\\
    &\vec{p}(t+dt)\simeq\vec{p}(t)+\frac{d\vec{p}}{dt}(q(t+dt/2))dt,\\
    &\vec{q}(t+dt)\simeq\vec{q}(t+dt/2)+\frac{d\vec{q}}{dt}(p(t+dt))dt/2
\end{align}
となる。

```
def leap_flog(dxdt: List[Callable[[ndarray], float]],
          init: ndarray, t_end: float) -> ndarray:
    # 次元チェック
    if len(dxdt) != len(init):
        raise Exception("Dimension error!")

    # ソルバーのパラメータ
    dt: float = 1e-2
    steps: List[int] = range(1, math.floor(t_end/dt), 1)

    # 初期状態
    x_list: ndarray = np.array([init])
    e_list: ndarray = np.array([(init[2]**2+init[3]**2)/2
                                -(init[0]**2+init[1]**2)**(-1/2)])
    # ソルバー実施
    for t in steps:
        # リープ・フロッグ法
        x_next: ndarray = np.zeros(len(init))

        x_next[0] = x_list[-1,0]+dxdt[0](x_list[-1])*dt/2
        x_next[1] = x_list[-1,1]+dxdt[1](x_list[-1])*dt/2

        x_next[2] = x_list[-1,2]+dxdt[2](x_next)*dt
        x_next[3] = x_list[-1,3]+dxdt[3](x_next)*dt

        x_next[0] = x_next[0]+dxdt[0](x_next)*dt/2
        x_next[1] = x_next[1]+dxdt[1](x_next)*dt/2

        # 結果を追加
        x_list = np.vstack((x_list,x_next))
        e_list = np.append(e_list,
            (x_next[2]**2+x_next[3]**2)/2-(x_next[0]**2+x_next[1]**2)**(-1/2))

    # 位置・運動量・エネルギーを戻り値として戻す
    return np.vstack((x_list.T, e_list))
```

```
ans: ndarray = leap_flog(dxdt, [1,0,0,1], 100)
plt.plot(ans[0], ans[1])
```

エネルギーの誤差は振動。エネルギーが近似的に保存される。

```
plt.ticklabel_format(useOffset=False)
plt.plot(ans[4])
```

ext: ipynb

max_stars_repo_name: applejxd/colaboratory

max_issues_repo_name: applejxd/colaboratory

max_forks_repo_name: applejxd/colaboratory


text:
**Author: Mack Sweeney**

This document is a review of the literature on methods for addressing contextual, non-stationary bandit problems. We'll start with a synopsis, which proposes extensions and a plan for developing and evaluating them, and then concludes with the bibliography. The rest of the document contains notes on the various sources.

# Synopsis

The standard bandit problem assumes rewards are stationary, the value of actions are based only on the actions themselves, actions are independent, and feedback is immediate. The contextual bandit formulation introduces the notion that we also observe additional context about each data point and the value of actions is dependent on the observed context vector. The non-stationary bandit problem drops the assumption that the underlying reward distributions do not change over time. There have been two main varieties of non-stationarity studied: piecewise stationary rewards (sometimes called switching bandits) and drifting rewards (sometimes called restless bandits).

A variety of methods have been studied for various combinations of non-stationarity and contextual settings. Methods for handling non-stationarity generally fall into one of three categories:

1.  discount/forget old data,
2.  learn only from data in a sliding window, with the sliding window size either fixed or adaptive, and
3.  detect changepoints or accumulated error thresholds and restart in some way.

These methods can usually be mixed and matched with the basic bandit design principles and adaptations of those principles for contextual settings. For instance, Thompson Sampling (TS) and Upper-Confidence Bound (UCB) principles have been applied to solve many simple bandit problems. However, they have also both been extended for the contextual setting via linear model formulations. Then these linear model formulations have been extended for the non-stationary setting by discounting, sliding window, and changepoint detection augmentations.

## Limitations and Opportunities

While there is a robust body of literature in this area, there are still very promising opportunities for improvement at the intersection of non-stationary and contextual. These are summarized here:

1.  **Discard data more judiciously.** While many of these existing methods are flexible and can be applied to a variety of problems, this flexibility often comes at the cost of discarding data indiscriminately or based on parameters that must be tuned. Further, few of these methods take advantage of the information structure of the action space to learn as efficiently as possible.
2.  **Experiment with feedback delay.** The vast majority of the bandit literature focuses on problems with immediate reward feedback. This neglects problems where optimizing for delayed feedback is critical to making good decisions. For instance, many marketers are moving away from optimizing for click-through rates towards deeper performance metrics such as conversions. The few papers studying delayed feedback problems indicate huge differences in efficacy between methods that otherwise look very similar in terms of empirical performance and theoretical bounds.
3.  **Experiment with finite-horizons.** Much of the bandit literature is devoted to regret analysis in an infinite-horizon setting. However, some authors have studied finite-time horizons and found that many methods that do well asymptotically may do very poorly in earlier time points. Most marketing-based experimentation is concerned with iterating as rapidly as possible, introducing new actions and removing older ones as information is gathered and lessons learned. So many of these theoretical studies are called into question for these sorts of problems.

There are several promising design principles to apply to these problems:

1.  Recently introduced bandit design principles, such as Information-Directed Sampling (IDS), propose to not just minimize regret, but also optimize for information gain. This leads to more efficient exploration of the action space even in the stationary setting. In non-stationary problems, the potential is even greater due to the additional uncertainty introduced by moving reward targets. This is a similar but perhaps more principled version of the notion of "exploration bonuses" which has been successful in other methods such as POKER.
2.  Bayesian methods for filtering and smoothing have been under-utilized for the non-stationary problem. The one application of Kalman Filtering to a stationary setting showed great promise for handling both restless and switching bandit cases without dependence on many sensitive hyperparameters. Old data is only discounted if there is evidence of high transition noise, and the method can separate observation noise from transition noise and handle changepoints.
3.  An early study of TS and UCB under delayed feedback showed that randomized policies such as TS handle delayed feedback significantly better. This is increasingly true as the delay increases.

## Proposed Methodological Improvement

These observations motivate a new approach: **implement a Bayesian hierarchial generalized linear model (GLM) using Kalman Filtering to estimate action rewards for contextual bandit problems.** Then combine principles from TS and IDS to optimize a utility function that combines expected gains and information gain to arrive at a randomized probability matching policy. This method has the added advantage of the interpretability provided by fully Bayesian methods, which is critical for applications where the actions are provided by humans. This is true of many marketing scenarios, such as UX optimization and testing of new products and incentive offerings.

### Incremental Development Plan

The proposed method involves several improvements which can be implemented in a sequence which will produce a variety of increasingly flexible methods:

1.  Extend PG-TS using a hierarchical model. The current implementation uses a simple non-hierarchical LR formulation that loses the opportunity to share information across predictors. For now, let's call this H-PG-TS.
2.  Apply the simple discounting procedure used in discounted TS to the PG-TS and H-PG-TS methods to extend them to a non-stationary setting: dPG-TS and dH-PG-TS.
3.  Incorporate a Kalman Filter on the coefficients of both methods to demonstrate superior performance relative to manually tuned discount factors: KF-H-PG-TS and KF-PG-TS.
4.  Adjust policy from TS to IDS by incorporating a measure of info gain into the policy calculation: will come up with better names later.  : )

### Evaluation Plan

Compare the developed methods to each other as well as several existing methods. Consider $\epsilon$-greedy, Bayes UCB, TS, dUCB, dTS, dLinUCB, dLinTS, linear bandit IDS, KG, and $\text{KG}^*$. Use simulated data and consider principles outlined in the following section.

## Experiments

Many authors do not actually conduct experiments, preferring to derive theoretical regret bounds. Of those who do conduct experiments, most use simulated data. A few do use real datasets to compare methods side by side. In general, it seems that simulated experiments are a useful approach for bandit algorithms. The most valuable lessons learned from the surveyed literature are:

1.  Experiment with a variety of action spaces. For non-contextual, use both few and many arms. For contextual, also use varying dimensions of context vectors.
2.  Be clear about assumed reward distributions and how those distributions are being parameterized to generate synthetic data (e.g. Gaussian with params 0, 1).
3.  Also run experiments for a variety of lengths, e.g. few iterations and many iterations. This can reveal trends not at all apparent from asymptotic guarantees many authors rely on. Also consider the relation of iterations and number of arms. Vermorel and Mori looked at cases corresponding to less rounds than levers, equal to, and more than.
4.  Evaluate algorithms using a variety of measures, including simple, interval, switching, and dynamic regret. Also consider other metrics such as gain that match closer to what experimenters tend to consider in realistic applications.
5.  Evaluate methods that are built for stationary rewards against those built for non-stationary rewards with data simulated for both environments. Sometimes surprising trends emerge. For instance, sometimes methods built for non-stationary rewards do best in both settings, Also for certain low levels of drift or settings with very few changepoints and short-duration experiments, simpler methods may do better than those specifically built for non-stationarity.
6.  Also simulate overdispersed data. This is often neglected in experiments. How well do the methods under consideration handle misaligned reward distribution assumptions?
7.  Experiment with a variety of feedback delays. We'd expect that non-randomized policies deteriorate quickly as delay increases.
8.  Also consider environments with seasonal trends and consider which methods can distinguish those from changepoints and/or drift.

## Methods

| Short Name | Long Name | Authors | Month and Year | Stationary | Contextual | Method Summary |
| ---------- | --------- | ------- | -------------- | ---------- | ---------- | -------------- |
| dUCB | discounted UCB | L. Kocsis and C. Szepesvari | 2006 | Drifting & Piecewise | No | Discount old responses by a factor $\gamma \in (0, 1]$, where $\gamma = 1$ is no discount. At each time point, combine data by discounting the previous data (which was previously discounted) by discount factor and summing. |
| dTS | discounted TS | Vishnu and Sheetal | July 2017 | Drifting | No | Apply discount factor to historical sufficient stats when combining those with stats from current observations. The rest is the same as the classic Bayesian Bernoulli bandit approach. Authors also implement an OTS variant. |
| dLinUCB | dynamic Linear UCB | Q. Wu, N. Iyer, and H. Wang | 2018 | Piecewise | Yes | Detect changepoints by maintaining plausible (UCB) models and evaluating each in terms of prediction error. Those beyond an error threshold are discarded and when no suitable models remain, a new one is introduced. |
| dLinTS | dynamic Linear TS (a Weighted Least Squares approach) | Burtini, Loeppky, and Lawrence | 2015 | Drifting | Yes | WLS method that weights observations by inverse of recency to discount older data in a regression model with detrending and autoregressive terms. Found that autoregressive terms weren't as important as linear deternding in a variety of drifting scenarios. Very promising method; experiments cover a variety of drifting scenarios and compare to a variety of other methods. |
| adTS | adaptive TS | Hariri, Mobasher, and Burke | 2015 | Piecewise | Yes | Changepoint detection through combo of CUSUM charts and bootstrap, using fixed-length sliding window. Change is detected based on Mahalonobis distance between model trained on current and previous window. Conjugate Bayes multivariate Gaussian  regression for estimating rewards and TS for policy. Represent users and items as features and factorize using PCA. Method looks a bit hacky -- 5+ parameters to tune and sparse mention of how that was done. |
| KF-MANB | Kalman Filtered Multi-Armed Normal Bandit | Granmo and Berg | June 2010 | Drifting | No | Estimate Gaussian rewards using Kalman Filter on all model parameters and use TS for policy. Must provide observation and transition noise parameters, but sensitivity analysis shows method is robust to settings. Very promising results compared to strong competitors. |
| POKER | Price of Knowledge and Estimated Reward | Vermorel and Mohri | 2005 | Yes | Yes | Choose arm with highest estimated price, where price is estiamted as sum of estimated mean reward and exploration bonus. This bonus is the product of: (1) the probability of an improvement, (2) expected reward mean, and (3) number of time steps remaining. This method assumes a known horizon but that seems to be an intuitive hyperparameter compared to other exploration params like $\epsilon$ for other methods. |
| SW-UCB | Sliding Window UCB | Garivier and Moulines | Apr. 2008 | Piecewise | No | Specify a fixed sliding window length $\tau$ and only estimate rewards from the last $\tau$ observations. |
| tuned SW-UCB | tuned Sliding Window UCB | Cheung, Simchi-Levi, Zhu | Oct. 2018 | Drifting | Yes | Window size is set using known variation budget $B_T$, feature dimension $d$, and time horizon $T$: $w = \lfloor d^{2/3}T^{2/3}B_T^{-2/3} \rfloor$. $B_T$ is a bound on the sum of $l_2$ differences of consecutive $\theta_t$'s, where $\theta_t$ are the effects of context variables, i.e. reward = linear combination of these effects. |
| BOB | Bandits-over-bandits | Cheung, Simchi-Levi, Zhu | Oct. 2018 | Drifting | Yes | An alternative to tuned SW-UCB for when the variation budget is unknown. Divide time horizon into some number of blocks. Restart SW-UCB for each block with a window size $w_i$ chosen by a higher-level bandit (implemented using Exp3). Exp3 chooses window size, then at the end of the block, observes feedback as total normalized reward from that block. After observing this, updates beliefs and chooses next window size. |
| LM-DSEE | Limited Memory Deterministic Sequencing of Exploration and Exploitation | Wei and Srivastava | Feb. 2018 | Piecewise & Drifting | No | Interleave blocks of exploration and exploitation. Use only rewards from current exploration epoch to estimate mean rewards. Increase length of exploitation epoch using a power law. Big weakness: 5 params, which authors state they "tuned based on environment." I'm skeptical of their results. |
| SW-UCB# | Sliding Window UCB# | Wei and Srivastava | Feb. 2018 | Piecewise & Drifting | No | Modify SW-UCB so width of sliding window varies over time. Width chosen as $\text{min}(\lambda t^\alpha, t)$ at time $t$, where $\alpha \in (0, 1]$ and $\lambda \in \mathbb{R}_{\ge 0}$ are "tuned based on environment." I'm skeptical about the results since the tuning procedure was glossed over. |
| Exp4.S, ADA-* Algorithms | Adaptive versions of Greedy, BinGreedy, ILOVETOCONBANDITS, and Exp3.S | Luo, Wei, and Langford | June 2018 | Piecewise | Yes | Extends each of the methods listed with its own "sophisticated" statistical tests for changepoint detection. Authors only conducted theoretical analyses for these methods, but their discussion provides interesting perpsective into various notions of regret (interval, switching, and dynamic). |
| AFF-{TS, OTS, UCB, $\epsilon$-greedy} Algorithms  | Augment existing techniques with Adaptive Forgetting Factors | Lu, Adams, and Kantas | Dec. 2017 | Drifting & Piecewise | No | The adaptive forgetting factor $\lambda = (\lambda_1, \ldots, \lambda_t)$ is an expanding sequence over time. Each element is computed via a single gradient descent step. The AFF was developed previously and these authors simply incorporate the idea into existing methods. AFF-OTS (Optimistic TS) and AFF-TS perform best in the authors' experiments. |
| $\gamma$-restart | $\gamma$-restart | Hartland et al. | Nov. 2006 | Piecewise | No | Run UCB-tuned (UCBT) method until changepoint is detected. Use "discounted inertia" version of Page-Hinckley test to only trigger in change-point case, rather than drifting case. Once triggered, discount old data by a factor $\gamma$ and restart on the new data with no discount. Specifically, the mean reward estimates are unchanged but the sample sizes are discounted by $\gamma \in (0, 1)$. |
| Adapt-EvE | Adapt-EvE | Hartland et al. | Nov. 2006 | Piecewise | No | Same as $\gamma$-restart until a changepoint is detected. Then a meta-bandit is initialized with two arms; one which continues using the trained version of UCBT, and the other which resets all params and instantiates a new instance of Adapt-EvE. Training then continues at both levels. Considered a **meta-bandit algorithm** because it uses a bandit algorithm to choose which lower-level bandit to use. Some other papers refer to this method as "Meta-Bandit." |
| IDS | Information-directed Sampling | Russo and Van Roy | July 2017 | Yes | Yes | The authors introduce a randomized probability matching method called IDS (Information-directed sampling) that myopically minimizes the squared regret incurred per bit of information acquired about the optimum. Specifically, IDS minimizes the ratio between squared expected single-period regret and the mutual information between the optimal action and next observation. This method outperforms UCB, TS, and the Knowledge Gradient (KG) algorithm. |
| KG | Knowledge Gradient | Ryzhov, Powell, Frazier | July 2008 | Yes | Yes | Selects actions by maximizing estimated mean reward ($\mu_{t,a}$) plus immediate improvement in decision quality due to sampling the action and observing its reward ($v_{t,a}^{KG}$): $\mu_{t, a} + (T - t)v_{t,a}^{KG}$. This places value on indirect information but may still have issues with cumulating info (according to the IDS authors Russo and Van Roy). |

## Bibliography

1.  Q. Wu, N. Iyer, and H. Wang, “Learning Contextual Bandits in a Non-stationary Environment,” The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval  - SIGIR ’18, pp. 495–504, 2018.
2.  N. Hariri, B. Mobasher, and R. D. Burke, “Adapting to User Preference Changes in Interactive Recommendation,” in IJCAI, 2015.
3.  O.-C. Granmo and S. Berg, “Solving Non-Stationary Bandit Problems by Random Sampling from Sibling Kalman Filters,” in Trends in Applied Intelligent Systems, 2010, pp. 199–208.
4.  J. Vermorel and M. Mohri, “Multi-armed Bandit Algorithms and Empirical Evaluation,” in Machine Learning: ECML 2005, 2005, pp. 437–448.
5.  L. Wei and V. Srivastava, “On Abruptly-Changing and Slowly-Varying Multiarmed Bandit Problems,” arXiv:1802.08380 [cs, stat], Feb. 2018.
6.  A. Garivier and E. Moulines, “On Upper-Confidence Bound Policies for Non-Stationary Bandit Problems,” arXiv:0805.3415 [math, stat], May 2008.
7.  W. C. Cheung, D. Simchi-Levi, and R. Zhu, “Learning to Optimize under Non-Stationarity,” arXiv:1810.03024 [cs, stat], Oct. 2018.
8.  H. Luo, C.-Y. Wei, A. Agarwal, and J. Langford, “Efficient Contextual Bandits in Non-stationary Worlds,” arXiv:1708.01799 [cs, stat], Aug. 2017.
9.  M. Speekenbrink and E. Konstantinidis, “Uncertainty and Exploration in a Restless Bandit Problem,” Top Cogn Sci, vol. 7, no. 2, pp. 351–367, Apr. 2015.
10. X. Lu, N. Adams, and N. Kantas, “On Adaptive Estimation for Dynamic Bernoulli Bandits,” arXiv:1712.03134 [cs, stat], Dec. 2017.
11. G. Burtini, J. Loeppky, and R. Lawrence, “Improving Online Marketing Experiments with Drifting Multi-armed Bandits,” in Proceedings of the 17th International Conference on Enterprise Information Systems - Volume 1, Portugal, 2015, pp. 630–636.
12. C. Hartland, S. Gelly, N. Baskiotis, O. Teytaud, and M. Sebag, “Multi-armed Bandit, Dynamic Environments and Meta-Bandits.” 14-Nov-2006.
13. M. Abeille and A. Lazaric, “Linear Thompson Sampling Revisited,” in Artificial Intelligence and Statistics, 2017, pp. 176–184.
14. O. Besbes, Y. Gur, and A. Zeevi, “Optimal Exploration-Exploitation in a Multi-Armed-Bandit Problem with Non-Stationary Rewards,” Social Science Research Network, Rochester, NY, SSRN Scholarly Paper ID 2436629, Nov. 2018.
15. A. Slivkins and E. Upfal, “Adapting to a Changing Environment: the Brownian Restless Bandits,” Jul. 2008.
16. C.-Y. Wei, Y.-T. Hong, and C.-J. Lu, “Tracking the Best Expert in Non-stationary Stochastic Environments,” in Advances in Neural Information Processing Systems 29, D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett, Eds. Curran Associates, Inc., 2016, pp. 3972–3980.
17. D. Russo and B. Van Roy, “Learning to Optimize via Information-Directed Sampling,” in Advances in Neural Information Processing Systems 27, Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, Eds. Curran Associates, Inc., 2014, pp. 1583–1591.
18. I. O. Ryzhov, W. B. Powell, and P. I. Frazier, “The Knowledge Gradient Algorithm for a General Class of Online Learning Problems,” Oper. Res., vol. 60, no. 1, pp. 180–195, Jan. 2012.
19. V. Raj and S. Kalyani, “Taming Non-stationary Bandits: A Bayesian Approach,” arXiv:1707.09727 [cs, stat], Jul. 2017.
20. K. Chaloner and I. Verdinelli, “Bayesian Experimental Design: A Review,” Statist. Sci., vol. 10, no. 3, pp. 273–304, Aug. 1995.
21. B. Dumitrascu, K. Feng, and B. E. Engelhardt, “PG-TS: Improved Thompson Sampling for Logistic Contextual Bandits,” arXiv:1805.07458 [cs, stat], May 2018.

# Learning Contextual Bandits in a Non-stationary Environment
By: Q. Wu, N. Iyer, and H. Wang
2018

**Environment: Piecewise stationary**

**Method: Multi-model UCB-based Changepoint detection (dLinUCB)**

-   General idea:
    -   Maintain multiple candidate models and estimate distribution of error for each.
    -   Prediction error for each candidate model is estimated based on last set of observations in a sliding window.
    -   The model chosen for the current policy is the one with highest LCB on its error distribution.
-   Other details:
    -   Model abandonment: If model error exceeds its UCB on error, it is discarded.
    -   Model creation: If no suitable models remain, a new one is created.
    -   Model updating: Each candidate model is updated with data from the time of creation to the time their error exceeds their UCB.
-   In this paper, the LinUCB model is used for each of the candidate models, but the method can be generalized to other models.

**Details:**

For each slave bandit model $m$, define a binary random variable $e_i(m)$. This indicates whether the slave model's prediction error at time $i$ exceeds its confidence bound:

$$e_i(m) = \mathbb{1}\{ |\hat{r}_i(m) - r_i| > B_i(m, a_i) + \epsilon \}$$

where:
-   $B_i(m, a_i)$ = upper confidence bound of model $m$ reward estimation (method-specific, see paper for LinUCB equation)
    -   note, the dependence on the arm $a_i$ assumes the arm has some set of features and users have preferences for those
-   $|\hat{r}_i(m) - r_i|$ is the error in the reward estimation, where $r_i$ is the observed reward at time $i$.
-   $\epsilon = \sqrt{2} \sigma erf^{-1}$. This represents the high probability bound of Gaussian noise in received feedback.
    -   $erf^{-1}$ is the inverse of the Gauss error function

According to theorem 3.1 in the paper, for a linear model $m$, if the underlying environment is stationary, $\mathbb{P}(e_i(m) = 1) \le \sigma_1$, where $\sigma_1 \in (0, 1)$ is a hyper-parameter in $B_i(m, a_i)$. In practice, the error is estimated as the mean prediction error on the last $\tau$ observations ($\tau$-sized sliding window). So $\hat{e}_i(m) = 1/\hat{\tau}(m) \sum_{i=t-\hat{\tau}(m)}^t e_i(m)$, where $\hat{\tau}(m) = \min{t - t_m, \tau}$, and $t_m$ is when the model was created.

**Regret Analysis: $O(\Gamma_T \sqrt{S_{max}} log S_{max})$**

-   $\Gamma_T$ = total number of ground-truth environments up to time $T$.
-   $S_{max}$ = longest stationary period till time $t$.
-   Authors state: "This arguably is the _lowest_ upper regret bound any bandit algorithm can achieve in such a non-stationary environment without further assumptions."
-   Not too impressive: many assumptions are unnecessary, and complete independence across environments seems to be one of those.

**Experiments:**

-   Compared against:
    1.  LinUCB
    2.  adaptive Thompson Sampling (adTS): changepoint detection module
    3.  Windowed mean-shift detection algorithm (WMDUCB1): UCB1-type with changepoint detection
    4.  Meta-Bandit algorithm: switches between two UCB1 models
    5.  Collaborative filtering with bandit learning (CLUB) -- on the real-world datasets
-   Datasets: 1 synthetic and 2 real-world
    -   Synthetic:
        -   Size-$K$ ($K$ = 1000) arm pool, with each arm associated to $d$-dimensional feature vector $x_a$.
        -   Also $\theta^*$, which is the global user preferences for each of the $d$ features.
        -   All parameters drawn from uniform (0, 1).
        -   All rewards generated by multiplying $x_a \theta^*$ are corrupted by Gaussian noise before returned to learner.
        -   After $S$ rounds, $\theta^*$ is randomized until a specified distance $\Delta$ between rewards returned is achieved for some specified proportion $\rho$ of the arms.
        -   All algorithms executed for 5000 iterations.
        -   Accumulated regret used to evaluate each.
        -   **Ranking:** (1) dLinUCB, (2) adTS, (3) LinUCB, (4) Meta-Bandit, (5) WMDUCB1
    -   Real-world 1: Yahoo! Today Module
        -   Clickstream dataset with 45,811,883 user visits in ten-day period in May 2009.
        -   For each visit, both user and 10 candidate articles are associated with feature vector of six dimensions.
        -   Optimizing for CTR.
        -   Replay is used to evaluate all methods, based on CTR normalized by logged random strategy's CTR.
        -   Two modes of testing:
            1.  Estimate latent user preferences for items with 2 variants: (1) non-personalized: assume all users share same preferences, (2) personalized: assume each user has own latent preferences.
                -   **Ranking:** (1) dLinUCB, (2) UCB, (3) adTS & CLUB basically tied
            2.  Estimate article popularity over time
                -   **Ranking:** (1) dLinUCB & LinUCB tied, (2) adTS
    -   Real-world 2: LastFM & Delicious
        -   HetRec 2011 workshop data
        -   LastFM contains 1892 users and 17632 items (artists) -- listened artists = positive feedback
        -   Delicious contains 1861 users and 69226 items (URLs) -- bookmarked URLs = positive feedback
        -   For both, tags preprocessed into TFIDF vectors, then PCA to reduce to 25 principle components
        -   $K$ fixed to 25, with 1 of 25 picked to be positive for a particular user
        -   Followed Hartland et al. (Multi-armed Bandit, Dynamic Environments, and Meta-Bandits) to simulate non-stationarity.
        -   **Ranking:** (1) dLinUCB, (2) adTS, (3) LinUCB

**Thoughts:**

-   Using a Bayesian approach, we could estimate these error distributions using the posterior predictive distributions.
-   Another option would be to use the likelihood on the past $\tau$ samples. Perhaps when this drops low enough, discard the model? I'm not sure what the bound would be on the likelihood that is equivalent to the bound given in theorem 3.1.
-   It's not entirely clear how much of the benefit comes from the improved changepoint detection and how much comes from maintaining multiple candidate models.
-   dLinUCB showed marked improvement over adTS on the first two datasets but nearly tied with it on the last two. The Yahoo! dataset actually seems closest to the domain I'm interested in, so this may not be concerning. It doesn't seem to be cherry-picking, but it would be interesting to try to tease out why this is the case.

# Adapting to User Preference Changes in Interactive Recommendation
By: Hariri, Mobasher, and Burke
2015

**Environment: Piecewise stationary**

**Method: TS-based Changepoint detection (adTS)**

-   General idea:
    -   Changepoint detection through combo of CUSUM charts and bootstrap (from Wayne, 2000 - change-point analysis...)
    -   Use two sliding intervals of fixed length $N$.
    -   Fit two models, one on each, and compare their distributions using Mahalanobis distance.
    -   Changepoint detection method produces confidence of changepoint. Pre-defined threshold used to trigger detection above certain confidence.
    -   Accept detected changepoint if there are at least L (look-ahead parameter) points after the change.
-   Other details:
    -   Start checking for changepoints after S (splitting threshold parameter) time points of observations.
    -   What happens when change is detected? Options are to combine data before and after in some way or just discard before. The method in this paper discards before.
    -   Use conjugate Bayesian multivariate Gaussian regression for estimating rewards and Thompson Sampling for policy.
    -   Represent both users and items as set of features. Factorize items using PCA.
    -   Learn latent user preferences and detect changes in these.

**Experiments:**

-   Dataset description:
    -   Yahoo! Music ratings dataset version 1.0. Over 10M ratings of musical artists over 1mo period prior to March 2004.
    -   Ratings of 0 to 100 given by 1,948,882 users to 98,213 artists.
-   Two goals: evaluate accuracy of changepoint detection and compared to conventional recommenders.
-   Evaluation method: 5-fold CV: one fold for testing, one for tuning, three for training.
    -   Artists not rated were given rating of overall mean rating across dataset.
-   Changes simulated by merging two users together, providing ratings from one, then switching to the other after $T = 30$ rounds.
-   Filtered test users so each had at least $T$ ratings in user's profile.
-   Compared against:
    1.  User-based kNN (k = 10)
    2.  Standard Beta-binomial TS
    3.  Optimal recommender (always serves item with highest rating, knows the changepoint exactly)
-   **Ranking:** (1) Optimal (2) adTS, (3) Standard TS, (4) User-based kNN
    -   As expected, adTS was no better than standard TS until after the changepoint.
    -   It probably would have looked much better after a few changepoints.

**Thoughts:**

-   The experiments in this paper were very limited. It would have been much more interesting to try out varying hybrid user sizes. For instance, merge 2-5 users and see how the method fares when there are more changepoints.
-   It also seems like the method of merging users into hybrid users should have ensured actual differences between users were present. Otherwise the simulated change might actually not have been a real change in preferences.
-   Overall, it seems like this method is a serious hack-job. It has _so_ many parameters.
    1.  Interval size $N$
    2.  Limits on both how many points ahead (Look-ahead $L$) and behind (splitting threshold $S$) before detecting
    3.  Number of latent features $d$
    4.  Change detection confidence (set to 0.95 for all experiments in paper)
    5.  Several model hyperparameters (though these can be set to be uninformative or weakly-informative in the usual Bayesian manner)
-   Also, the authors conduct no sensitivity analyses, so it's unclear how much their choices were cherry-picked. Relatively competitive performance in experiments by Wu, Iyer, and Wang show good performance and alleviate these concerns somewhat.
-   It is interesting to think of segments/audiences as "users." Then we can try to factorize the set of all experiences, or "items."
Grouping users into segments provides a means to use collaborative filtering. So this approach may be a feasible alternative to content tagging.

# Solving Non-Stationary Bandit Problems by Random Sampling from Sibling Kalman Filters
By: Granmo and Berg
June 2010

**Environment: Constantly Changing (Drifting) with Normally Distributed Rewards**

**Method: Kalman Filtered Conjugate Gaussian Estimation (KF-MANB)**

-   Presents TS algorithm for normal rewards that incorporated Kalman Filter on all value model parameters.
    -   Although authors don't appear to be aware that they're doing TS after estimating the distributions.
    -   They state in a footnote "To the best of our knowledge, the concept of having automata choose actions based on the _order of statistics_ of instances of estimate distributions, has been unreported in the literature."
-   Must provide observation and transition noise parameters, but sensitivity analysis shows method is robust to settings.
-   **Not contextual**; estimates single mean and variance for each arm at each time point.

**Experiments:**

-   Summary:
    -   On both stationary & non-stationary simulated data, this method outperforms all others.
    -   Next contenders are POKER and $\epsilon^n$-greedy ($\epsilon$ reduced on some schedule).
-   Compared to:
    1.  UCB1-Normal
    2.  Interval Estimation (UCB)
    3.  Pursuit (top-performing method from Learning Automata field)
    4.  $\epsilon^n$-greedy
    5.  POKER
-   Evaluation in terms of regret
-   Report results for 10-armed stationary and non-stationary simulated datasets
    -   Ensemble of 1000 independent replications with different RNGs used
    -   In each replication, 2000 arm pulls conducted
    -   Set difference in reward distribution means to be 50.0 and space them evenly
    -   Simulated several variations with increasing observation noise and fixed transition noise (stationary) for wide range of signal-to-noise ratio (12.5 to 200 grid)
        -   All other methods do much worse with higher observation noise settings
    -   Also simulated several variations with increasing transition noise and fixed observation noise (0 to 200 grid)
        -   Even the next-best method (POKER) is much worse with even a small amount of transition noise.
    -   **Ranking:** (1) KF-MANB (2) POKER, (3) $\epsilon^n$-greedy with $c=0.3$, other settings not much worse

**Thoughts:**

-   Very promising method for continuously changing environments.
-   Interesting to note how poorly the other methods which are designed for changepoint detection do in drifting environments.
-   The term "sibling" in "Sibling Kalman Filters" in the title is a bit nebulous. I wonder if the siblings are the filters for the different arms?
-   Clearly there's a gap in the literature with applying Kalman Filters to the more interesting linear models capable of working in contextual settings.

# Multi-armed Bandit Algorithms and Empirical Evaluation
By Vermorel and Mohri
2005

**Environment: Stationary and Non-Contextual**

**Summary:**

-   Evaluate many methods on two stationary datasets: (1) simulated mixed-duration with normal rewards and (2) CDN latency minimization network dataset.
-   Authors posit that main contributions are empirical valuation of many methods and introduction of theirs, which, obviously, outperforms all the others.
-   Also provides a nice summary of existing methods.

**Method: Price of Knowledge and Estimated Reward (POKER)**

-   Based on 3 ideas:
    1.  Natural way of balancing exploration and exploitatoin is to assign a price to knowledge gained while pulling a particular lever. This notion of "value of information" has been studied in other domains and is often referred to as "exploration bonuses" in bandit literature. "The objective is to quantify the uncertainty in the same units as the rewards."
    2.  Properties of unobserved arms could potentially be estimated, to some extent, from those already observed. This idea motivates the use of other arm estimates in determining the expected magnitude of improvement from pulls that yield improvements.
    3.  Specification of the expected time horizon remaining provides an intuitive way of tuning exploration vs. exploitation. The value of exploration goes down when there is less time to exploit what is learned through it.
-   POKER chooses the arm with the highest estimated price.
    -   The price is the sum of the estimated mean reward and the exploration bonus.
    -   The exploration bonus is the product of the probability of an improvement, the expected reward mean improvement, and the horizon remaining.

**Details:** sketch of method:

1.   Estimate means (normally distributed) for each arm: $\hat{\mu}_i$
2.   Sort them largest to smallest: $\hat{\mu}_{i_1} \ge \ldots \ge \hat{\mu}_{i_q}$
3.   Compute index of benchmark as $\sqrt{q}$, where $q$ is the number of arms with nonzero rewards so far.
4.   Compute estimated reward improvement as $\delta_\mu = \frac{ \hat{\mu}_{i_1} - \hat{\mu}_{i_{\sqrt{q}}} }{q}$.
5.   Compute probability of improvement as $P[\mu_i \ge \hat{\mu}_i^* + \delta_\mu]$, where $\hat{\mu}_i^*$ is the highest estimated mean reward.
6.   Compute price of each arm $i$ as: $p_i = \hat{\mu}_i + P[\mu_i \ge \hat{\mu}_i^* + \delta_\mu]\delta_\mu H$.
7.   Serve arm with highest price and observe feedback.
8.   Update model mean estimates and repeat.

**Experiments**

-   Simulated data
    -   all arms normally distributed; means and stdevs drawn from uniform (0, 1)
    -   data generated for 1000 arms and 10K rounds.
    -   3 configurations: 100, 1000, and 10K rounds, corresponding to less rounds than levers, equal to, and more than.
    -   idea of these configs is to evaluate in non-asymptotic cases as well as closer to asymptotic
-   Real-world: URLs Retrieval Latency
    -   Agent selects one source and waits until the data are retrieved (latency feedback); objective is minimization of total latency
    -   Home pages (arms) of 700+ universities, retrieved roughly every 10 min for about 10 days.
    -   Evaluated by running for 130 and 1300 rounds
-   Ranking: POKER is best in both settings. IntEstim (UCB) second best. $\epsilon$-first and $\epsilon$-greedy also do quite well. GaussMatch (TS) does well with enough rounds, but poorly when there are few rounds.

**Thoughts:**

-   Interesting tidbit in the lit review section: "A different version of the bandit problem has been studied by [10, 23, 9, 8] where the reward distributions are assumed to be known to the player. This problem is not about balancing exploration and exploitation, it admits an optimal solution based on the so-called Gittins indices."
-   Another interesting tidbit: "This paper considers the _opaque_ bandit problem where a unique reward is observed at each round, in contrast to the _transparent_ one where all rewards are observed." I hadn't heard of this distinction clarified by this opaque/transparent terminology before; it's a nice way to do it.
-   The notion of estimating properties of unobserved arms from those observed seems similar to a common-prior assumption. Another way they could be estimated is by factorizing the user and item spaces (i.e. contextual and multivariate extensions). Is there a unifying way to look at information gained that spans these problem settings?
-   **It's a good idea** to run experiments of multiple durations, as performance with less rounds can be _way_ different than performance after many rounds.

# On Abruptly-Changing and Slowly-Varying Multiarmed Bandit Problems
By Lai Wei and Vaibhav Srivastava
Feb. 2018

**Environment: Piecewise stationary and Drifting**

**Summary:**

-   Introduced two new methods, performed regret analysis of each in both environments, and conducted simulated experiments with just these methods.

**Method: LM-DSEE and SW-UCB#.**

-   **Limited Memory Deterministic Sequencing
of Exploration and Exploitation (LM-DSEE)**: interleave blocks of exploration and exploitation.
    -   Has five parameters, which are tuned based on environment.
    -   Motivated by desire to bias less towards old data in a constantly changing world, modified DSEE algorithm in two ways:
        1.  Use only the rewards from the current exploration epoch to estimate the mean rewards
        2.  Increase the length of exploitation epoch using a power law
-   **Sliding-Window Upper Confidence Bound# (SW-UCB#)**: modify SW-UCB so width of sliding window varies over time.
    -   Width is chosen as $\tau(t, \alpha) = \text{min}\{\lambda t^\alpha, t\}$ at time $t = \{1, \ldots, T\}$, where $T$ is the horizon.
    -   Parameters $\alpha \in (0, 1]$ and $\lambda \in \mathbb{R}_{\ge 0}$ are tuned based on environment.

**Interesting quotes:**

-   The MAB problem is a prototypical example of the explore-versus-exploit tradeoff: choosing between the most informative and seemingly the most rewarding alternative.
-   In a non-stationary environment, achieving logarithmic expected cumulative regret may not be feasible and the focus is on design of algorithms that achieve sublinear expected cumulative regret...
-   A sliding-window of observations is used in the SW-UCB algorithm and tuning of the fixed width of this sliding-window requires the knowledge of the the horizon length of the problem. In the SW-UCB# algorithm, we relax this requirement by considering a time-varying length of the sliding-window.
-   In an abruptly-changing environment, the mean rewards from arms switch to unknown values at unknown time-instants.
-   In a slowly-varying environment, the change in the mean reward at each arm between any two subsequent time-instants is small and is upper bounded...

**Experiments**

-   10-armed bandit simulation for both environments
-   Each method was evaluated for a few different settings of parameter values used to simulate data.
-   The point of these experiments was to empirically demonstrate the upper bounds from the analysis.

**Thoughts:**

-   Both methods have multiple parameters that are "tuned based on environment characteristics," but the authors don't mention how they tune them in their experiments.
-   In addition, neither method is compared against the other or against any others.
-   The first is a serious limitation and the second an annoyance for practitioners. These limitations make me doubt the utility of either method.

# On Upper-Confidence Bound Policies for Non-Stationary Bandit Problems
By Garivier and Moulines
Apr. 2008

**Environment: Piecewise Stationary**

**Summary:**

-   Introduce SW-UCB and analyze it and dUCB (previously introduced) in terms of regret. Conduct two simple Monte Carlo experiments to support findings.

**Method: dUCB and SW-UCB**

-   **Discounted UCB (dUCB)**
    -   Discount old responses by a discount factor $\gamma \in (0, 1]$, where $\gamma = 1$ is no discount.
    -   Assuming binary rewards:
        -   mean for arm $i$ is estimated as $\bar{X}(\gamma, i) = \frac{1}{N_t(\gamma, i)} \sum_{s=0}^t \gamma^{t-s} X_s(i) \mathbb{I}\{a_s = i\}$
        -   sample size is estimated as $N_t(\gamma, i) = \sum_{s=0}^t \gamma^{t-s} \mathbb{I}\{a_s = i\}$.
    -   Since this is UCB, there is also a discounted padding function that is added to the estimated means to get the upper bounds.
    -   At each time point, combine data by discounting the previous cumulative data by discount factor and summing.
-   **Sliding Window UCB (SW-UCB)**
    -   Specify a sliding window length $\tau$ and only estimate rewards from the last $\tau$ observations.

**Details:**

-   Both methods have cumulative regret bounds of $O(\sqrt{T log(T)})$.

**Experiments:**

-   Two different simulated experiments with $K=2$ and $K=3$ arms. Compares dUCB, SW-UCB, UCB1, and Exp3.S
-   In this first, the winner changes abruptly at some time points in the horizon $T=10^4$.
    -   dUCB and SW-UCB perform similarly, followed by Exp3.S then UCB1
-   In the second ($K=2$), the rates are constantly evolving. One is set to 0.5 and the other evolves periodically so that the winner changes periodically.
    -   Same rank ordering as first experiment

# Learning to Optimize under Non-Stationarity
By Cheung, Simchi-Levi, and Zhu
Oct. 2018

**Environment: Drifting, linear bandits**

**Summary:**

-   Introduces tuned SW-UCB method and Bandits-over-bandits (BOB) framework to alleviate need to choose window size; conducts regret analysis and numerical experiments to evaluate in terms of _dynamic regret_.

**Method: tuned SW-UCB and BOB**

-   Choose SW-UCB if we know variation budget $B_T$, else BOB.
    -   "...a drifting environment, which is a hybrid of a stochastic and an adversarial environment... can be dynamically and adversarially changed, the total change (quantified by a suitable metric) in a T step problem is upper bounded by $B_T (= \Theta(T^\rho))$ for some $\rho \in (0, 1)$, the variation budget."
-   **tuned SW-UCB**
    -   It's "tuned" in the sense that it uses $B_T$ and $T$, if known to set the window size $w$
    -   If $B_T$ is known, the learner can set $w = \lfloor d^{2/3}T^{2/3}B_T^{-2/3} \rfloor$.
    -   Else, set $w = \lfloor (dT)^{2/3} \rfloor$.
    -   $d$ is the dimension (number of predictors), $T$ is the horizon
-   **Bandits-over-bandits (BOB)**
    -   Divide time horizion into $b = \lceil T / H \rceil$ blocks
    -   Specify a set $J$ from which each $w_i$ is drawn from
    -   Restart SW-UCB with window size $w_i$ for each block $b_i$
    -   Choose window sizes using EXP3 algorithm, with total reward from each block as feedback.
        -   EXP3 chooses window size $w_i$, then at end of block $b_i$, observes feedback as total reward from that block
        -   Then updates beliefs and chooses next window size.
        -   If Exp3 chooses window size > $H$, effective size is $H$ since that's the max block size

**Details:**

-   Dynamic regret incorporates the idea that the winner may change over time and considers the optimal choice at any given point in time.
-   Define drifting environment as "the sum of $l_2$ differences of consecutive $\theta_t$'s should be bounded by some variation budget $B_T = \Theta(T^\rho)$ for some $\rho \in (0, 1)$, i.e."

$$\sum_{t=1}^T ||\theta_{t+1} - \theta_t || \le B_T$$

**Experiments: 3 with simulated data**

1.  2-arm bandit with sinusoidal reward functions and varying time horizons.
    -   Known variation budget: SW-UCB achieves regret of about 20% of that achieved by Modified EXP3.S
    -   Unknown variation budget: BOB outperforms SW-UCB by a smaller but still meaningful margin
2.  piecewise linear reward functions with 2 arms
    -   BOB performs about the same as EXP3.S, while the realistic SW-UCB with unknown variation budget is much worse
3.  piecewise linear in "linear bandit" setting -- full factorial of arms resulting in many more
    -   BOB outperforms UCB, but only by a small amount. Authors hypothesize that given the many arms, UCB happens to choose a good one.

**Interesting quotes:**

-   Traditionally, most MAB problems are studied in the stochastic (Auer et al., 2002a) and adversarial (Auer et al., 2002b) environments. In the former, the model uncertainty is static and the partial feedback is corrupted by a mean zero random noise. The learner aims at estimating the latent static environment and converging to a static optimal decision. In the latter, the model is dynamically changed by an adversary. The learner strives to hedge against the changes, and compete favorably in comparison to certain benchmark policies. While assuming a stochastic environment could be too simplistic in a changing world, sometimes the assumption of an adversarial environment could be too pessimistic.
-   Instead of assuming the stochastic environment, where reward function remains stationary across the time horizon, we allow it to change over time. Specifically, we consider the general _drifting_ environment..." It is interesting that a "stochastic" environment is considered to be different than a "drifting" environment, even though both are actually stochastic, and in fact, the drifting environment is _more_ stochastic.

**Thoughts:**

-   Problem is motivated through advertisting CTR optimization setting, which is similar to ours.
-   The section in the intro that outlines "identifiable challenges" is a nice way to lay out why the problem is interesting. This paper listed (1) Uncertainty in reward distributions, (2) Non-stationarity, and (3) Partial/bandit feedback.

# Efficient Contextual Bandits in Non-stationary Worlds
By Luo, Wei, Agarwal, and Langford
June 2018

**Environment: Piecewise Stationary, Contextual**

**Summary:**

-   Equip existing methods with changepoint detection methods. Evaluate in terms of interval regret, switching regret, and dynamic regret -- useful notions for non-stationary environments. No experiments, just theoretical analysis.

**Method: Exp4.S, ADA-Greedy, ADA-ILTCB, ADA-BinGreedy**

-   First, extend Exp3.S to contextual setting, calling this Exp4.S
-   Extends the following methods, each with its own custom non-stationarity test:
    -   $\epsilon$-greedy, called ADA-Greedy
    -   ILOVETOCONBANDITS, called ADA-ILTCB (ILTCB = I Love To Con Bandits)
-   Finally, develop ADA-BinGreedy, which is a parameter-free version of ADA-Greedy that achieves good dynamic regret. The differences are:
    1.  Each block is further divided into bins with equal length
    2.  In addition to the small probability of exploration $\mu_t$ at each round, some bins are randomly selected for pure exploration
    3.  The non-stationarity test is only executed in exploration bins, and only checks for intervals within the bin
    4.  Parameters L and v are removed and the exploration probability $\mu_t$ is set adaptively.

**Details:**

-   Discussion of various notions of regret:
    -   "A high-level outcome of our analysis is that the algorithms enjoy a regret bound on any time interval that is sufficiently stationary (called interval regret), compared with the best fixed policy for that interval."
    -   "For example, if the data-generating process is typically i.i.d., except there are hard switches in the data distribution every so often, then our algorithms perform as if they knew the change points in advance, up to a small penalty in regret (called switching regret)."
    -   "More generally, if the data distribution is slowly drifting, we can still provide meaningful regret bounds (called dynamic regret) when competing to the best policy at each time (instead of a fixed policy across all rounds)."

**Interesting quotes:**

-   The simplest oracle-efficient contextual bandit algorithm is the EPOCH-GREEDY method (Langford and Zhang, 2008) which assumes i.i.d. data.

**Thoughts:**

-   Are all of these regret calculations premised on immediate feedback? What happens to these methods in the presence of delayed feedback?
-   What are simple Bayesian changepoint detection methods that could be incorporated into TS-based methods?

# Uncertainty and Exploration in a Restless Bandit Problem
By Speekenbrink and Konstantinidis
Dec. 2014

**Environment: Drifting**

**Summary:**

-   Investigate whether humans performing a restless multi-armed bandit task use a "probability of maximum reward strategy" (e.g. TS) to make their decisions.

**Method:**

-   This isn't really a bandit study, but more of a psychological survey. The authors wanted to see how people balance explore-exploit tradeoffs in games. Turns out the most similar model to humans is Bayesian updating with Kalman Filtering.
-   Considered three mechanisms by which expected values were updated:
    1.  Bayesian updating with Kalman filter
    2.  Delta rule: popular model-free alternative to Bayesian updating (Gluck & Bower, 1988; Yechiam & Busemeyer, 2005)
        -   Learning rate is fixed in delta rule, while the "learning rate" (Kalman gain) in Bayesian updating depends on the current level of uncertainty
    3.  Decay rule: discount old data, progressively moving towards no data over time if nothing new is observed.
-   Considered several different policies (they call them choice rules):
    1.  $\epsilon$-greedy
    2.  Softmax
    3.  Softmax with an exploration bonus
    4.  Thompson Sampling (they call it: probability of max utility) only used for Bayesian updating

**Details:**

-   Considers normally distributed rewards that drift over time.
-   Formulates an ideal learner using Kalman Filters

**Interesting quotes:**

-   "In a restless bandit task, the longer an arm has not been played, the higher the (subjective) probability that it may now provide better rewards. If exploration is based on this probability, the probability of exploration increases with the time that an arm has not been played." **So restless = drifting**
-   "While Daw et al. (2006) did not find evidence that exploration is related to uncertainty, they only considered a heuristic decision strategy with an exploration bonus which increased linearly with the standard deviation of the prior distribution."

# On Adaptive Estimation for Dynamic Bernoulli Bandits
By Lu, Adams, and Kantas
Dec. 2017

**Environment: Drifting and Piecewise Stationary, Non-contextual**

**Summary:**

-   Extend $\epsilon$-greedy, UCB, and TS using "adaptive forgetting factors" and evaluate them empirically for both drifting and piecewise stationary environments.

**Method: AFF-TS, AFF-OTS, AFF-UCB, AFF-$\epsilon$-greedy**

-   The adaptive forgetting factor $\mathbf{\lambda} = (\lambda_1, \ldots, \lambda_t)$ is a expanding sequence over time, and
the forgetting factor $\lambda_t$ is computed via a single gradient descent step (see paper for details).
    -   For more on adaptive forgetting factors (AFFs), see: "Continuous monitoring for changepoints in data streams using adaptive estimation"
-   **The sections on extending TS and Optimistic TS (OTS) show how to use AFFs for updating sufficient stats.** Probably useful.

**Details:**

-   "In this work, we look at the problem of dynamic bandits where the expectation of the reward distribution changes over time, focusing on the Bernoulli reward distribution because of its wide relevance in real applications... emphasise on cases where
the changes of the reward distribution can really have an effect on the decision making."

**Interesting quotes:**

-   "Restless bandits are regarded as intractable, i.e., it is not possible to derive an optimal strategy even if the transitions are deterministic (Papadimitriou and Tsitsiklis, 1999)."

**Experiments (on simulated data)**

1.  Piecewise stationary (abruptly changing)
    -   AFF-UCB and AFF-TS perform similarly, with AFF-$\epsilon$-greedy lagging slightly behind
    -   The OTS variation performed best by a small margin
2.  Drifting
    -   Similar to piecewise stationary
3.  Large number of arms (repeat first 2 experiments with 50-100 arms)
    -   TS-based methods do better than UCB-based methods here
4.  Also conducted sensitivity analysis to learning rate -- this is only parameter, and is used in gradient descent update in AFF
    -   Conclude methods aren't too sensitive to it for Bernoulli rewards

# Improving Online Marketing Experiments with Drifting Multi-Armed Bandits
By Burtini, Loeppky, and Lawrence
2015

**Environment: Drifting, Contextual**

**Summary:**

-   Introduces weighted least squares (WLS) method that weights observations by inverse of recency to discount older data in a regression model with detrending and autoregressive terms. Found that autoregressive terms weren't as important as linear detrending in a variety of simulated drifting scenarios.

**Method: Weighted Least Squares (WLS) -- author calls it LinTS**

-   Fit a model $Y_{t,i} = \alpha_t + AR_i(p) + Trend_i(t) + A_{t,i} + \epsilon_{t,i}$
    -   $\alpha_t$ is the intercept term
    -   $Trend(t)$ is a function representing expected time trend
    -   $AR(p)$ is the autoregressive term of order $p$
    -   A_{t,i} is binary dummy variables indicating which arm is selected
    -   $\epsilon_{t,i}$ is the standard Gaussian noise term
    -   All these terms are combined into the design matrix X
-   This model is then used in two algorithms:
    1.  Weighted least squares: weight observations by inverse of recency; discount confidence of old data
        -   experiment with various weighting schemes: logarithmic, polynomial, linear, exponential, sinusoidal
    2.  Optimistic Thompson Sampling
        -   partially Bayesian treatment. Draw only values above the mean of Normal with:
        -   mean $\sum_i (\hat{\beta}_i X_{i,t})$
        -   variance $\sum_i (\hat{\text{Var}}(\hat{\beta}_i) X_{i,t}^2)$

**Details:**

-   The OTS variant seems to be inspired by Minka's 2001 thesis on approximate Bayesian inference via, e.g. Expectation Propogation

**Interesting quotes:**

-   In time-series analysis, stochastic drift is used to refer to two broad classes of non-stationarity in the population parameter being estimated: (1) cyclical or model-able drift that arise because of model misspecification and (2) the random component.
-   Often it is possible to _detrend_ nonstationary data by fitting a model that includes time as a parameter. Where the function of time is well-formed and appropriate for statistical modelling, a _trend stationary_ model can be found with this detrending process.
-   For models where detrending is not sufficient to make a process stationary, _difference stationary_ models may fit, where the differences between values in time $Y_t$ and $Y_{t-n}$ can be represented as a well-formed function appropriate for statistical modelling. Difference stationary models are represented with autoregressive models.
-   If these two detrending strategies are not sufficient to make a given process stationary, more complex filters such as a band-pass or Hodrick-Prescott filter may be applied.
-   A number of results have shown improvements by performing optimistic Thompson sampling (Chapelle and Li, 2011; May et al., 2012) where one only considers the positive uncertainty surrounding an arm estimate. Unlike UCB-based policies, traditional Thompson sampling both increases (if the draw is above the point estimate of the mean) and decreases (if the draw is below the point estimate of the mean) a prediction, depending on the sample draw; for the purpose of maximizing reward (minimizing regret), the decrease appears to have no benefit.

**Experiments:**

-   Simulated a variety of "true worlds," including "arm distribution type and parameters, arm count, and drift type from a set of functional forms including random walk, exponential random walk, logarithmic, linear (in varying degree), exponential and periodic drift (sinusoidal over varying periods)."
-   Showed results for four of these: (1) no drift, (2) linear drift, (3) logarithmic drift, and (4) random walk
-   Compared SW-UCB, UCB-tuned, $\epsilon$-greedy, UCB, and discounted UCB
-   UCB-tuned and SW-UCB were the best of the competitor techniques and outperformed their OTS technique across the board
-   "Across all true worlds, we find in general that a detrending term congruent with the _true drift_ form... outperforms all other strategies in the long run, producing a zero-regret strategy..."
-   "Surprisingly, we find that linear detrending is an effective technique for handling the random walk..."
-   Found linear WLS technique most robust across experiments, even in no-drift scenario

**Thoughts:**

-   **This paper is highly relevent!**
-   Also, the style and conciseness are something to aspire to. Very clear and easy to read.

# Multi-armed Bandit, Dynamic Environments and Meta-Bandits
By Hartland, Gelly, Baskiotis, Teytaud, and Sebag
Nov. 2006

**Environment: Piecewise Stationary**

**Summary:**

-   Introduced Adapt-EvE and $\gamma$-restart techniques, utilizing Page-Hinckley changepoint detection method.

**Method: Adapt-EvE and $\gamma$-restart**

-   Both methods start out the same: Run UCB-tuned (UCBT) method until changepoint is detected. Use "discounted inertia" version of Page-Hinckley test to only trigger in change-point case, rather than drifting case.
-   **Adapt-EvE**: Once triggered, a meta-bandit is initialized with two arms: one, which continues using the trained version of UCB-Tuned, and the other which resets all parameters and instantiates a new instance of Adapt-EvE. Training then continues at both levels.
    -   Considered a **meta-bandit algorithm** because it uses a bandit algorithm to choose which lower-level bandit to use. Some other papers refer to this method as "Meta-Bandit."
-   **$\gamma$-restart**: Once triggered, discount the old data by a factor $\gamma$ and restart on the new data with no discount. Specifically, the mean reward estimates are unchanged but the sample sizes are discounted by $\gamma \in (0, 1)$.

**Experiments:**

-   Both methods compared favorably to UCBT and discounted UCBT on simulated data experiments.

# Linear Thompson Sampling Revisited
By Abeille and Lazaric
2017

**Summary:**
Purely theoretical analysis. Just listing interesting contributions for now:

1.  "...show that the TS does not need to sample from an actual Bayesian posterior distribution and that any distribution satisfying suitable concentration and anti-concentration properties guarantees a small regret. In particular, we show that the distribution should _over-sample_ w.r.t. the standard least-squares confidence ellipsoid by a factor $\sqrt{d}$ to guarantee a constant probability of being optimistic."
2.  Showed some regret bounds for TS and "Finally, we show how our proof can be easily adapted to regularized linear optimization (with arbitrary penalty) and to the generalized linear model (GLM), for which we derive the first frequentist regret bound for TS..."

# Chasing Demand: Learning and Earning in a Changing Environment
By Keskin and Zeevi
June 2016

**Environment: Drifting slowly and rapidly and piecewise stationary, Price setting algorithms**

**Summary:**

-   Analyze regret bounds for a variety of environments for the price-setting problem. Looks at similar sorts of approaches as seen in the bandits literature -- moving windows, gradually decaying old data, and simultaneous changepoint detection and parameter estimation. All approaches are premised on a weighted least squares estimator.

**Method:**

-   "Our study presents three families of dynamic pricing policies designed to perform well in changing demand environments. The moving window and decaying weights policies in Section 3 are based on a weighted least squares estimator that discounts older observations at a certain rate. The detection policy in Section 4 uses the same weighted least squares estimator but can reduce the weight of all past observations to zero upon detecting a change. All of these policies have near-optimal performance in their respective settings, but at the same time, they use quite distinct rules for weighing past observations, which suggests that successful pricing policies in presence of smooth and bursty changes can have very different structures."
-   "unlike Besbes et al., we study the case where the variation budget B is unknown to the seller, which allows us to quantify the additional cost of adapting to an unknown variation budget"

**Thoughts:**

-   Interesting that there is another closely related field of literature (price-setting problems). Might be good to dig into this more later.
-   These authors use very similar techniques as in other papers and also consider the notation of **variation budget**. I think it would be good to dig into this idea further to make sure you understand it.

# Optimal Exploration-Exploitation in a Multi-Armed-Bandit Problem with Non-stationary Rewards
By Besbes, Gur, Zeevi
May 2014

**Environment: Drifting, non-contextual**

**Summary:**

-   Main contribution: ..."fully characterizing the (regret) complexity of a broad class of MAB problems with non-stationary reward structure by establishing a direct link between the extent of reward "variation" and the minimal achievable worst case regret."

**Method: Restart procedure for other methods**

-   Proposes a simple restart procedure that uses other methods as a subroutine.
-   Restart method every epoch, with epoch length defined as an input that remains fixed throughout execution.

**Details:**

-   "The main constraint that we impose on the evolution of the mean rewards is that their variation over the relevant time horizon is bounded by a variation budget $V_T$. The idea of **Variation Budget** was first introduced in a paper by the same authors in 2013 ("Non-stationary Stochastic Optimization"). This is an upper bound on the scale of **variation in reward functions** $f_1, \ldots, f_T$ over the time horizon $T$. This variation is defined as:

$$Var(f_1, \ldots, f_T) = \sum_{t=2}^T ||f_t - f_{t-1}||,$$

where for any bounded functions $g$ and $h$ mapping actions $\mathbb{X}$ to rewards $r \in \mathbb{R}$, we denote $||g - h|| = \text{sup}_{x \in \mathbb{X}} |g(x) - h(x)|$.

**Interesting quotes:**

-   In particular, our results highlight the tension that exists between the need to "remember" and "forget..." In a nutshell, the fewer past observations one retains, the larger the stochastic error associated with one's estimates of the mean rewards, while at the same time using more past observations increases the risk of these being biased.

# Adapting to a Changing Environment: the Brownian Restless Bandits
By Slivkins and Upfal (Microsoft Research)
July 2007

**Environment: Drifting, non-contextual**

**Summary:**

-   The main contribution in this paper is a theoretical analysis relating minimal dynamic regret to volatility (variance of Brownian motion).
-   "Consider the state-oblivious dynamic MAB problem with $k$ arms such that each arm $i$ has volatility at most $\sigma_i$. Then there exists a MAB algorithm whose steady-state regret is $O(k \sigma_{av})$, where $\sigma_{av}^2 = \frac{1}{k} \sum_{i=1}^k \sigma_i^2$."

**Method: Theoretical analysis**

**Details:**

-   "We are interested in systems that exhibit a stationary, steady-state behavior. For this reason instead of the usual Brownian motion on a real line (which diverges to infinity) we consider a Brownian motion on an interval with reflecting bounds... Our goal here is to characterize the long-term average cost of adapting to such changing environment in terms of the stochastic rate of change – the volatility of Brownian motion."
-   "The paradigmatic setting for us is one in which each arm’s state has the same stationary distribution and, therefore, all arms are essentially equivalent in the long term. In such setting the standard benchmark – the best time-invariant policy – is uninformative. Instead, we optimize with respect to a more demanding (and also more natural) benchmark – a policy that at each step plays an arm with the currently maximal expected reward."
-   The authors premise their formulation on the idea that each arm has an underlying state this is either observed or not. The state changes dynamically, and the rewards depend on the state. For our purposes, we never observe a state, so we can just think of the reward function as changing gradually over time.

**Interesting quotes:**

-   A number of models have been proposed for capturing the dynamic aspect of the MAB problem. Motivated by task scheduling, Gittins [13] considered the case where only the state of the active arm (the arm currently being played) can change in a given step, giving an optimal policy for the Baysean formulation with time discounting. This seminal result gave rise to a rich line of work (e.g. [11, 12, 32, 31, 30, 6, 29]), a proper review of which is beyond the scope of this paper.
-   In particular, **Whittle [33] introduced an extension termed restless bandits** [33, 7, 24], where the states of all arms can change in each step according to a known (but arbitrary) stochastic transition function. Restless bandits are notoriously intractable: e.g. even with deterministic transitions the problem of computing an (approximately) optimal strategy is PSPACE-hard [26]. [emphasis mine]

# Tracking the Best Expert in Non-stationary Stochastic Environments
By Wei, Hong, and Lu
2016 (NIPS)

**Environment: Drifting, non-contextual**

**Summary:**

-   The authors relate various regret results from adversarial and stochastic bandit literature, providing a characterization of degrees of non-stationarity using a set of existing and newly introduced parameters.

**Method: Theoretical analysis**

**Details:**

-   "We base ourselves in the stochastic world with non-stationary distributions, characterized by the parameters $\Gamma$ and V. In addition, we introduce a new parameter $\Lambda$..."
    -   The horizon can be divided into $\Gamma$ stationary intervals, such that $\Gamma - 1$ is the number of changepoints.
    -   V is the sum of differences between the means of consecutive distributions.
    -   $\Lambda$ is the total statistical variance of the loss distributions over the time horizon.

# Learning to Optimize Via Information-Directed Sampling
By Russo and Van Roy
July 2017

**Environment: Stationary, contextual and non-contextual**

**Summary:**

-   The authors introduce a randomized probability matching method called IDS (Information-directed sampling) that myopically minimizes the squared regret incurred per bit of information acquired about the optimum. Specifically, IDS minimizes the ratio between squared expected single-period regret and the mutual information between the optimal action and next observation. This method outperforms UCB, TS, and the Knowledge Gradient (GD) algorithm.

**Method: IDS (Information-directed sampling)**

-   See formulation below.

**Details:**
Introduce three kinds of information alternative methods fail to account for:

1.  **Indirect**: selecting an action provides useful info about other actions even if it doesn't provide useful feedback about the selected action.
2.  **Cumulating**: feedback does not immediately enable higher expected reward but can eventually do so when combined with info from subsequent actions.
3.  **Irrelevant**: some info does not help determine which actions should be selected. IDS avoids investing in acquiring this info.

## Basic Problem Formulation

-   actions ($A_t)_{t \in \mathbb{N}}$ are selected from finite action set $\mathcal{A}$
-   outcomes from selection action $A_t$ at time $t \in \mathbb{N}$ are $(Y_{t, A_t})_{t \in \mathbb{N}}$, with $Y_{t, A_t} \in \mathcal{Y}$
-   $Y_t = (Y_{t, a})_{a \in \mathcal{A}}$ is the vector of outcomes for all actions at time $t$
-   conditioned on RV $\theta$, $Y_t$ for all $t$ is an iid sequence
-   agent associates reward $R(y)$ with each outcome $y \in \mathcal{Y}$, where reward function is fixed and known
-   $R_{t,a}$ is shorthand for the reward of action $a$ at time $t: R(Y_{t, a})$
-   uncertainty about $\theta$ induces uncertainty about the true optimal action, denoted $A^* = \text{argmax}_{a \in \mathcal{A}} \mathbb{E}[R_{t,a} | \theta]$

The **$T$-period regret** of the sequence of actions $A_1, \ldots, A_T$ is: $$\text{Regret}(T) = \sum_{t=1}^T (R_{t, A^*} - R_{t, A_t})$$

-   Since we don't know $\theta$, we measure the expected regret: $\mathbb{E}[\text{Regret}(T)]$
    -   the expectation is taken over the randomness in actions $A_t$, outcomes $Y_t$, and the prior distribution over $\theta$
    -   this is commonly called _Bayesian Regret_ or _Bayes Risk_
-   action $A_t$ is based on history of observations $\mathcal{F}_t = ((A_1, Y_{1, A_1}), \ldots, (A_t, Y_{t, A_t}))$ up to time $t$

**A _randomized policy_** $\pi = (\pi_t)_{t \in \mathbb{N}}$ is a sequence of deterministic functions; $\pi_t(\mathcal{F}_t)$ specififes a probability distribution over the action set $\mathcal{A}$

-   $\mathcal{D}(\mathcal{A})$ denotes the set of probability distributions over $\mathcal{A}$
-   notation is abused to write $\pi_t(\mathcal{F}_t)$ as $\pi_t$, allowing us to say $\pi_t(a) = P(A_t = a | \mathcal{F}_t)$
    -   this is the probability under our policy of selecting action $a$ given the observed history
-   the expected regret is often written as $\mathbb{E}[\text{Regret}(T, \pi)]$ to denote the expected value when the actions are selected according to $\pi$
-   $\alpha_t(a) = P(A^* = a | \mathcal{F}_t)$ is the posterior distribution of $A^*$
    -   Note: this is the TS probability of selecting an action

### Information Gain $g_t(a)$

The **Kullback-Leibler divergence** for two probability measures $P$ and $Q$ over a common measurable space, if $P$ is absolutely continuous with respect to $Q$: $$D_{KL}(P || Q) = \int log(\frac{dP}{dQ})dP$$

The **Shannon Entropy** of probability distribution $P$ over a finite set $\mathcal{X}$ is defined as: $$H(P) = -\sum_{x \in \mathcal{X}} P(x) log(P(x))$$

The **Mutual Information** under the posterior distribution between two RVs is denoted by: $$I_t(X_1; X_2) = D_{KL}(\ P((X_1, X_2) | \mathcal{F}_t)\ || \ P(X_1 | \mathcal{F}_t) P(X_2 | \mathcal{F}_t)\ ),$$

or in other words, the KL divergence between the joint posterior and the product of the marginals.

The **Information Gain** from an action $a$ is defined as $g_t(a) = I_t(A^*; Y_{t, a})$. This is equal to the expected reduction in entropy of the posterior distribution of $A^*$ due to observing $Y_{t,a}$: $$g_t(a) = \mathbb{E}[H(\alpha_t) - H(\alpha_{t+1})\ | \ \mathcal{F}_t, A_t = a]$$

To clarify this a bit, let's write out $H(\alpha_t)$:
\begin{align}
H(\alpha_t) &= -\sum_{a \in \mathcal{A}} \alpha_t(a) log(\alpha_t(a))  \\
            &= -\sum_{a \in \mathcal{A}} P(A^* = a | \mathcal{F}_t) log P(A^* = a | \mathcal{F}_t)
\end{align}

So this info gain calculation involves computing the TS probabilities for each action for the:

1.  current sequence of observations, and
2.  the sequence we may end up with after taking action $a$ (marginalizing over the uncertainty in $Y_{t,a}$)

### Instantaneous Regret $\Delta_t(a)$

The **Instantaneous Regret** of action $a$ at time $t$ is $$\Delta_t(a) = \mathbb{E}[R_{t, A^*} - R_{t,a} | \mathcal{F}_t]$$

### Gain and Regret of Policies

Overloaded notation is used for denoting the gain and regret of a policy $\pi \in \mathcal{D}(\mathcal{A})$:

\begin{align}
g_t(\pi) &= \sum_{a \in \mathcal{A}} \pi(a) g_t(a) = \pi^T \vec{\Delta}  \\
\Delta_t(\pi) &= \sum_{a \in \mathcal{A}} \pi(a) \Delta_t(a) = \pi^T \vec{g}
\end{align}

$g_t(\pi)$ is the expected info gain when actions are selected according to policy $\pi$ and $\Delta_t(\pi)$ is defined analogously

### Information-directed Sampling (IDS)

Minimizes squared regret per bit of information acquired about optimum over all action sampling distributions:

$$
\pi_t^{IDS} = \text{argmin}_{\pi \in \mathcal{D}(\mathcal{A})} \left\{ \Psi_t(\pi) = \frac{\Delta_t(\pi)^2}{g_t(\pi)} \right\}
$$

# The knowledge gradient algorithm for a general class of online learning problems
By Ryzhov, Powell, and Frazier
July 2008

**Environment: Stationary, contextual and non-contextual**

**Summary:**

-   Selects actions by maximizing estimated mean reward ($\mu_{t,a}$) plus immediate improvement in decision quality due to sampling the action and observing its reward ($v_{t,a}^{KG}$): $\mu_{t, a} + (T - t)v_{t,a}^{KG}$. This places value on indirect information but may still have issues with cumulating info (according to the IDS authors Russo and Van Roy).

**Method: KG (Knowledge Gradient)**

-   $v_{t,a}^{KG} = \mathbb{E}[V_{t+1} | \mathcal{F}_t, A_t = a] - V_t$
-   $V_t = \max_{a'} \mu_{t, a'}$

**Details:**

-   See paper for detailed mathematics.

**Interesting quotes:**

-   We allow several variations of this basic setup: the rewards may be discounted over time, the time horizon may be finite or infinite, and our beliefs about the unknown rewards may be correlated. Correlated beliefs are not handled by the traditional bandit literature, but are significant in practice... Correlated problems, where our prior belief about the mean value of one alternative is correlated with our beliefs about mean values of other alternatives, also arise frequently in applications... A major problem class is subset selection, where we have to investigate a treatment (of one or more drugs), a portfolio (of stocks or different energy saving technologies for a building) or a set of features to sell a product. Choices are correlated because they may have common elements.

# Taming Non-stationary Bandits: A Bayesian Approach
By Vishnu and Sheetal
July 2017

**Environment: Drifting, non-contextual**

**Summary:**

-   The authors apply a discount factor to historical data for a Bernoulli bandits problem.

**Method: dTS (discounted TS)**

-   Basically implemented by applying discount factor to historical sufficient stats when combining those with stats from current observations. This is a nice general-purpose approach which can be used for any Bayesian method where inference is based on sufficient stats.
-   The authors also implement an optimistic version in the tradition of OTS.

# Bayesian Experimental Design: A Review
By Chaloner and Verdinelli
1995

This isn't actually a study of bandits, so the format differs here.

**Summary:**

-   Goal: arrive at credible interval bands that are sufficiently narrow to allow us to detect effects of interventions with confidence. A Bayesian decision-theoretic approach is used. The utility function assigns some value to the case where we detect gain when gain is present and negative value to not detecting, while also incorporating a cost for each sample used.

**Interesting quotes:**

-   The basic idea in experimental design is that statistical inference about the quantities of interest can be improved by appropriately selecting the values of the control variables. These values should be chosen to achieve small variability for the estimator of interest

**Formulation:**

-   a design $\eta$ must be chosen from set $\mathcal{H}$
-   data $y$ from a sample space $\mathcal{Y}$ will be observed
-   based on $y$ a decision $d$ will be chosen from set $\mathcal{D}$, composed of two parts:
    -   selection of $\eta$
    -   choice of terminal decision $d$
-   the unknown parameters are $\theta$ from parameter space $\Theta$

A general utility function is of the form $U(d, \theta, \eta, y)$. For any design $\eta$, the expected utility of the best decision is given by:

$$U(\eta) = \int_{\mathcal{Y}} \max_{d \in \mathcal{D}} \int_{\Theta} U(d, \theta, \eta, y) p(\theta | y, \eta) p(y | \eta) d\theta d y$$

The Bayesian solution to the experimental design problem is provided by the design with max utility:

$$U(\eta^*) = \max_{\eta \in \mathcal{H}} U(\eta)$$

# PG-TS: Improved Thompson Sampling for Logistic Contextual Bandits
By Dumitrascu, Feng, and Engelhardt
May 2018

This isn't so much a new method as a way of performing more efficient inference for Logistic Regression based TS.

**Summary:**

-   Use polya-gamma data augmentation to derive efficient Gibbs sampler for a Bayesian logistic regression model. Authors compare to the commonly used Laplacian approximation method used for Logistic Regression TS to demonstrate superior performance.