機械学習の幅広い分野への応用が進むにつれ，機械学習がその予測の根拠などを理解できない「ブラックボックス」となることが問題視されており，機械学習の解釈性や説明性が注目されています．今回のテーマであるSHAP（SHapley Additive exPlanations）は，機械学習モデルへの特定の入力に対する予測の根拠を提示する代表的な手法の一つです．SHAPには用途に応じていくつかのアルゴリズムがありますが，その中でも今回はあらゆる機械学習モデルに適用可能（Model-Agnostic）なKernel SHAPという手法についてまとめました．

構成としては，まずKernel SHAPとは何かについての概要を述べた後に， Kernel SHAPを理解する上で必要な要素である「シャープレイ値」と「SHAP」について説明します．さいごに，Kernel SHAPについて「理論」と「実装」に分けて書いていきます，「理論」の部分は基本的にはSHAPの提案論文をベースにしていますが，論文に書かれていることだけではKernel SHAPのアルゴリズムを完全に理解することはできず，理解のためにはSHAPの提案者が開発しているライブラリslundberg/shapのコードリーディングが必要でした．そのため，「実装」の部分で実際にKernel SHAPがどのような流れで計算されているのかをコードから追い，重要だと思ったポイントを私なりにまとめました．

Kernel SHAPの概要

まずは，理論や実装などの細かい話は後回しにして，Kernel SHAPとはどういったものかを見ていく．SHAP（SHapley Additive exPlanations）は，協力ゲーム理論のシャープレイ値に基づき機械学習モデルの予測を解釈する統一的なフレームワークであり，2017年にLundbergらにより提案された．予測を解釈する方法として，SHAPではモデルの予測値に対する特徴量の寄与をその度合いとともに提示することができる．

論文の中では，いつくかのSHAPのアルゴリズム（Max SHAP，Deep SHAPなど）が提案されているが，その中の一つであるKernel SHAPは解釈を与える機械学習モデルを限定しないModel-Agnosticな手法である． Lundbergらによるライブラリ実装であるslundberg/shapを用いると，以下のように簡単に使うことができる．

import shap
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC

# Train a SVM classifier
X_train, X_test, y_train, y_test = train_test_split(*shap.datasets.iris(), test_size=0.2, random_state=1)
model = SVC(kernel='rbf', probability=True)
model.fit(X_train, y_train)

# Compute SHAP values
explainer = shap.KernelExplainer(model.predict_proba, X_train, link="logit")
shap_values = explainer.shap_values(X_test, nsamples=100)

以下のように予測に対する特徴量の寄与の度合いを可視化することもできる．

shap.initjs()
shap.force_plot(explainer.expected_value[0], shap_values[0][0,:], X_test.iloc[0,:], link="logit")

f:id:hirotsuru314:20210715141054p:plain — 図1 Kernel SHAPによる解釈（force plot）

上の例では，SVMを用いてアヤメの特徴量（がくと花弁の長さ・幅）からアヤメの品種（setosa，versicolor，virginica）を予測している．SVMのモデルがある入力データに対してsetosaである確率を約97%だと予測しており，Kernel SHAPの結果は，その予測に最も寄与しているのは，「petal length（がくの長さ）= 1.2 cm」であると解釈を与えていることになる．

SHAPのライブラリのREADMEを見ると，他のSHAPアルゴリズムやデータセットへの適用事例や豊富な可視化の例を見ることができる．

シャープレイ値

Kernel SHAPはシャープレイ値を近似するための手法であるため，Kernel SHAPを理解する前にシャープレイ値への理解が必要である．シャープレイ値とは，協力ゲーム理論において複数プレイヤーの協力によって得られた利得（報酬）を各プレイヤーの貢献度に応じて公正に分配するための手段の一つである．

具体的な計算方法：特性関数形ゲーム

例として，3人のプレイヤー（1，2，3）が協力してゲームに挑戦し，利得として以下の賞金が得られるとする．

f:id:hirotsuru314:20210715145855p:plain:w400 — 表1 協力ゲームの賞金表の例

