Kernel SHAP

モデル非依存（model-agnostic）、つまり任意の機械学習アルゴリズムで作った予測モデルに対して適用できるSHAP values推定アルゴリズム。

SHAPの提案論文（Lundberg & Lee, 2017）のTheorem2あたりで出てくる話

数式の概要

Notation

• \(M\)時点の特徴量ベクトル\(x\)があるとし、特徴量が存在するかどうかを表すようにsimplifiedした二値変数のベクトルを\(s \in\{0,1\}^M\)とする。

• \(s\)を\(x\)に戻す関数を\(x = h_x(s)\)とおく

• 予測値を特徴\(i\)の寄与度\(\phi_i\) に分解する関数を\(g\)とおく： \(f(x)=g(x)=\phi_0+\sum_{i=1}^M \phi_i x_i\)

Shapley kernel

重み関数

\[ \pi_x(s) =\frac{M-1}{\binom{M}{|s|}|s|(M-|s|)} \]

をShapley kernelと呼ぶ。なお、ここで \(|s|\)は\(s\)の非ゼロ要素の数（\(|s|=\sum_{j=1}^M s_j\)）である。

重み付き二乗損失

モデルの予測値（実測値）\(f(h_x(s))\) と 説明モデルの出力\(g(s)\)を近づけるようにすると損失が小さくなる関数を定義し、Shapley kernelで重みをつける。

\[ L(f, g, \pi_x) =\sum_{s \in \{0,1\}^M}[f(h_x(s))-g(s)]^2 \pi_x(s) \]

Theorem2

上記の \(L(f, g, \pi_x)\) を最小化する解はShapeley values（モデルの説明について望ましい性質1~3を満たす値）である

（そしてこの重み付き回帰がKernel SHAPと呼ばれる）

実装例で確かめる

全体像

Simple Kernel SHAP — SHAP latest documentation

を参考にしつつ改変したものが次

（referenceの値の重要性はいまいちわかっておらずひとまず0にしている）

import itertools
import numpy as np
import scipy.special


def shapley_kernel(M, s):
    """重みカーネル"""
    if s == 0 or s == M:
        return 10000 # infが返るのを防ぐ
    return (M - 1) / (scipy.special.binom(M, s) * s * (M - s))


def powerset(M: int):
    """特徴量のすべての組み合わせを列挙する"""
    s = list(range(M))
    return itertools.chain.from_iterable(itertools.combinations(s, r) for r in range(len(s) + 1))


def create_data(x):
    M = x.shape[0] # 特徴ベクトルの次元数

    # 特徴量の出現の有無（0,1）を表す
    S = np.zeros((2**M, M + 1))
    S[:, -1] = 1 # 切片項 φ_0 のためのカラムを追加

    weights = np.zeros(2**M)

    # 特徴量の値を入れる行列
    reference = 0
    X = np.zeros((2**M, M))
    for i in range(2**M):
        X[i, :] = reference

    ws = {}
    for i, s in enumerate(powerset(M)):
        s = list(s)
        X[i, s] = x[s]
        S[i, s] = 1
        w = shapley_kernel(M, len(s))
        ws[len(s)] = ws.get(len(s), 0) + w
        weights[i] = w
    return X, S, weights


def solve_wls(X, W, y):
    """重み付き最小二乗法（WLS）を解く"""
    A = X.T @ W @ X
    b = X.T @ W @ y
    return np.linalg.solve(A, b)


def kernel_shap(x, f):
    """Kernel SHAPの簡単な実装例"""
    X, S, weights = create_data(x)
    y = f(X) # 予測モデルの出力値
    W = np.diag(weights)  # 重み行列
    phis = solve_wls(S, W, y) # NOTE: 入力の行列はX（特徴量の値）ではなくS（特徴量の出現の有無）
    return phis


# ---------------------------------
# 例：線形回帰モデルを説明する
# ---------------------------------
x = np.array([10, 15]) # 特徴量ベクトル

def f(x):
    """線形回帰モデル"""
    alpha = 10 # 切片
    beta = np.array([5, 10]) # 回帰係数
    return alpha + x @ beta

