【論文紹介】Interpretable Predictions of Tree-based Ensembles via Actionable Feature Tweaking【KDD 2017】
どうも，接点QB です．
KDD 2017の論文を漁っていたら，「Interpretable Predictions」というタイトルに目を惹かれて表題の論文を読んだので紹介したいと思います． 機械学習（教師あり）の手法は「予測・分類」を目的としていて，「説明性」に関してはあまり重視されない印象があります． しかし，実際の問題を解決しようとする際には説明性を求められるケースが多いと思います． たとえば，「この人は商品Aを購入しないグループに属する事はわかりました．じゃあ，どうすればこの人に商品Aを買わせることができますか？」といったケースがあります．
今回紹介する論文で提案されている手法は，上記のような例でも対応することが出来ます．
論文の概要
問題設定と前提条件は，
- 2値分類
- 使用するのは決定木ベースのアンサンブル分類器
です．そして，論文の目的は
- negativeに分類された$\renewcommand{\vec}[1]{\boldsymbol{#1}}\vec{x}$を変換して，positiveに分類される$\boldsymbol{x}'$を作成する
ことです．
Notation
- 特徴量ベクトル
- $\boldsymbol{x}\in \mathcal{X}\subset \mathbb{R}^n$．$\mathbb{E}[\vec{x}]=\vec{0},\, V[\vec{x}]=\vec{1}$とする．
- ラベル
- $\boldsymbol{x}$に対して$y\in \left\{-1, +1 \right\}=\mathcal{Y}$が対応する．
- 特徴量ベクトルとラベルの対応を表す写像
- $f:\mathcal{X}\to \mathcal{Y}$．分類器はこの写像を推定したもの$\widehat{f}$．
- アンサンブル学習器
- 弱学習器（決定木）$T_1,\cdots, T_K$のアンサンブル学習器を$\mathcal{T}$とする．
手法
まず注意してほしいのは，本稿で紹介する手法は「分類器を作成すること」ではなく「作成した分類器を使って特徴量ベクトルをより良いものに変換すること」です． なので，分類器$\mathcal{T}$は既に学習されているとします．
$\boldsymbol{x}$がtrue negativeすなわち$f(\boldsymbol{x})=\widehat{f}(\boldsymbol{x})=-1$を満たすとします． このとき，取り組むタスク \begin{equation} \boldsymbol{x}'\in \mathcal{X} \quad \text{such that} \quad \widehat{f}(\boldsymbol{x})=+1 \end{equation} を作成することです．さらに言えば，変換されたベクトルの中でも何かしらの意味で最良のベクトルを選択するという問題に落とし込むと，求めたいベクトルは一意に定まります：
\begin{align} \newcommand{\argmin}{\mathop{\rm arg\, min}\limits} \boldsymbol{x}' = \argmin_{x^{*}} \left\{ \delta \left( \boldsymbol{x}, \boldsymbol{x}'\right) \left| \widehat{f}(\boldsymbol{x})=-1 \text{ and } \widehat{f}\left( \boldsymbol{x}^{*}\right)=+1 \right. \right\}. \tag{1} \end{align}
(1)での$\delta$は，「何かしらの基準」を測るためのコスト関数です．ユークリッド距離でもコサイン類似度でも，目的によって使い分けるのが良いと思います．
Positive and Negative Paths
決定木の根から葉ノードへの経路は，一つの特徴量$x_i$に関する不等号条件の積み重ねによって形成されます． そして，葉ノードでnegativeかpositiveを分類します．
決定木$T_k$において$\vec{x}$をpositiveと分類する$j$番目の経路を$p_{k,j}^+$で表します． 簡単のため，$p_{k,j}$は次のように表せると仮定します： \begin{equation} p_{k,j} = \left\{ \left( x_1 \substack{\geq \\ \leq}\theta_1\right), \cdots ,\left( x_n \substack{\geq \\ \leq}\theta_n\right) \right\}. \end{equation} さらに，$P_{k}^{+}=\bigcup_{j\in T_k} p_{k,j}$とします．negativeの場合の$P_{k}^{-}$も同様に定義します．
さて，このパスを指定することで入力空間から特定の領域（$p^+_{k,j}$を指定すればpositiveに分類される領域）を切り出す事ができます． たとえば，$\vec{x}=[x_1, x_2]$という特徴量ベクトルを分類する事を考えます． このとき，１つ目のノードの条件が$x_1\leq \theta_1$，２つ目のノードでの条件が$x_2\leq \theta_2$だとすると， これらの条件から$\mathcal{X}_{\rm posi}=\left\{ (x_1, x_2) \left| x_1\leq \theta_1,\, x_2\leq \theta_2\right. \right\}$という領域がpositiveと分類される$\vec{x}$を含む領域である事がわかります． 図1で緑になっている領域が$\mathcal{X}_{\rm posi}$に対応します．
特徴量ベクトルの構成方法
[Tolomei et al., 2016] で提案されている手法は，$\vec{x}$から近くて，$\widehat{f}\left( \vec{x}'\right)=1$となるような$\vec{x}'$を発見する方法です．つまり(1)の$\vec{x}'$を見つけるという問題になるわけですね． その具体的な方法を説明していきます． まずは$p_{k,j}^{+}$の$\varepsilon$-satisfactory instance \begin{equation} \vec{x}^{+}_{j(\varepsilon)}[i] = \begin{cases} \theta_i - \varepsilon & x_i \leq \theta_{i} \\ \theta_i + \varepsilon & x_i > \theta_{i} \end{cases} \tag{2} \end{equation} を定義します． (2) において，$\vec{x}^{+}_{j(\varepsilon)}[i]$は$\vec{x}^{+}_{j(\varepsilon)}$の第$i$要素です． また，決定木$T_k$における$\vec{x}_{j(\varepsilon)}^{+}$の集合を$\Gamma_{k}=\bigcup_{j\in P_{k}^{+}}\vec{x}_{j(\varepsilon)}^{+}$で定義します． 注意してほしいのは，$\vec{x}_{j(\varepsilon)}^{+}\in \Gamma_{k}$は必ずしも$\widehat{f}\left( \vec{x}_{j(\varepsilon)}^{+}\right)=+1$とはならないことです． あくまで$\widehat{h}_k\left(\vec{x}_{j(\varepsilon)}^{+}\right)=+1$ということしか言えません（$\widehat{h}$は決定木$k$での予測ラベル）．
そして，$\Gamma=\bigcup_{k=1}^{K}\Gamma_k$（i.e. すべての決定木（弱学習器）での$\varepsilon$-satisfactory instance の集合）に属し $\widehat{f}\left(\vec{x}_{j(\varepsilon)}^{+}\right)=+1$を満たす$\vec{x}_{j(\varepsilon)}^{+}$の中から， $\delta\left(\vec{x}, \vec{x}_{j(\varepsilon)}^{+}\right)$を最小にするものを$\vec{x}'$とします． つまり， \begin{equation} \vec{x}'=\argmin_{\substack{\vec{x}_{j(\varepsilon)}^{+}\in \Gamma \\ \widehat{f}\left(\vec{x}_{j(\varepsilon)}^{+}\right)=+1}} \delta\left(\vec{x}, \vec{x}_{j(\varepsilon)}^{+}\right) \end{equation} によって，$\vec{x}'$を決定します．
$\vec{x}^{+}_{j(\varepsilon)}$を決定するために探索する範囲は，決定木の数$K$とpositiveパスの本数の$|P_{k}^+|$の積で書けるので $O(K|P_{k}^+|)$になります．
実装 ―Pythonでやってみる―
（注）2017/10/25に下記内容を変更しました．デバッグが出来ていないので，もしかしたらバグが含まれているかもしれません
擬似コードは論文に記載されているので，ここではPythonで本手法を使ってみたいとおもいます． とりあえずはirisデータでやってみます． これまでの説明は2クラスの分類問題で説明してきましたが，簡単に3クラス問題に拡張できるので，3クラス分類でやりたいと思います．
まずはデータを読み込んで，Random Forestで分類器を作ります．
import numpy as np import pandas as pd import scipy.stats import copy from sklearn import datasets from sklearn.ensemble import RandomForestClassifier import copy import scipy.stats iris = datasets.load_iris() x_arr = iris['data'] mean_x = x_arr.mean(axis=0) std_x = x_arr.std(axis=0) x_arr = scipy.stats.zscore(x_arr) y_arr = iris['target'] rfc = RandomForestClassifier() rfc.fit(x_arr, y_arr)
rfcに各弱学習器（分類器）の情報が含まれていて取り出すことが出来ます．
たとえば，rfc[0]には
DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=None,
max_features='auto', max_leaf_nodes=None,
min_impurity_decrease=0.0, min_impurity_split=None,
min_samples_leaf=1, min_samples_split=2,
min_weight_fraction_leaf=0.0, presort=False,
random_state=2072047594, splitter='best')
が入っています．そして，この決定木に関する情報を取り出して目的のラベルに至るパスや，各ノードの閾値を得なければなりません．
決定木に関する情報は，オブジェクトの.tree_メンバが持っています．
ここで，まずはwc = rfr[0]の葉ノードから根ノードに至るパスを所得する方法を考えましょう．
ここで使えそうなtree_メンバでchildren_leftとchildren_rightがあります．
children_leftの第$i$要素には，ノードIDが$i$の左側子ノードのIDが入っています．
children_rightも同様で，こちらには右側子ノードのIDが入っています．
もちろん，葉ノードには子ノードが存在しないので，ノード$i$が葉ノードになっているときは$-1$が入っています．
試しにwc.tree_.children_leftを実行してみると
array([ 1, -1, 3, -1, 5, -1, 7, 8, -1, 10, -1, -1, -1])
が帰ってきます．しがたって，子ノードが$-1$になっている要素$i$を取ってくれば，子ノードのIDが分かります．
children_left = wc.tree_.children_left # information of left child node # select the top-maxj prediction node IDs leaf_nodes = np.where(children_left==-1)[0]
今回ほしいのはpositiveラベル（ここではラベル$3$にします）を出力する葉ノードへのパスだけなので，
leaf_nodesでラベル$3$を出力するものだけを選択します．
class_labels = [0, 1, 2] aim_label = 2 leaf_values = wc.tree_.value[leaf_nodes].reshape(len(leaf_nodes), len(class_labels)) leaf_nodes = np.where(leaf_values[:, aim_label] != 0)[0]
子ノードIDがわかったので，適当に子ノードのIDを一つ固定して話を進めます；
leaf_node = leaf_nodes[0]．
では，この葉ノードから根ノードまでのパスを求めます．
children_right = wc.tree_.children_right feature = wc.tree_.feature threshold = wc.tree_.threshold leaf_node = leaf_nodes[0] # 葉ノードを固定 child_node = leaf_node parent_node = -100 # initialize parents_left = [] # 左側親ノード parents_right = [] # 右側親ノード while (parent_node != 0): if (np.where(children_left == child_node)[0].shape == (0, )): parent_left = -1 # 左側親ノードが存在しない場合は-1 parent_right = np.where( children_right == child_node)[0][0] parent_node = parent_right elif (np.where(children_right == child_node)[0].shape == (0, )): parent_right = -1 # 右側親ノードが存在しない場合は-1 parent_left = np.where(children_left == child_node)[0][0] parent_node = parent_left parents_left.append(parent_left) parents_right.append(parent_right) """ 次のステップへの処理 """ child_node = parent_node
上記プログラムですが，まずはchild_node = leaf_nodeによって子ノードIDを初期化しています．
また，parent_node = -100で親ノードIDを初期化しています．
parent_leftは少しややこしくて，「ノード$i$を左側子ノードにもつノードIDが，第$i$要素に格納」されています（これを左側親ノードと呼ぶことにします）;
parent_right = np.where(children_right == child_node)[0][0]．
parent_rightも同様で，「ノード$i$を右側子ノードにもつノードIDが，第$i$要素に格納」されています．
parent_left, parent_right共に，条件を満たすノードが存在しなければ$-1$が格納されます；parent_left=-1．
決定木ではparent_leftの第$i$要素とparent_rightの第$i$要素に非$-1$である値が同時に入ることはありません
（i.e. 特定のノードへ2つ以上のノードから矢印が伸びることはない）．
なので，親ノードはwhile文の中のif文によって一意に決まります；parent_node=parent_right or parent_node=parent_left．
この親ノードを辿っていく処理をparent_nodeが$0$になるまで繰り返すことによって（根ノードのIDは必ず$0$），
葉ノードから根ノードまでのパスを求める事ができます．これらの処理を全ての葉ノードに対して行います．
""" search the path to the selected leaf node """ paths = {} for leaf_node in leaf_nodes: """ correspond leaf node to left and right parents """ child_node = leaf_node parent_node = -100 # initialize parents_left = [] # 左側親ノード parents_right = [] # 右側親ノード while (parent_node != 0): if (np.where(children_left == child_node)[0].shape == (0, )): parent_left = -1 # 左側親ノードが存在しない場合は-1 parent_right = np.where( children_right == child_node)[0][0] parent_node = parent_right elif (np.where(children_right == child_node)[0].shape == (0, )): parent_right = -1 # 右側親ノードが存在しない場合は-1 parent_left = np.where(children_left == child_node)[0][0] parent_node = parent_left parents_left.append(parent_left) parents_right.append(parent_right) """ 次のステップへの処理 """ child_node = parent_node # nodes dictionary containing left parents and right parents paths[leaf_node] = (parents_left, parents_right)
これで，pathsはkeyが葉ノードIDで，valueが（左側親ノードID，右側親ノードID）になっているディクショナリになります：
{1: ([0], [-1]),
3: ([2, -1], [-1, 0]),
5: ([4, -1, -1], [-1, 2, 0]),
8: ([7, 6, -1, -1, -1], [-1, -1, 4, 2, 0]),
10: ([9, -1, 6, -1, -1, -1], [-1, 7, -1, 4, 2, 0]),
11: ([-1, -1, 6, -1, -1, -1], [9, 7, -1, 4, 2, 0]),
12: ([-1, -1, -1, -1], [6, 4, 2, 0])}
では次に，パス内のノードでの分岐条件を記述する情報を取り出します． 必要な情報は，条件分岐に使われる
- 特徴量のインデックス；
features - 閾値；
thresholds - 不等号の向き；
inequality_symbols
です．
path_info = {}
for i in paths:
node_ids = [] # node ids used in the current node
# inequality symbols used in the current node
inequality_symbols = []
thresholds = [] # thretholds used in the current node
features = [] # features used in the current node
parents_left, parents_right = paths[i]
for idx in range(len(parents_left)):
if (parents_left[idx] != -1):
""" the child node is the left child of the parent """
node_id = parents_left[idx] # node id
node_ids.append(node_id)
inequality_symbols.append(0)
thresholds.append(threshold[node_id])
features.append(feature[node_id])
elif (parents_right[idx] != -1):
""" the child node is the right child of the parent """
node_id = parents_right[idx]
node_ids.append(node_id)
inequality_symbols.append(1)
thresholds.append(threshold[node_id])
features.append(feature[node_id])
path_info[i] = {'node_id': node_ids,
'inequality_symbol': inequality_symbols,
'threshold': thresholds,
'feature': features}
各パスに対する処理の中（for文の中）で，ノードが左側親ノードを持つ場合と右側親ノードをもつ場合に分けて考えています．
idxの初期値は葉ノードで，その親ノードに関する情報からnode_ids, inequality_symbols, thresholds, featuresに格納されていきます．
そして，最終的にディクショナリpath_infoに葉ノードをkeyとしたディクショナリが格納されます．
一旦，ここまでの処理を関数にします．
def search_path(estimator, class_labels, aim_label): """ return path index list containing [{leaf node id, inequality symbol, threshold, feature index}]. estimator: decision tree maxj: the number of selected leaf nodes """ """ select leaf nodes whose outcome is aim_label """ children_left = estimator.tree_.children_left # information of left child node children_right = estimator.tree_.children_right feature = estimator.tree_.feature threshold = estimator.tree_.threshold # leaf nodes ID leaf_nodes = np.where(children_left == -1)[0] # outcomes of leaf nodes leaf_values = estimator.tree_.value[leaf_nodes].reshape(len(leaf_nodes), len(class_labels)) # select the leaf nodes whose outcome is aim_label leaf_nodes = np.where(leaf_values[:, aim_label] != 0)[0] """ search the path to the selected leaf node """ paths = {} for leaf_node in leaf_nodes: """ correspond leaf node to left and right parents """ child_node = leaf_node parent_node = -100 # initialize parents_left = [] # 左側親ノード parents_right = [] # 右側親ノード while (parent_node != 0): if (np.where(children_left == child_node)[0].shape == (0, )): parent_left = -1 # 左側親ノードが存在しない場合は-1 parent_right = np.where( children_right == child_node)[0][0] parent_node = parent_right elif (np.where(children_right == child_node)[0].shape == (0, )): parent_right = -1 # 右側親ノードが存在しない場合は-1 parent_left = np.where(children_left == child_node)[0][0] parent_node = parent_left parents_left.append(parent_left) parents_right.append(parent_right) """ for next step """ child_node = parent_node # nodes dictionary containing left parents and right parents paths[leaf_node] = (parents_left, parents_right) path_info = {} for i in paths: node_ids = [] # node ids used in the current node # inequality symbols used in the current node inequality_symbols = [] thresholds = [] # thretholds used in the current node features = [] # features used in the current node parents_left, parents_right = paths[i] for idx in range(len(parents_left)): if (parents_left[idx] != -1): """ the child node is the left child of the parent """ node_id = parents_left[idx] # node id node_ids.append(node_id) inequality_symbols.append(0) thresholds.append(threshold[node_id]) features.append(feature[node_id]) elif (parents_right[idx] != -1): """ the child node is the right child of the parent """ node_id = parents_right[idx] node_ids.append(node_id) inequality_symbols.append(1) thresholds.append(threshold[node_id]) features.append(feature[node_id]) path_info[i] = {'node_id': node_ids, 'inequality_symbol': inequality_symbols, 'threshold': thresholds, 'feature': features} return path_info
次に，$\varepsilon$-satisfactory instance を計算する部分を作ります． これはそんなに難しくないので一気に書いて関数としてまとめておきます．
def esatisfactory_instance(x, epsilon, path_info): """ return the epsilon satisfactory instance of x. """ esatisfactory = copy.deepcopy(x) for i in range(len(path_info['feature'])): # feature index feature_idx = path_info['feature'][i] # threshold used in the current node threshold_value = path_info['threshold'][i] # inequality symbol inequality_symbol = path_info['inequality_symbol'][i] if inequality_symbol == 0: esatisfactory[feature_idx] = threshold_value - epsilon elif inequality_symbol == 1: esatisfactory[feature_idx] = threshold_value + epsilon else: print('something wrong') return esatisfactory
最後に，上で定めた2つの関数を使って提案手法を実装します．
def feature_tweaking(ensemble_classifier, x, class_labels, aim_label, epsilon, cost_func): """ This function return the active feature tweaking vector. x: feature vector class_labels: list containing the all class labels aim_label: the label which we want to transform the label of x to """ """ initialize """ x_out = copy.deepcopy(x) # initialize output delta_mini = 10**3 # initialize cost for estimator in ensemble_classifier: if (ensemble_classifier.predict(x.reshape(1, -1)) == estimator.predict(x.reshape(1, -1)) and estimator.predict(x.reshape(1, -1) != aim_label)): paths_info = search_path(estimator, class_labels, aim_label) for key in paths_info: """ generate epsilon-satisfactory instance """ path_info = paths_info[key] es_instance = esatisfactory_instance(x, epsilon, path_info) if estimator.predict(es_instance.reshape(1, -1)) == aim_label: if cost_func(x, es_instance) < delta_mini: x_out = es_instance delta_mini = cost_func(x, es_instance) else: continue return x_out
まとめ
よく「ブラックボックス」と言われる機械学習ですが，今回紹介した手法を使うと （どうすればpositiveなラベル方向へ改善出来るのかという事がわかるという意味で）解釈性を上げる事ができます． また，Decision Treeのアンサンブル手法であるRandom ForestやGradient Boostingは非常に手軽に試せる手法で， かつ比較的高精度な分類器を用意に作ることが出来ます．そのような分類器において解釈性の向上が得られるということは， 実際のデータマイニング応用にとても役立つと思われます．