例えば，プレイヤー1だけでゲームをプレイした場合は，4万円，プレイヤー1と2が二人で協力してプレイした場合は16万円となる．ここで，協力ゲームの一つの表現方法である特性関数形ゲームを考える．

特性関数形ゲームは，プレイヤーの集合 $N= \left\{1,2,...,n\right\}$ と特性関数 $v$ の組 $(N, v)$ で表される．ここで，特性関数とは，プレイヤーの部分集合 $S \in N$ （提携とも呼ばれる）に対して，そのプレイヤーが協力したときに得られる利得を与える関数である．例えば，上の例では， $v(\left\{1,2\right\})=16$ となる．シャープレイ値は，全プレイヤーの協力により得られる利得，つまり $v(\left\{1,2,3\right\})=60$ を各プレイヤーの貢献度に応じて分配する解の一つである．

ここで，限界貢献度と呼ばれるものを導入する．限界貢献度とは，ある部分集合 $S$ に対してプレイヤー $i$ が参加したときの利得の増加分 $v(S\cup\left\{i\right\})-v(S)$ である．プレイヤー $i$ の限界貢献度は，参加する部分集合 $S$ に依存する．すなわち，限界貢献度はプレイヤーが参加した順番に依存する．そのため，以下のように，全てのプレイヤーの参加順に対して，各プレイヤーの限界貢献度を計算する．

f:id:hirotsuru314:20210715161131p:plain:w550 — 表2 各プレイヤーの限界貢献度

例えば，プレイヤーの参加順「1→2→3」のときのプレイヤー3の限界貢献度は， $v(\left\{1,2,3\right\})-v(\left\{1,2\right\})=60-16=44$ のように計算できる．

シャープレイ値とは，全ての部分集合に対するプレイヤーの限界貢献度の期待値であるため，各プレイヤーのシャープレイ値は以下のようになる．

プレイヤー1： $(4+4+10+30+12+30)/6=15$

プレイヤー2： $(12+38+6+6+38+20)/6=20$

プレイヤー3： $(44+18+44+24+10+10)/6=25$

このように計算されたシャープレイ値を全プレイヤーで総和を取ると，全プレイヤーの協力により得られる利得に一致しており，シャープレイ値がそのまま利得の分配の解になっていることがわかる．

$15+20+25=60=v(\left\{1,2,3\right\})=v(N)$

ここまででシャープレイ値の計算の仕方を具体的に見てきたので，最後にシャープレイ値の数式表現を見ておく．

ゲーム $(N, v)$ において，全ての部分集合（提携）の形成が等しい確率で起こる場合，各プレイヤーの限界貢献度の期待値をシャープレイ値と呼び，プレイヤー $i$ のシャープレイ値は以下のように表される．

$\displaystyle \phi_i=\sum_{S\subseteq N\backslash\{i\}}\frac{|S|!(n-|S|-1)!}{n!}(v(S\cup\{i\})-v(S)) \tag{1}$

ここで， $|S|$ は部分集合 $S$ に含まれるプレイヤー数である．式の意味を理解するために，式の要素を分解して見ていく．式の後半の $v(S\cup\left\{i\right\})-v(S)$ は，限界貢献度そのものであり， $S$ にプレイヤー $i$ が参加した場合の利得の増加分を表している．分数で表されている部分の分母は，全プレイヤーの順列の総数 $n!$ である．分子は，プレイヤー $i$ が加わる前の部分集合 $S$ の順列の数 $|S|!$ とプレイヤー $i$ が加わった後に， $S$ と $i$ 以外の残りのプレイヤーが参加する順列の数 $(n-|S|-1)!$ をかけ合わせたものである．そしてこれらの計算をプレイヤー $i$ を含まない全ての部分集合 $S$ に対して足し合わせたものがシャープレイ値となる．

機械学習への応用

シャープレイ値は近年，機械学習モデルの予測結果に対する各特徴量の重要性の指標として，広く注目されている．協力ゲームと機械学習という一見繋がりのないものがどのように関係しているかを以下の図に示している．

f:id:hirotsuru314:20210715215701p:plain — 図2 協力ゲームと機械学習の対比

