Tree SHAP

決定木ベースの手法で高速にSHAPを計算する手法。

論文上のアルゴリズム名はTreeSHAP、shapパッケージの実装としては TreeExplainer

アルゴリズム¶

Shapley values¶

$N$ は全特徴量集合、 $M$ は特徴量数とする。特徴量集合 $S$ が観測されているときの条件付き期待値 $f_x(S)=E[f(x)\mid x_S]$ をとすると、Shapley valuesは

\phi_i = \sum_{S \subseteq N \setminus \{i\}} \frac{|S|!(M-|S|-1)!}{M!} [f_x(S \cup \{i\}) - f_x(S)]

(1)

となる。

Tree SHAPアルゴリズムは木構造に特化してShapley valuesを厳密に計算し、計算量を従来の $O(TL2^M)$ から $O(TLD^2)$ に削減できる。ここで $T$ は木の数、 $L$ は各木の最大葉数、 $D$ は最大の深さである。

上の式のうち $f_x(S)=E[f(x)\mid x_S]$ の推定が課題となる。

シンプルだが $O(TL2^M)$ のアルゴリズム¶

$f_x(S)=E[f(x)\mid x_S]$ を推定する単純で計算量の多いアルゴリズム。

Algorithm 1 Estimating

E[f(X) | do(X_S = x_S)]

def EXPVALUE(x, S, tree = {v, a, b, t, r, d}):
    def G(j):                # 再帰関数 G の定義開始
        if v_j != internal:  # ノードjが葉かどうか判定
            return v_j       # 葉の値を返す
        else:
            if d_j in S:     # この特徴量で条件付けしてるか確認
                # この木の経路の子ノードを使う
                return G(a_j) if x_{d_j} <= t_j else G(b_j)
            else
                # カバレッジrで重み付け
                return [G(a_j) * r_{a_j} + G(b_j) * r_{b_j}] / r_j
    return G(1) # root nodeから実行

ここで、

$v$ はノード値のベクトル。内部ノード（葉じゃないノード）に対しては internal という値を取る
ベクトル $a$ と $b$ は、それぞれ各内部ノードにおける左子ノードと右子ノードのインデックスを表す
ベクトル $t$ は各内部ノードにおける閾値
$d$ は、内部ノードで分割に使われる特徴量のインデックスのベクトル
ベクトル $r$ は各部分木にどれだけのデータサンプルが落ちるか（cover）を表している。

出所：Lundberg et al. (2020)

例¶

次のような木を考える

root: x_1 <= 5 ?
├── left leaf:  v = 10, cover = 30
└── right node: cover = 70
    ├── x_2 <= 3 leaf: v = 20,  cover = 50
    └── x_2 > 3  leaf: v = 100, cover = 20

説明対象のインスタンスが $x_1 = 7, x_2 = 2$ だとする

$S = \{1, 2\}$ の場合：
$x_1 <= 5$ → right node
→ $x_2 <= 3$ で $v=20$ → $E[f(X) \mid X_1=7, X_2=2]=20$

$S = {1}$ の場合：
$x_1 <= 5$ → right node
→ $x_2$ が未知なのでカバー比率 $r$ で子ノードを荷重和

E[f(X) \mid X_1=7] =\frac{50}{70} \cdot 20+\frac{20}{70} \cdot 100 = \frac{3000}{70} \approx 42.86

(2)

特徴量の数 $M$ に対し、部分集合 $S$ は $2^M$ 個あるので計算量が多いのが問題

$O(TLD^2)$ のアルゴリズム¶

決定木では、ある葉の値に影響するのは、基本的にはその葉に至るパス上に現れた特徴量だけなので、ほかの特徴についてのパスは無視できる → 深さ $D$ に対し $O(TLD^2)$

Algorithm 2 Tree SHAP with path dependent feature perturbation

