colabでLIMEとSP-LIMEを動かす。

### はじめに

本記事では、機械学習を説明する手法LIMEおよびSP-LIMEについて解説して、`Google Colaboratory`で実際に動かします。
本記事は自分が某セミナーでLIMEの初出論文`"Why Should I Trust You?": Explaining the Predictions of Any Classifier`について発表したスライドをベースにしています。
https://irisu-inwl.github.io/reveal/lime

ただし、このスライドは機械学習を知らない人のために作った感があるので、改めて記事としてまとめました。
この記事でのコードは以下です。
https://colab.research.google.com/drive/16Dl5KHJ8ER3XvG4MMSMxQlE3VoaRWhp9

### google colabの準備

1. google colabを開いて新しいnotebookを作る。  
  https://colab.research.google.com/notebooks/welcome.ipynb?hl=ja
2. LIMEをインストール

```
!pip install lime
```

以上、これだけ！

### LIMEのモチベーション

機械学習（特にここでは教師あり学習）を利用するとき、学習モデルは予測した結果とその確率しか出してくれない。

たとえば、sklearnのwineデータを考える。

```python

from sklearn.datasets import load_wine
wine_data = load_wine()

from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score

# 正解の比率を変えずにデータをX_train:X_test = 4:1に分割（訓練するデータとテストするデータに分割）
X_train, X_test, y_train, y_test = train_test_split(wine_data.data, wine_data.target, test_size=0.2, stratify = wine_data.target)
# 学習
clf = RandomForestClassifier().fit(X_train, y_train)
# 予測
predict_data = clf.predict(X_test)
score = accuracy_score(predict_data, y_test)
print(f'正解率: {score}')
# 予測した値
print(predict_data)
```

- out

```
正解率: 1.0
[2 0 1 2 1 1 1 2 2 0 1 0 0 2 1 0 1 0 1 0 0 2 1 2 2 1 1 0 0 1 2 1 0 1 0 2]
```

教師あり学習では与えられたデータに対して学習モデルが予測したラベルを返してくれる。
このとき、機械学習を実用する上で以下の疑問を抱く。

1.  学習モデルを作り、そこから得られた予測結果は正しいのか？
2.  訓練データから学習した学習モデル自体は正しいのか？

これらの疑問について、1は**LIME(Local Interpretable Model-agnostic Explanations)**、2は**SP-LIME(Submodular Pick ...)**という手法で解決する。

- LIMEは予測するデータの周辺をサンプリングし、学習モデルにおける一つのラベルを予測する確率 $f(x)$ を、線形関数$g(x)=w^{T} \cdot x$で線形回帰することで、どの特徴が予測結果に寄与しているかを見せることで説明する。（寄与率は$w$となる）
  - このとき、$f(x)$は確率であればよいので、モデルに依存することなく説明できる。（model-agnostic）
- SP-LIMEは訓練データすべてをLIMEし、なるべく寄与率が大きく特徴を網羅するようなデータをいくつか訓練データから抜き出すことで、モデル自体がどのような判断をするかを説明する。
  - LIMEを使って訓練データの中からいくつか選び出すからPickである（Submodularは選択するために最適化する関数が劣モジュラ関数であるため）


### LIMEの解説

#### コード例
- アルゴリズムを理解できなくてもコードを動かすことはできる（推奨しないけど）ので、まずはそれで
- 先ほどのwineデータに対して予測の説明を行う

```python

from lime.lime_tabular import LimeTabularExplainer

train_data = wine_data.data
class_names = wine_data.target_names
feature_names = wine_data.feature_names

explainer = LimeTabularExplainer(train_data ,class_names=class_names, 
                                feature_names = feature_names, kernel_width=3, verbose=False)

# 正解の比率を変えずにデータをX_train:X_test = 4:1に分割
X_train, X_test, y_train, y_test = train_test_split(wine_data.data, wine_data.target, test_size=0.2, stratify = wine_data.target)
# 学習
clf = RandomForestClassifier().fit(X_train, y_train)

explain_data = X_test[0]
explain_result = explainer.explain_instance(explain_data, clf.predict_proba, num_features=6, top_labels=2)

explain_result.show_in_notebook(show_table=True, show_all=True)
```

- 結果

![lime_result.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299097/5b6d684f-b9c0-742a-5158-9c5f3bfd5636.png)

このように予測に対して、どの特徴量が寄与しているかを図示してくれる。

#### 予測モデルを説明するとは

学習モデルが予測した結果の説明を以下のように行う。

- $f(x):\mathbb{R}^d\to\mathbb{R}$をデータ$x$があるクラスに入ってる確率
  - 例えばwineデータを予測した際の確率関数は以下のようになる。
  
```python

print(clf.predict_proba(X_test))
```

- out

```
[[0.  0.1 0.9]
 [1.  0.  0. ]
 [0.  1.  0. ]
 [0.  0.  1. ]
 [0.  0.5 0.5]
 [0.  0.6 0.4]
 [0.1 0.9 0. ]
```