協力ゲームでは，複数のプレイヤーが協力してゲームに取り組んだ結果，特定の利得が獲得できる．そして，この得られた利得を各プレイヤーに分配することを考える．これを機械学習に置き換える場合，プレーヤーが特徴量，利得が予測値に対応する．つまり，複数の特徴量がモデルに入力され，それぞれの特徴量が互いに影響を及ぼし合った結果，モデルから特定の予測値が得られる．そして，得られた予測値を各特徴量に分配することで予測値に対する特徴量の貢献度を求めようというわけである．

シャープレイ値の計算上の課題

前述のように特徴量をゲームのプレイヤーに見立てて，予測値に対する貢献度を計算するというのは直感的にも理解できると思うが，実際にシャープレイ値を機械学習モデルの予測の解釈に適用する場合に考えなければならない課題がある．

1．計算コスト

1つ目は，「計算コスト」の課題であり，これはシャープレイ値を計算する際の一般的な課題である．各特徴量（プレイヤー）のシャープレイ値を計算する場合，表2に示したように全ての特徴量の順列のパターンを考慮しなければならない．特徴量の数を $n$ とすると，その順列の数は $n!$ である．前述の例のように対象が3つの場合は， $3! = 6$ であるため問題とならないが，例えば， $n=10$ のとき，順列の数は約362万となる．したがって， $n$ が大きくなるに伴い，現実的な時間内でシャープレイ値を厳密に計算することは不可能になる．

2．特徴量の不在の再現

2つ目は，「特徴量の不在の再現」の課題であり，これは，シャープレイ値を機械学習モデルの予測の解釈に適用する際の固有の課題である．各特徴量のシャープレイ値を計算する場合，表1に示したような各部分集合 $S$ に対する機械学習モデルの予測値を得る必要がある．これは一般的な機械学習モデルでは難しい．なぜなら，ほとんどの機械学習モデルは，全ての特徴量が入力された場合の予測値を得ることはできるが，ある特徴量が存在しない場合の予測値を得ることはできない．そのため，このような部分集合，つまり一部の特徴量が不在である場合の機械学習モデルの予測値をどのように得るかを考えなければならない．これは，シャープレイ値を機械学習モデルの予測の解釈に適用する場合に，特性関数形ゲーム $(N, v)$ の特性関数 $v$ をどのように設計するかという問題と捉えることもできる．

後ほど，これらの問題をKernel SHAPがどのように解決しているのかを説明する．

SHAP (SHapley Additive exPlanations)

2017年にLundbergらによりSHAPという機械学習の予測を解釈するためのフレームワークが提案された．SHAPでは，モデルの予測に関する説明をモデルそのものとして見るという視点を導入し，これを説明モデルと呼ぶ．これにより，当時の先行手法であるLIMEやDeepLiftなどの手法で用いられている説明モデルについて共通点を見出し，以下で定義する「Additive Feature Attribution Methods」のクラスとして統一的に表現できることを示した．また，そのクラスにおいてある条件を満たした場合，説明モデルの解はシャープレイ値のみに定まることが理論的に証明されている．

Additive Feature Attribution Methods

解釈したい予測モデルを $f$ ，説明モデルを $g$ として，単一の入力 $x$ に対する予測値 $f(x)$ を解釈することに焦点を当てる．説明モデルは入力として，写像関数 $h_x$ を介して元の入力 $x$ に写像する簡略化した入力（論文中ではsimplified inputs） $x'$ を用いることが多い．

Additive Feature Attribution Methodsは，簡略化した入力をバイナリ変数 $z' \in \left\{0,1\right\}^M$ （Mは入力の特徴量数）で表現し，説明モデル $g$ を以下のように定義する．