def TREESHAP_PATH(x, tree={v, a, b, t, r, d}):
    # x: SHAP値を計算する対象サンプル
    # tree: 木構造。v=ノード値, a/b=子ノード, t=閾値, r=カバー数, d=分割特徴量
    φ = array of len(x) zeros  # 各特徴量のSHAP値

    def RECURSE(j, m, pz, po, pi):
        """メインの処理となる再帰関数
        
        引数：
          j: 現在のノード
          m: ルートから現在ノードまでのパス情報
          pz: pi が「ない」とき、この経路を通る割合
          po: pi が「ある」とき、この経路を通る割合
          pi: 今回パスに追加する特徴量
        """
        m = EXTEND(m, pz, po, pi)  # mにpiの情報を追加

        if v_j != internal:  # 葉ノードなら寄与を計算
            for i ← 2 to len(m):
                # i: パス上の特徴量インデックス。1番目はダミーなので飛ばす

                w = sum(UNWIND(m, i).w)
                # w: 特徴量 m_i.d を除いたときのShapley重み

                φ_{m_i} ← φ_{m_i} + w * (m_i.o - m_i.z) * v_j
                # m_i.d: 寄与を更新する特徴量
                # m_i.o - m_i.z: 特徴量が「ある/ない」で葉への到達割合がどれだけ変わるか
                # v_j: 葉の出力値

        else:  # 内部ノードなら子ノードへ再帰
            h, c = (a_j, b_j) if x_{d_j} <= t_j else (b_j, a_j)
            # h: hot child。x が実際に進む子ノード
            # c: cold child。x が実際には進まない子ノード

            iz = io = 1  # 既存の zero/one fraction の初期値

            k = FINDFIRST(m, d_j)
            # k: 特徴量 d_j がすでにパス上にある場合の位置

            if k != nothing:
                iz, io = (m_k.z, m_k.o)  # 既存の割合を引き継ぐ
                m = UNWIND(m, k)         # 重複する特徴量の古い拡張を巻き戻す

            RECURSE(h, m, iz * r_h / r_j, io, d_j)
            # hot child:
            #   特徴量 d_j がない場合はカバー比 r_h/r_j で分岐
            #   特徴量 d_j がある場合は x が進むので one fraction を維持

            RECURSE(c, m, iz * r_c / r_j, 0, d_j)
            # cold child:
            #   特徴量 d_j がない場合はカバー比 r_c/r_j で分岐
            #   特徴量 d_j がある場合は x が進まないので one fraction は0

    def EXTEND(m, pz, po, pi):
        """m に特徴量 pi を追加し、部分集合サイズごとの重みを更新する"""
        l, m = len(m), copy(m)
        # l: 追加前のパス長

        m_{l+1}.(d, z, o, w) = (pi, pz, po, (1 if l == 0 else 0))
        # d: 特徴量番号
        # z: zero fraction
        # o: one fraction
        # w: Shapley重み

        # 後ろから更新して、更新済みの重みを再利用しないようにする
        for i ← l to 1 do:
            # pi を部分集合に含める場合の重み
            m_{i+1}.w = m_{i+1}.w + po * m_i.w * (i / l)
            # pi を部分集合に含めない場合の重み
            m_i.w = pz * m_i.w * ((l - i) / l)
        return m

    def UNWIND(m, i):
        """m から i 番目の特徴量を取り除いた状態の重みを復元する"""
        l, n = len(m), m_l.w, copy(m_{1...l-1})
        # l: 現在のパス長
        # n: 重み復元用の一時変数
        for j ← l - 1 to 1 do:
            # 後ろから前へ重みを復元する
            if m_i.o != 0:
                t = m_j.w   # 更新前の重みを退避
                m_j.w = n * l / (j * m_i.o)
                n = t - m_j.w * m_i.z * ((l - j) / l)
            else:
                # po=0 の場合は one fraction で割れないので、zero fraction 側から復元
                m_j.w = (m_j.w * l) / (m_i.z * (l - j))
            for k ← j - 1 to 1 do:
                m_j.(d, z, o) = m_{j+1}.(d, z, o)
                # 取り除いた特徴量の分だけ、後続要素を前にずらす
        return m

    # 根ノードから開始。空パス、zero/one fraction はどちらも1
    RECURSE(1, [], 1, 1, 0)
    return φ

