決定株のラデマッハ複雑度をモンテカルロ法により数値的に計算してみる

　MLPシリーズの金森敬文『統計的学習理論』（講談社，2015）を読んでいて，ラデマッハ複雑度のイメージが掴みにくかったので，決定株という簡単なモデルとモンテカルロ法でラデマッハ複雑度を数値計算することで，ラデマッハ複雑度のイメージを掴むことを目指します．

# 1 経験ラデマッハ複雑度

## 1.1 経験ラデマッハ複雑度

　**経験ラデマッハ複雑度**の定義は，入力空間を$\mathcal{X}$，$\mathcal{X}$上の実数値関数からなる集合を$\mathcal{G} \subset \\{f:\mathcal{X}\rightarrow \mathbb{R} \\}$として，

> 入力点の集合を$S = \\{x_1,...,x_n\\} \subset \mathcal{X}$とします．また，$+1$と$-1$を等確率でとる独立な確率変数を$\sigma_1,...,\sigma_n$とします．このとき，$\mathcal{G}$の経験ラデマッハ複雑度は，
>```math
>\hat{\mathfrak{R}}_{S}(\mathcal{G}) = \mathbb{E}_{\sigma}\left[ \sup_{g \in \mathcal{G}} \frac{1}{n} \sum_{i=1}^{n} \sigma_i g(x_i) \right]
>```
>と定義されます．ここで，$\mathbb{E}_{\sigma}[\cdot]$は，$\sigma_1,...\sigma_n$に関する期待値を意味します．

です[1]．それほど難しい定義ではないですが，イメージを掴みにくいので，もう少しかみ砕いてみます．

　経験ラデマッハ複雑度は，$\mathcal{X}$上の実数値関数からなる集合に対して定義されるものですが，以下では簡単な二値分類問題について考えてみます．$g(x)$を$+1$か$-1$を返す関数とすると，$g(x_i) = \sigma_i$ならば$\sigma_i g(x_i) = 1$で，$g(x_i) \neq \sigma_i$ならば$\sigma_i g(x_i) = -1$となるので，最も精度よく$\sigma_1,...,\sigma_n$を予測できる$g(x)$によって，$\frac{1}{n} \sum_{i=1}^{n} \sigma_i g(x_i)$の上界は達成されます．この例では，$\frac{1}{n} \sum_{i=1}^{n} \sigma_i g(x_i)$は（予測の成功数ー予測の失敗数）／（全体の数）になります．

　経験ラデマッハ複雑度は，この$\sup_{g \in \mathcal{G}} \frac{1}{n} \sum_{i=1}^{n} \sigma_i g(x_i)$をラベル$\sigma_i$の取り方について期待値をとったものです．つまり，二値分類問題では，経験ラデマッハ複雑度は，適当にふったラベルに対して，どれだけフィッティングできるかの指標です．考えているモデルが複雑になり表現力があがればあがるほど，適当にふったラベルに対してもフィッティングできるだろうと考えられるので，経験ラデマッハ複雑度も大きくなると推測できます．

（回帰問題のように，$g(x)$が実数値を返す関数とすると，$\sigma_i = +1$のとき，$g(x)$が出来る限り大きい値をとって，$\sigma_i = -1$のとき，$g(x)$が出来る限り小さい値をとるときに，$\frac{1}{n} \sum_{i=1}^{n} \sigma_i g(x_i)$の上界が達成されると考えられます．つまり，経験ラデマッハ複雑度は，適当にふったラベルに対して，どれだけ分離できるかの指標になると思います．）


## 1.2 決定株を用いた実験

　と言っても，まだイメージを掴みにくいと思うので，簡単なモデルである決定株を用いて，ラデマッハ複雑度をモンテカルロ法により算出してみます．決定株とは，深さ1の決定木のことで，1つの次元を選んで閾値より大きいか小さいかで分類するモデルです．

　計算の流れは以下です．

1. 適当な分布に従う入力のサンプル$S = \\{x_1,...,x_n\\} \subset \mathcal{X}$を定める．
2. 以下の処理を$T$回繰り返す．
    2.1 $\\{x_1,...,x_n\\}$に対して，$+1$か$-1$をとる適当なラベル$\sigma_1,...,\sigma_n$を定める．
    2.2 $(x_1, \sigma_1),...,(x_n, \sigma_n)$を用いて，最も$\sigma_i$を予測できる決定株$\hat{g}(x)$を作る．
    2.3 $\sup_{g \in \mathcal{G}} \frac{1}{n} \sum_{i=1}^{n} \sigma_i g(x_i) = \frac{1}{n} \sum_{i=1}^{n} \sigma_i \hat{g}(x_i)$を計算する．
3. $\frac{1}{n} \sum_{i=1}^{n} \sigma_i \hat{g}(x_i)$の$T$回の平均をとる．

$\frac{1}{n} \sum_{i=1}^{n} \sigma_i \hat{g}(x_i)$の$T$回の平均は，経験ラデマッハ複雑度をよく近似しているだろうと考えられます．

　以下では，入力が従う分布を$d$次元標準正規分布としました．