$\displaystyle g(z')=\phi_0+\sum_{i=1}^M\phi_iz_i' \tag{2}$

つまり，説明モデルをバイナリ変数の線形モデルとして定義している．この説明モデルを $g(z′)≒f(h_x(z′))$ となるようにして元のモデルを近似する． $\phi_i$ は各特徴量の予測値に対する効果を表しており，すべての特徴量の効果を合計することで，元のモデルの予測値 $f(x)$ を近似する．

Additive Feature Attribution Methodsの性質

Additive Feature Attributionのクラスは，以下に示す３つの望ましい性質を満たす唯一の解が存在する．

1．Local accuracy

簡略化された入力 $x'$ （ $x＝h_x(x')$ ）に対する説明モデルの出力が，元のモデル $f$ の予測値 $f(x)$ と一致すること．別の表現をすると，全ての特徴量が存在しないの場合のモデルの出力 $\phi_0=f(h_x(\mathbf{0}))$ に各特徴量の貢献度 $\phi_i$ の合計値を加えたものがモデル $f$ の予測値 $f(x)$ と一致している．

$\displaystyle f(x)=g(x')=\phi_0+\sum_{i=1}^M\phi_ix_i'$

2．Missingness

存在しない特徴量は予測値に影響を与えない．

$\displaystyle x_i'=0$ のとき $\displaystyle \phi_i=0$

3．Consistency

$f_x(z')=f(h_x(z'))$ とおき， $z'\backslash i$ を $z_i'=0$ を表すものとする．もし，二つのモデル $f$ と $f'$ が，

$\displaystyle f_x'(z')-f_x'(z'\backslash i)\geq f_x(z')-f_x(z'\backslash i)$

を全ての入力 $z' \in \left\{0,1\right\}^M$ に対して満たす場合，下の式が成り立つ．