phi = kernel_shap(x, f)
base_value = phi[-1]
shap_values = phi[:-1]
print("          x =", x)
print("shap_values =", shap_values)
print(" base_value =", base_value)
print("   sum(phi) =", np.sum(phi))
print("       f(x) =", f(x))

          x = [10 15]
shap_values = [ 50. 150.]
 base_value = 10.000000000000002
   sum(phi) = 210.0
       f(x) = 210

詳細解説

予測モデル

単純な線形回帰を考える

x = np.array([10, 15]) # 特徴量ベクトル

def f(x):
    """線形回帰モデル"""
    alpha = 10 # 切片
    beta = np.array([5, 10]) # 回帰係数
    return alpha + x @ beta

f(x)

210

Shapley kernel（重み関数）

\[ \pi_x(s) =\frac{M-1}{\binom{M}{|s|}|s|(M-|s|)} \]

import scipy.special

def shapley_kernel(M: int, s: int) -> float:
    # sはここではlen(s)
    if s == 0 or s == M:
        return 10000.0 # infが返るのを防ぐ
    return (M - 1) / (scipy.special.binom(M, s) * s * (M - s))

# 例
shapley_kernel(M=2, s=1)

0.5

データ行列を作る

def powerset(M: int):
    """特徴量のすべての組み合わせを列挙する"""
    s = list(range(M))
    return itertools.chain.from_iterable(itertools.combinations(s, r) for r in range(len(s) + 1))

list(powerset(M=2))

[(), (0,), (1,), (0, 1)]

def create_data(x):
    M = x.shape[0] # 特徴ベクトルの次元数

    # 特徴量の出現の有無（0,1）を表す
    S = np.zeros((2**M, M + 1))
    S[:, -1] = 1 # 切片項 φ_0 のためのカラムを追加

    weights = np.zeros(2**M)

    # 特徴量の値を入れる行列
    reference = 0
    X = np.zeros((2**M, M))
    for i in range(2**M):
        X[i, :] = reference

    ws = {}
    for i, s in enumerate(powerset(M)):
        s = list(s)
        X[i, s] = x[s]
        S[i, s] = 1
        w = shapley_kernel(M, len(s))
        ws[len(s)] = ws.get(len(s), 0) + w
        weights[i] = w
    return X, S, weights

重み付き最小二乗法を解く

重み付き最小二乗法（weighted least squares: WLS）の推定量は、重みを持った対角行列\(W\)と入力の行列\(X\)、目的変数ベクトル\(y\)について

\[ \phi = (X^\top W X)^{-1} X^\top W y \]

と定義される。

def solve_wls(X, W, y):
    A = X.T @ W @ X
    b = X.T @ W @ y
    return np.linalg.solve(A, b)

Kernel SHAP

これらを全部やればKernel SHAPとなる

実際のライブラリの実装は大規模データの場合にランダムサンプリングするなど高速化の工夫があるようだが、スモールデータでは一致する

def kernel_shap(x, f):
    """Kernel SHAPの簡単な実装例"""
    X, S, weights = create_data(x)
    y = f(X) # 予測モデルの出力値
    W = np.diag(weights)  # 重み行列
    phis = solve_wls(S, W, y) # NOTE: 入力の行列はX（特徴量の値）ではなくS（特徴量の出現の有無）
    return phis

phi = kernel_shap(x, f)
base_value = phi[-1]
shap_values = phi[:-1]
print("          x =", x)
print("shap_values =", shap_values)
print(" base_value =", base_value)
print("   sum(phi) =", np.sum(phi))
print("       f(x) =", f(x))

          x = [10 15]
shap_values = [ 50. 150.]
 base_value = 10.000000000000002
   sum(phi) = 210.0
       f(x) = 210

パッケージとの比較

import shap

reference = np.zeros(len(x))
explainer = shap.KernelExplainer(f, np.reshape(reference, (1, len(reference))))
shap_values = explainer.shap_values(x)
print("shap_values =", shap_values)
print("base value =", explainer.expected_value)

shap_values = [ 50. 150.]
base value = 10.0

参考文献

• Simple Kernel SHAP — SHAP latest documentation

  • 公式ドキュメントでの実装つき説明