　1つ計算例を示します．$n=100, d=2$とします．標準正規分布から乱数をとって，適当にラベルをふったものが次の図です．

![2label.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670764/9649cbe8-7784-a29e-0b5a-ae9e71a5dc6b.png)

　このときの最適な決定株の境界は次のようになります．

![decision_stump.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670764/03e1d5ee-7bf0-1e95-6a67-e19451475820.png)


　標準正規分布の乱数に適当にラベルをふっているので全くうまく行っているようにみえませんが，一応これが最適だと思います．このときの$\sup_{g \in \mathcal{G}} \frac{1}{n} \sum_{i=1}^{n} \sigma_i g(x_i) = \frac{1}{n} \sum_{i=1}^{n} \sigma_i \hat{g}(x_i)$は，$0.22$で、ラベルをふりなおして計算を繰り返し，平均をとった値は，$0.203$でした．この値は経験ラデマッハ複雑度とほぼ一致しているだろうと考えられます．
　走り書きですが，経験ラデマッハ複雑度を計算するコードを下に示します．

```python
# ある入力サンプルとラベルに対して，最適な決定株を学習し，スコアを計算する関数
def score_decision_stump(x, sigma):
    n, d = x.shape
    divide_idx_all = np.zeros((d, 2))
    max_cumsum_all = np.zeros((d, 2))
    for i in range(d):
        sigma_sorted_by_i = sigma[x[:, i].argsort()]
        cumsum_sigma_sorted_by_i = np.cumsum(sigma_sorted_by_i)
        divide_idx_all[i, 0] = np.argmax(cumsum_sigma_sorted_by_i)
        divide_idx_all[i, 1] = np.argmin(cumsum_sigma_sorted_by_i)
        max_cumsum_all[i, 0] = cumsum_sigma_sorted_by_i[int(divide_idx_all[i, 0])]
        max_cumsum_all[i, 1] = abs(cumsum_sigma_sorted_by_i[int(divide_idx_all[i, 1])])
    argmax_max_cumsum_all = np.argmax(max_cumsum_all)
    divide_dim = argmax_max_cumsum_all // 2
    plus_minus = argmax_max_cumsum_all % 2
    divide_idx = int(divide_idx_all[divide_dim, plus_minus])
    thres = np.sort(x[:, divide_dim])[divide_idx]
    sigma_pred = (2 * plus_minus - 1) * np.ones(n)
    sigma_pred[x[:, divide_dim] <= thres] = 1 - 2 * plus_minus
    return sum(sigma_pred * sigma) / n

# ラベルをT回ふりなおして経験ラデマッハ複雑度をモンテカルロ法により計算する関数
def empirical_rademacher(x, T):
    n, d = x.shape
    score_each = [score_decision_stump(x, 1 - 2 * np.random.binomial(1, 0.5, n)) for t in range(T)]
    return np.mean(score_each)
```

　経験ラデマッハ複雑度を$d$と$n$を変えて算出した結果が以下です．

```python
for d in [1, 2, 10, 100]:
    for n in [100, 1000, 10000]:
        x = np.random.multivariate_normal([0]*d, cov = np.eye(d), size = n)
        print('d:', d, 'n:', n)
        print('emiprical rademacher:', empirical_rademacher(x, 100))
        print('upper bound:', np.sqrt(2*np.log(2*(n+1)*d)/n))
```

```python
d: 1 n: 100
emiprical rademacher: 0.1696
upper bound: 0.3258302532731178
d: 1 n: 1000
emiprical rademacher: 0.0537
upper bound: 0.12330370602601663
d: 1 n: 10000
emiprical rademacher: 0.018840000000000003
upper bound: 0.044505252605813755
d: 2 n: 100
emiprical rademacher: 0.1862
upper bound: 0.34645100311475935
d: 2 n: 1000
emiprical rademacher: 0.0637
upper bound: 0.1288025554128109
d: 2 n: 10000
emiprical rademacher: 0.019871999999999997
upper bound: 0.04603636546926007
d: 10 n: 100
emiprical rademacher: 0.2566
upper bound: 0.39015004268602227
d: 10 n: 1000
emiprical rademacher: 0.08428
upper bound: 0.1407443572785013
d: 10 n: 10000
emiprical rademacher: 0.02648
upper bound: 0.04940885070618524
d: 100 n: 100
emiprical rademacher: 0.32899999999999996
upper bound: 0.4452738007875446
d: 100 n: 1000
emiprical rademacher: 0.106
upper bound: 0.15625026173330564
d: 100 n: 10000
emiprical rademacher: 0.033699999999999994
upper bound: 0.053867908319377975
```

　結果から，次元$d$が大きくなると，経験ラデマッハ複雑度は大きくなり，サンプル数$n$が大きくなると，経験ラデマッハ複雑度は小さくなっていることが分かると思います．次元$d$が大きくなると，境界のひき方の選択肢が増え，モデルが複雑になるため，経験ラデマッハ複雑度が大きくなると考えられます．また．サンプル数$n$が大きくなると，綺麗に分類するのが難しくなるため，経験ラデマッハ複雑度が小さくなると考えられます．