Algorithm 2 は、決定木の各葉に到達する確率が、ある特徴量を条件付けるかどうかでどれだけ変わるかを計算し、その差に葉の出力値を掛けて特徴量ごとの寄与に分配するアルゴリズム

木を根から葉まで再帰的にたどる
各分岐で、対象サンプル x が進む方向を hot child とする
特徴量が「ある」場合は hot child にだけ進む
特徴量が「ない」場合は、訓練データのカバー比 r_child / r_parent で左右に分配する
パス上では、特徴量ごとの zero fraction, one fraction, Shapley重みを管理する
葉に到達したら、w * (one fraction - zero fraction) * 葉の値 を各特徴量のSHAP値に加算する

変化

φ_{m_i} ← φ_{m_i} + w * (m_i.o - m_i.z) * v_j
# w: 重み
# (m_i.o - m_i.z): 特徴の条件付けでの葉の到達確率の変化
# v_j: 葉ノードでの出力値

TreeSHAP¶

解説¶

SHAP interaction values¶

TreeSHAPと同じ論文で提案された、特徴量同士の交互作用の強さを評価する手法

主効果から相互作用効果を分離することで、さらに追加の洞察を得ることができる。ペアごとの相互作用を考えると、あるモデル予測に対して、すべての特徴量ペアが与える影響を表す寄与度行列が得られる。SHAP値はゲーム理論における古典的なShapley値に基づいているため、相互作用効果への自然な拡張は、より現代的なShapley interaction index によって得られる。

\Phi_{i,j} = \sum_{S \subseteq N \setminus {i,j}} \frac{|S|!(M-|S|-2)!}{2(M-1)!} \nabla_{ij}(S)

(3)

ただし $i \ne j$ であり、

\begin{aligned} \nabla_{ij}(S) &= f_x(S \cup {i,j}) - f_x(S \cup {i}) - f_x(S \cup {j}) + f_x(S)\\ &= f_x(S \cup {i,j}) - f_x(S \cup {j}) - [f_x(S \cup {i}) - f_x(S)] \end{aligned}

(4)

式3では、特徴量 $i$ と特徴量 $j$ の間のSHAP相互作用値は、それぞれの特徴量に等しく分配される。そのため、

\Phi_{i,j} = \Phi_{j,i}

(5)

となり、全体の相互作用効果は

\Phi_{i,j} + \Phi_{j,i}

(6)

で表される。ある予測に対する主効果は、その特徴量のSHAP値から、その特徴量に関するSHAP相互作用値を差し引いたものとして定義できる。

\Phi_{i,i} = \phi_i - \sum_{j \ne i} \Phi_{i,j}

(7)

SHAP相互作用値は直接計算することもできるが、木モデルではAlgorithm 2を利用することで、計算コストを大幅に削減できる。SHAP相互作用値は、特徴量 $j$ が存在している場合の特徴量 $i$ のSHAP値と、特徴量 $j$ が存在していない場合の特徴量 $i$ のSHAP値との差として解釈できる。

このため、Algorithm 2を2回使うことができる。1回目は特徴量 $j$ が存在しているものとして固定して無視し、2回目は特徴量 $j$ が存在しないものとして計算する。各特徴量についてこの処理を繰り返すため、実行時間は $O(TMLD^2)$ となる。

アルゴリズム¶

Shapley values¶

シンプルだが O(TL2M)O(TL2^M)O(TL2M) のアルゴリズム¶

例¶

O(TLD2)O(TLD^2)O(TLD2) のアルゴリズム¶

TreeSHAP¶

解説¶

SHAP interaction values¶

シンプルだが $O(TL2^M)$ のアルゴリズム¶

$O(TLD^2)$ のアルゴリズム¶