入力したデータがどのラベルであるかの確率がベクトル値として帰ってくるが、$f$はベクトルのある成分に固定したものだと考える。

- データのベクトル$x\in \mathbb{R}^d$に対して人が解釈しやすい空間へと変換した点$x^\prime \in \lbrace 0,1 \rbrace ^{d^\prime}$を考える
- 例えば、アルコール度数($\in\mathbb{R}$)を、`ある範囲からある範囲に入っていれば1となる`バイナリベクトルに変換する。

$x = 14.38 \in \mathbb{R}$を

<table>
<tr>
<th>$x \leq 12.36$</th>
<th>$12.36 < x \leq 13.05$</th>
<th>$13.05 < x \leq 13.68$</th>
<th>$13.68 < x$</th>
</tr>
<tr>
<td>0</td>
<td>0</td>
<td>0</td>
<td>1</td>
</tr>
</table>

というバイナリベクトル$(0,0,0,1)\in \lbrace 0,1 \rbrace^4$に変換する。これをすべての成分に行う。（要はBinning）

- $g: \mathbb{R}^{d^\prime}\to\mathbb{R}$: 解釈可能モデル, 　$G$: 解釈可能モデルの集合
- $\pi_x(z)$:$x$からどれくらい離れているかの尺度。$z$が$x$に近ければ小さい値を取る。
- $\Omega(g)$: 解釈するためのモデル$g$の複雑度
  - 例えば、線形モデル$g(x^\prime) = w\cdot {x^{\prime}}^T$だったら係数の非ゼロ成分数
- $L(f,g,\pi_x)$: 学習モデル$f$を説明する$g$が信頼できない尺度(loss)
- $Z$: $x,x^\prime$の周辺からサンプリングされた点の集合
- $f$を説明するモデルを求めるのは以下の式によって求められる。

$$
\xi(x)=argmin_{g\in G} \left\[ L(f,g,\pi_x)+\Omega(g)\right\]
$$

#### LIMEのアルゴリズム

LIMEでは上記した説明方法において特殊な条件で行っている。

- $g$は線形モデルのみを考えてる。つまり、$g(x^\prime)=w\cdot {x^{\prime}}^T$
  - $g$を求めることは即ちパラメータ$w$を求めることである。
- 離れている尺度は$\pi_x(z)=\exp(-distance(x,z)^2/\sigma^2)$
- 信頼の損失は$L(f,g,\pi_x)=\sum_{z,z^\prime\in Z} \pi_x(z)(f(z)-g(z^\prime))^2$
- $w$を求めるために指定された$K$成分以外は0とするK-Lassoという手法を用いる
  - gitのコード見るとRidge回帰してるんだけど…  
    https://github.com/marcotcr/lime/blob/master/lime/lime_base.py#L161

実際のアルゴリズムは以下（論文からの抜粋）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299097/c218bef7-840e-fdaa-7176-707e917f3af4.png)


サンプリングして線形回帰している様子がLIMEの解説でよく使われる以下の図で表されている（論文からの抜粋）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299097/49ba077d-10ce-6b82-9096-9a46d1c38337.png)

コードの例で出てきた特徴量の寄与率が、$w$となっている。+であればclassに所属している寄与率で、-であればclassに所属していないことを表すようになっている。

### SP-LIMEの解説

コードの例

```python

from lime import submodular_pick

clf = RandomForestClassifier().fit(wine_data.data, wine_data.target)
sp_obj = submodular_pick.SubmodularPick(explainer, wine_data.data, clf.predict_proba, sample_size=50, num_features=13, num_exps_desired=5)

[exp.as_pyplot_figure() for exp in sp_obj.sp_explanations];
```

- 結果

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299097/628bae8a-0d68-0b4d-3b10-58eb4955f542.png)

これだけ見せられても何もわからんって感じだが、データの中から代表的なものを選んで学習モデルを説明している。
でも、正直、分からない。

#### SP-LIMEのアルゴリズム
- $n$件のデータの集合$X$すべてにLIMEを行う。
  - すると説明した結果の$n\times d^\prime$行列$W=(\mathcal{W} _ {ij}) _ {i\leq n,j\leq d^\prime}$が得られる（$i$番目のデータについてLIMEを行って得られた$w_j$を並べたもの）
- $X$の中から、よく使われる特徴(LIMEのアウトプット$w$が大きいもの)を持っているデータを$B$個抜き出す。
  （実は実装では計算量のためなのか$X$からサンプリングした中から抜き出している。）
  - $I_j=\sqrt{\sum_{i\leq n}\mathcal{W} _ {ij}}$とし
$$
c(V,W,I)=\sum_{j \leq d^\prime} \mathbb{I} _ {\[\exists i\in V: \mathcal{W} _ {ij}>0\]} I_j
$$
  - すると、SP-LIMEでは以下の最適解を求めることになる。
$$
argmax_{|V|\leq B}c(V,W,I)
$$