$\displaystyle \phi_i(f',x) \geq \phi_i(f,x)$

この性質は少しわかりづらいので，言葉で説明すると，モデルの変更（ここでは， $f$ → $f'$ ）により，任意の特徴量がモデルの予測値に与える影響が増加または同じになった場合，その特徴量の貢献度は減少しない（増加または同じになる）ということである．

論文中では，Additive Feature Attribution Methodsの定義である式(2)に従う説明モデル $g$ が上記の三つの性質を満たすための $\phi_i$ の唯一の解はシャープレイ値であることが示されている．

SHAP Values

SHAPでは，特徴量の重要性の統一的な尺度としてSHAP Values（以下，SHAP値）を導入する．SHAP値とは，元のモデルの条件付き期待関数のシャープレイ値である．下の図は，特徴量を何も知らない場合に予測されるベース値 $E[f(z)$ ]から入力 $x$ に対するモデルの予測値 $f(x)$ に至る過程を示している．

f:id:hirotsuru314:20210716155453p:plain — 図3 出典：参考文献[1]のFigure1

重要なこととして，上図の条件付き期待値をどのように計算するかという問題がある．これに対する具体的な計算方法は，後述の「Kernel SHAPの理論」の中で述べる．

また，上の図は一つの順列のケースを示している．したがって，厳密にSHAP値を計算する場合は，全ての順列を考慮して $\phi_i$ を平均化する必要がある．このことは，前述のシャープレイ値の計算上の課題の一つ目で述べたことと同様に計算コスト上の問題となる．次に説明する本記事の本題であるKernel SHAPはSHAP値を近似する効率的な方法を提供している．

Kernel SHAPの理論

定式化

Kernel SHAPではSHAP値の計算を重み付きの最小二乗線形回帰問題として定式化する．重み $\pi_{x'}$ と損失関数 $L$ は以下のように定義される．

$\displaystyle \pi_{x'}(z')=\frac{M-1}{(M\ choose\ |z'|)|z'|(M-|z'|)} \tag{3}$

$\displaystyle L(f,g,\pi_{x'})=\sum_{z'\in Z}[f(h_x(z'))-g(z')]^2 \pi_{x'}(z') \tag{4}$

ここで， $|z'|$ は $z'$ の中の非ゼロ要素の数であり，重み $\pi_{x'}$ はShapley kernel（シャープレイカーネル）と呼んでいる．式(2)の説明モデル $g$ の定義が成り立つ下で，この損失関数 $L$ を最小化する説明モデル $g$ のパラメータ $\phi_i$ を求めた場合，その解は前述のAdditive Feature Attribution Methodsの性質1〜3を満たす，すなわち解がSHAP値の近似値であることが証明されている．（証明は論文のSupplementary Material参照） $g(z')$ は線形モデルであり， $L$ は二乗誤差であるので，結果としてSHAP値は重み付き線形回帰を用いて計算することができる．

シャープレイカーネルの性質

下の図は，部分集合をcardinalityで並べた場合の重みの大きさを示している．ここで，cardinalityとは集合の濃度のことで，有限集合の場合は集合の元の数に相当する．

f:id:hirotsuru314:20210716213519p:plain:w550 — 図4 出典：参考文献[1]のFigure2 (A)

図からわかるように重みを部分集合のcardinalityで並べた場合，重みは対称になっている．また，グラフの形状から，バイナリベクトル $z'$ の非ゼロの要素が一つもしくはゼロの要素が一つの場合は重みが大きく，非ゼロもしくはゼロの要素が増えていくにつれて重みが小さくなり，非ゼロとゼロの要素が同程度のときに重みが最も小さくなることがわかる．

Kernerl SHAPは先行手法であるLIMEに強く影響を受けている手法であるため，図4では比較対象としてLIMEを挙げている．しかし，本記事ではLIMEについての説明は割愛する．

条件付き期待値の計算

SHAP値を計算する上で，条件付き期待値（図3の $E[...$ ]で表現されている部分）をどのように計算するかは，特徴量の不在をどのように再現するかを考える問題が含まれているので重要なポイントである． SHAPの論文中でも条件付き期待値の計算についての説明があるが（論文中の式(9)〜(12)）実際の計算過程が少しイメージしづらく感じた．なぜなら，実際に条件付き期待値を計算する上で重要となる「バックグラウンドデータセット」についての考えが明示されていないからである．

ライブラリslundberg/shapのKernelExplainerの実装のコメントの中で以下のような記述がある．

To determine the impact of a feature, that feature is set to "missing" and the change in the model output is observed. Since most models aren't designed to handle arbitrary missing data at test time, we simulate "missing" by replacing the feature with the values it takes in the background dataset.

（翻訳）ある特徴量の影響を調べるために，その特徴量を「missing（不在）」に設定し，モデル出力の変化を観察する．ほとんどのモデルは，テスト時に任意のmissingデータを処理するようには設計されていないため，バックグラウンドデータセットでその特徴が取る値に置き換えることで「missing」を再現する．

これを数式で表現すると，解釈したいデータを $x$ ，バックグラウンドデータセットを $X$ ， $\bar{S}$ を部分集合 $S$ に存在しない特徴量の集合（ $N\backslash\ S$ ）として以下のように表せる．

$\displaystyle E[f(x_S,X_{\bar{S}})] \tag{5}$

このようにすることで，バックグラウンドデータセットのサイズを適切に選択するとSHAP値の条件付き期待値の計算量を削減できる．一方，これは特徴量間の独立を想定していることに注意しなければならない．

シャープレイ値の計算上の課題の解決

理論の話の最後に，先に述べた「シャープレイ値の計算上の課題」をKernel SHAPがどのように解決しているのかを述べる．

1．計算コスト

特徴量数 $n$ が大きくなった際に考慮すべき場合の数が膨大になる問題について，Kernel SHAPではモンテカルロ近似を採用している．重み付き最小二乗問題を解くために必要なデータ点を一定の個数サンプリングすることで計算量を抑えている．必要なサンプル数やサンプリングの優先順位などの細かい話は実装のところで述べる．

2．特徴量の不在の再現

特徴量が不在であることをどのように再現するかについては，「条件付き期待値の計算」で述べた通りである． Kernel SHAPでは不在の特徴量の値をバックグラウンドデータセットの値と入れ替えたのちに期待値を計算している．

Kernel SHAPの実装

SHAPの提案論文の著者によるライブラリslundberg/shapを用いると，以下のように非常に簡単にKernel SHAPを用いたSHAP値の計算が実行できる．（model.predictは任意の機械学習モデルの予測値を返す関数，X_train，X_testはそれぞれ学習データ，テストデータが用意されているとする．）

explainer = shap.KernelExplainer(model.predict, X_train)
shap_values = explainer.shap_values(X_test)

このコードの内部で何が起きているかを知り，Kernel SHAPのアルゴリズムについてより理解を深めるためにKernel SHAPの実装の実装を追っていく．

Kernel SHAPの実行の流れは，以下の三つのステップで構成される．

部分集合（バイナリベクトル）のサンプリング
バイナリベクトルを元の特徴空間に変換し，モデルの予測値を得る
重み付き最小二乗法を解く

以下，それぞれのステップについて詳しく見ていく．（所々コードを引用しているがコードは2021年7月時点のもの）

1. 部分集合（バイナリベクトル）のサンプリング

特徴量の集合 $N= \left\{1,2,...,n\right\}$ としたとき，部分集合 $S \in N$ をサンプリングする．以降，部分集合をバイナリベクトル $z'$ で表現する．例えば， $n=5$ で，部分集合 $S=\left\{1,3\right\}$ のとき， $z' = [1, 0, 1, 0, 0$ ]のように部分集合内に存在する特徴量を1，存在しない特徴量を0としたバイナリベクトルを考える．本来この部分集合は，バイナリベクトルの定義から明らかなように $2^n$ 個の組み合わせの数がある．しかし，これでは $n$ が大きくなった場合に計算が現実的ではなくなるため，Kernel SHAPでは一定のサンプル数のサンプリングを行う．

サンプル数

このサンプル数をどのように決定するかという問題がある．ライブラリの実装では，shap_values関数にnsamplesという引数があり，この引数でユーザが任意のサンプル数を指定できるようになっている．しかし，適切なサンプル数を指定することは困難であるため，ライブラリでは推奨値を持っている．引数nsamplesを指定しない場合，基本的には以下のように「 $2 * n + 2048$ 」でサンプル数が決定される．

# self.Mは特徴量数
if self.nsamples == "auto":
    self.nsamples = 2 * self.M + 2**11

本来は，特徴量数 $n$ に対して， $2^n$ で増えていく部分集合の数を推奨値では $n$ の線形に置き換えているので，特徴量数 $n$ が大きい場合は，推奨値ではサンプル数が小さすぎるかもしれない．これは計算時間と精度のトレードオフになる部分なので，自分たちの計算時間に対する要求も含めて，場合によっては独自に設定した方が良いと思われる．

サンプリングの優先順位

ライブラリの実装ではバイナリベクトルを完全にランダムでサンプリングしているわけではなく，サンプリングの優先順位がある．実装では，以下のように特徴量数の半分（の切り上げ）をfor文で回していくようになっているが，これはバイナリベクトル中の非ゼロの要素数（subset_size）を決めている．

num_subset_sizes = np.int(np.ceil((self.M - 1) / 2.0))
for subset_size in range(1, num_subset_sizes + 1):
    ・・・

そして，要素を非ゼロとする位置の各組み合わせに対してサンプルを追加（addsample）していくのだが，ここで重要なのが，同時にバイナリベクトル中の0と1の要素の反転したベクトルも追加している．

# subset_size=非ゼロ要素数
for inds in itertools.combinations(group_inds, subset_size):
    mask[:] = 0.0 # mask=バイナリベクトル
    mask[np.array(inds, dtype='int64')] = 1.0
    self.addsample(instance.x, mask, w)
    if subset_size <= num_paired_subset_sizes:
        mask[:] = np.abs(mask - 1) # 0と1の要素の反転
        self.addsample(instance.x, mask, w)

この意図は，図4を見るとわかる．図4の部分集合のcardinalityが小さいケースと大きいケースから真ん中に向かってサンプリングをしている．つまり，重み付き線形回帰において重みが大きいサンプルから優先にサンプリングしていることがわかる．

2. バイナリベクトルを元の特徴空間に変換し，モデルの予測値を得る

これは「Kernel SHAPの理論」の「条件付き期待値の計算」で述べた部分である．バイナリベクトルの要素が1（特徴量が存在する）である場合，解釈したいデータの値をそのまま入れるだけで良い．バイナリベクトルの要素が0（特徴量が存在しない）場合は前述の式(5)ようにバックグラウンドデータセットの値で置き換えて期待値をとる．

バックグラウンドデータセットは，KernelExplainerクラスのインスタンス作成時に渡しているX_trainで設定している．

explainer = shap.KernelExplainer(model.predict, X_train)

バックグラウンドデータセットは独自に定義することができるが，学習データの一部をそのまま渡す，もしくはクラスタリングしてクラスタの代表値を渡すといったやり方がよく用いられる．

バックグラウンドデータセットとして，複数のデータを渡すのではなく，データを一つのみ渡すこともできる．例えば，学習データの平均値を計算して，一つのデータのみをバックグラウンドデータセットとして設定しても良い．

計算のボトルネック

Kernel SHAPの計算の中で，式(5)の期待値の中のモデル $f$ による予測値の計算がボトルネックになる．なぜなら，モデルの予測は「サンプル数×バックグラウンドデータセットのサイズ」個のデータを予測しなければならないためである．例えば，特徴量数が10の場合，そのサンプル数の推奨値は「2*10+2048=2068」，またバックグラウンドデータセットのサイズを100とすると，206,800個のデータに対してモデルの予測を行わなければならない．この計算は，当然ながら用いるモデルの計算時間に依存するが，ほとんどの場合でこの計算は後述の重み付き最小二乗法を解くより時間がかかる．

ちなみに，ライブラリ実装では，バックグラウンドデータセットのサイズが100を超えると以下のように警告を発するようになっている．

# warn users about large background data sets
if len(self.data.weights) > 100:
    log.warning("Using " + str(len(self.data.weights)) + " background data samples could cause " +
                "slower run times. Consider using shap.sample(data, K) or shap.kmeans(data, K) to " +
                "summarize the background as K samples.")

特徴量の不在を再現するための良い参照値が定義できるのであれば，バックグラウンドデータセットのサイズはなるべく小さくする方が計算時間的に望ましい．

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

ここまでで，サンプリングしたバイナリベクトルと，そのバイナリベクトルごとのモデルの予測値が得られているため，あとは式(4)の損失関数を最小化するように重み付きの最小二乗線形回帰問題を解くだけである．サンプルであるバイナリベクトルごとに与えられる重みは，式(3)から特徴量数とバイナリベクトル中の非ゼロ要素数で計算ができる．

ライブラリ実装では，shap_values関数にl1_reg引数を持っている．これは，特徴選択のためのL1正則化のための引数で，適切に設定した場合，通常の線形回帰ではなく，sklearnのLassoLarsICやLassoを用いたラッソ回帰が適用されるようになっている．

デフォルト（l1_reg="auto"）では，サンプル数が，下のコードのように計算されるmax_samplesの20％以下の場合に正則化にAICが使用され，それ以外の場合は正則化は行われない．

self.max_samples = 2 ** 30
if self.M <= 30:
    self.max_samples = 2 ** self.M - 2
    if self.nsamples > self.max_samples:
        self.nsamples = self.max_samples

しかし，ライブラリのコメントの中で，デフォルトのときの挙動を変更する予定があるとの内容があったので，今後のバージョンアップの際は注意が必要である．

THE BEHAVIOR OF "auto" WILL CHANGE in a future version to be based on num_features instead of AIC.

最小二乗法により求められた回帰モデルの係数はそのままSHAP値となっているため，計算はこれで終わりである．通常，この後にSHAP値の可視化を行うことになるが，それはライブラリのREADMEを見ると色々な可視化ができることがわかる．

さいごに

今回は，協力ゲーム理論のシャープレイ値に基づきあらゆる機械学習モデルの予測を解釈できるKernel SHAPについて，その理論と実装をまとめた． SHAPは2017年に提案されて以降も様々なアルゴリズムの改善や派生アルゴリズムが提案されている．特に，2017年の提案論文のラストオーサーであるS. Lee氏らにより，Kernel SHAPの改善が2020年に以下の論文で提案されており，非常に興味深いので，今後はこの辺りも追っていきたいと思う．

arxiv.org