　実は，決定株の経験ラデマッハ複雑度は，入力の従う分布によらず，$n \geq d$のとき，

$$
\hat{\mathfrak{R}}_{S}(\mathcal{G}) \leq \sqrt{\frac{2}{n}\log(2(n+1)d)}
$$

と上から抑えられます．上の結果にupper boundとして示しています．上から抑えられていることが分かると思います．この例では，そこまで細かくない評価ですね．また，この式からも次元$d$が大きくなると，経験ラデマッハ複雑度は大きくなり，サンプル数$n$が大きくなると，経験ラデマッハ複雑度は小さくなることが分かります．

# 2 ラデマッハ複雑度

## 2.1 ラデマッハ複雑度

　経験ラデマッハ複雑度では，サンプルは固定していましたが，経験ラデマッハ複雑度をサンプルの現れ方で期待値をとったもの

```math
\mathfrak{R}_n(\mathcal{G}) = \mathbb{E}_{S}\left[ \hat{\mathfrak{R}}_{S}(\mathcal{G}) \right]
```

が**ラデマッハ複雑度**です．

## 2.2 決定株を用いた実験

　$d$次元標準正規分布から，サンプルを繰り返し取り直して経験ラデマッハ複雑度を平均して計算したラデマッハ複雑度が以下です．

```python
for d in [1, 2, 10, 100]:
    for n in [100, 1000, 10000]:
        rademacher = 0
        empirical_rademacher_each = [empirical_rademacher(np.random.multivariate_normal([0]*d, cov = np.eye(d), size = n), 100)]
        print('d:', d, 'n:', n)
        print('rademacher:', np.mean(empirical_rademacher_each))
        print('upper bound:', np.sqrt(2*np.log(2*(n+1)*d)/n))
```

```python
d: 1 n: 100
rademacher: 0.16879999999999992
upper bound: 0.3258302532731178
d: 1 n: 1000
rademacher: 0.054120000000000015
upper bound: 0.12330370602601663
d: 1 n: 10000
rademacher: 0.017466
upper bound: 0.044505252605813755
d: 2 n: 100
rademacher: 0.20300000000000004
upper bound: 0.34645100311475935
d: 2 n: 1000
rademacher: 0.06462
upper bound: 0.1288025554128109
d: 2 n: 10000
rademacher: 0.019666000000000003
upper bound: 0.04603636546926007
d: 10 n: 100
rademacher: 0.26700000000000024
upper bound: 0.39015004268602227
d: 10 n: 1000
rademacher: 0.08355999999999998
upper bound: 0.1407443572785013
d: 10 n: 10000
rademacher: 0.02667599999999998
upper bound: 0.04940885070618524
d: 100 n: 100
rademacher: 0.324
upper bound: 0.4452738007875446
d: 100 n: 1000
rademacher: 0.10586
upper bound: 0.15625026173330564
d: 100 n: 10000
rademacher: 0.03384999999999999
upper bound: 0.053867908319377975
```

　ラデマッハ複雑度は経験ラデマッハ複雑度と同じ値により上から抑えることができます．つまり，

$$
\mathfrak{R}_{n}(\mathcal{G}) \leq \sqrt{\frac{2}{n}\log(2(n+1)d)}
$$

です．結果にupper boundとして示しました．上から抑えられていることが分かると思います．



# 3 なぜラデマッハ複雑度を考えるか

　なぜこのように定義が面倒なラデマッハ複雑度を考える必要があるのでしょうか．
　それは**一様大数の法則**という定理があり，この定理によりラデマッハ複雑度を用いて経験損失と予測損失の差を上から抑えられることが最も重要な理由だろうと思います．経験損失は手元にあるサンプルで評価した損失で，予測損失はサンプルの現れ方で損失の期待値をとったものです．本当に小さくしたいのは予測損失ですが，我々の手元には普通，有限個のサンプルしかないので経験損失しか評価できません．経験損失と予測損失の差を上から抑えることができれば，予測損失がどのぐらいか見積もることができます．

　ラデマッハ複雑度が小さいと経験損失と予測損失の差も小さくなります．決定株のラデマッハ複雑度が$\sqrt{\frac{2}{n}\log(2(n+1)d)}$で上から抑えられたことを思い出すと，$n$が大きくなると経験損失と予測損失の差が小さくなり，$d$，つまり決定株における複雑度が大きくなると経験損失と予測損失の差も大きくなるだろうと考えられます．複雑度が大きくなると経験損失と予測損失の差が大きくなることは，複雑度が大きくなるとバリアンスが大きくなることと合っていますね．

　上では決定株という非常に簡単なモデルについて考えましたが，『統計的学習理論』の中ではカーネル法やサポートベクタ―マシン，ブースティングのラデマッハ複雑度について説明されています．近年では，敵対的学習などについてもラデマッハ複雑度が調べられているそうです[2]．

# 参考文献

[1] 金森敬文．統計的学習理論．講談社，2015．
[2] https://www.slideshare.net/masakatooz/ss-177134363