- その結果、SP-LIMEのアルゴリズムは以下である。（論文からの抜粋）

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/299097/a9b70acb-b824-ac53-28fc-3ac71cc50912.png)


- 実は$argmax_{|V|\leq B}c(V,W,I)$はNP-hardである。
- 一般に非負単調劣モジュラ関数$f$について$argmax_{|V|\leq k}f(V)$を求めるNP-hardな問題は、greedy algorithmで近似解を得られることが分かっている。（参考: https://las.inf.ethz.ch/files/krause12survey.pdf）

とくに以下の定理(上記pdfのTheorem 1.5)が成り立つ。
非負単調劣モジュラ関数$f:P(V)\to \mathbb{R}_{+}$とし、$\lbrace S_i \rbrace$をgreedyに選択した集合とする。つまり、

$$
S_i := S_{i-1}\cup argmax_e \[f(S_{i-1}\cup e)-f(S_{i-1})\]
$$

このとき任意の$k,l$について

$$
f(S_l)\geq (1-e^{^\frac{l}{k}})max_{|S|\leq k}f(S).
$$

とくに$l=k$として、$f(S_k)\geq (1-\frac{1}{e})max _ {|S|\leq k}f(S)$となる。
つまり、greedy algorithmでの結果と最適解の比が
$$\frac{S_k}{max _ {|S|\leq k}f(S)}\geq 1-\frac{1}{e}\approx 0.63212055882$$
である。

- 劣モジュラ関数の結果を使って$c(V,W,I)$が$V$に対する非負単調劣モジュラ関数であれば、greedy algorithmで近似解を求められる。明らかだが論文では示されていないので示す。
- $c$の単調性: $A\subset B$について

```math
\begin{align*}
\sum_{j\leq d^\prime} \mathbb{I} _ {[\exists i\in A: \mathcal{W} _ {ij}>0]} 
&\leq \sum_{j\leq d^\prime} \mathbb{I} _ {[\exists i\in A: \mathcal{W} _ {ij}>0]}
 + \sum_{j\leq d^\prime, \forall i\in A, \mathcal{W} _ {ij}=0} \mathbb{I} _ {[\exists i\in B-A: \mathcal{W} _ {ij}>0]}\\
&=\sum_{j\leq d^\prime} \mathbb{I} _ {[\exists i\in B: \mathcal{W} _ {ij}>0]}
\end{align*}
```

なので単調である。

- $c$の劣モジュラ性

劣モジュラ関数は非負線形性を持つ。すなわち、劣モジュラ関数$g_j(V)$,$a_j\geq 0$に対して$\sum_j a_jg_j(V)$は劣モジュラである。
$c(V)=\sum_{j \leq d^\prime} \mathbb{I} _ {\[\exists i\in V: \mathcal{W} _ {ij}>0\]} I_j$に対して、$g_j(V)=\mathbb{I} _ {\[\exists i\in V: \mathcal{W} _ {ij}>0\]}$として、$I_j\geq 0$なので、$g_j(V)$が劣モジュラあれば、$c$の劣モジュラ性が示される。
いま、任意の$A,B$に対して、

|$\exists i\in A\cap B: \mathcal{W} _ {ij}>0$|$\exists i\in A\cup B: \mathcal{W} _ {ij}>0$|$\exists i\in A: \mathcal{W} _ {ij}>0$|$\exists i\in B: \mathcal{W} _ {ij}>0$|
|:--|:--|:--|:--|
|True|True|True|True|
|False|True|少なくともどちらかがTrue|少なくともどちらかがTrue|
|False|False|False|False|

となるので、

```math
\begin{align}
g_j(A\cap B)+g_j(A\cup B) &= \mathbb{I} _ {[\exists i\in A\cap B: \mathcal{W} _ {ij}>0]}+\mathbb{I} _ {[\exists i\in A\cup B: \mathcal{W} _ {ij}>0]} \\
&\leq \mathbb{I} _ {[\exists i\in A: \mathcal{W} _ {ij}>0]} + \mathbb{I} _ {[\exists i\in B: \mathcal{W} _ {ij}>0]} = g_j(A)+g_j(B)
\end{align}
```

よって劣モジュラ性は示された。

- $c$の非負単調劣モジュラであることは分かったので、SP-LIMEアルゴリズムの正当性は示された。

### 所感
- LIMEとSP-LIMEの中身を知ることができた。
  - LIMEはアイデア勝利なところがあるが、SP-LIMEはアルゴリズミックで読んでいて楽しかった。
- LIMEは線形モデルであるが、その他のモデル（決定木とか）の場合どうなるか気になる（初出論文のfuture workとなっている）
- また、$dom(g)=\lbrace 0,1\rbrace ^ {d^\prime}$が解釈可能なバイナリベクトルというデウスエクスマキナ感あるので、実用する上でその設計をどうするのか気になった。
  - githubにあるLIMEの実装では4種類の離散化クラスを用意している。  
    https://github.com/marcotcr/lime/blob/master/lime/discretize.py