Fire Engine

化学の修士号→消防士→ITエンジニア(2016年11月〜)

機械学習の予測を解釈するKernel SHAPの高速性と拡張性の向上を目指したライブラリを開発した

先日,協力ゲーム理論のシャープレイ値に基づき機械学習モデルの予測を解釈するKernel SHAPという手法の理論と既存のライブラリの実装についてのブログを書いた.

blog.tsurubee.tech

既存のSHAPライブラリであるslundberg/shap(以下,単にSHAPライブラリと呼ぶ)は,SHAPの提案論文*1のファーストオーサーにより開発されており,多くのSHAPのアルゴリズムの実装や可視化の機能が非常に充実している素晴らしいライブラリである.

しかし,私が自身の研究の中でSHAPライブラリの中のKernel SHAPを使っている際に,計算速度と拡張のしやすさの観点で改善したいポイントがいくつか出てきた.今回は,まだ絶賛開発中であるが,Kernel SHAPの高速性と拡張性の向上を目指したShapPackというライブラリのプロトタイプが完成したので,それについて紹介する.

目次

ShapPackの概要

github.com

ShapPackではSHAPライブラリのKernel SHAPと比較して,以下の三つの新たな機能を実装している.

  1. マルチプロセスで並列処理できる機能
  2. 特性関数を独自に実装して組み込める機能
  3. SHAP値を計算しない特徴量を指定できる機能

1と3は高速性に関わる機能で,2は拡張性に関わる機能である.
わかりやすい結果として,1の並列処理の機能により8コアのサーバでSHAP値の計算を実行すると,以下のように計算時間が短縮できる.

SHAPライブラリ ShapPack
5.54 s 0.684 s

上の結果の実験条件や,他の二つの機能の詳細については後述する.
下図は,SHAPライブラリとShapPackの実装の違いの概要をまとめたものである.

f:id:hirotsuru314:20210720165023p:plain
図1 SHAPライブラリとShapPackの実装の違いの概要図

ShapPackでは,サンプリングされた部分集合に対する特性関数の計算がボトルネックになることに着目し,部分集合を分割してマルチプロセスで並列処理することで計算速度を改善している.また,SHAPライブラリでは特性関数としてライブラリに実装済みのものしか利用できないことに着目し,利用者が外から独自に実装した特性関数を組み込める仕組みをとることで拡張性を改善している.

ShapPackの使い方

データとモデルの準備

データとして,scikit-learnに付属しているボストンの住宅価格のデータセットを用いた.データは,特徴量の尺度を揃えるために標準化している.
モデルは,動径基底関数(RBF)カーネルのサポートベクトル回帰(SVR)を用いた.

import numpy as np
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVR

# Fix seed to reproduce results
SEED = 123
np.random.seed(SEED)

# Prepare dataset
boston = load_boston()
X_train, X_test, y_train, y_test = train_test_split(boston["data"], boston["target"], test_size=0.2, random_state=SEED)
scaler = StandardScaler()
X_train_std = scaler.fit_transform(X_train)
X_test_std = scaler.transform(X_test)

# Prepare model
model = SVR(kernel="rbf")
model.fit(X_train_std, y_train)

SHAP値の計算

ShapPackを用いてSHAP値を計算するためのコードは以下の通りである.

import shappack
i = 2
explainer = shappack.KernelExplainer(model.predict, X_train_std[:100])
shap_value = explainer.shap_values(X_test_std[i])

ShapPackでは,既に広く使われているSHAPライブラリとなるべく同じ使い方ができるようにしており,実際に上のコードは「shappack」を「shap」と置き換えるとそのまま実行できる.
SHAPライブラリとの重要な違いとして,ShapPackのshap_values関数は,前述の三つの機能に紐づく新たな三つの引数(n_workerscharacteristic_funcskip_features)を渡せるようになっており,これらの引数が高速性や拡張性に寄与する.
以下,それぞれの引数の使い方を説明する.

1. n_workers:プロセス数を指定
shap_value = explainer.shap_values(X_test_std[i], n_workers=-1)

n_workersには,SHAP値の計算に使用するプロセス数を指定する.指定しない場合は,n_workers=1となる. n_workers=-1は,実行するサーバのコア数をn_workersに設定することを意味する. この引数により,プログラムをマルチコアのサーバーで実行する場合は,計算時間の短縮が期待できる.

2. characteristic_func:特性関数を指定
def my_characteristic_func(instance, subsets, model, data):
    # own implemented characteristic function

shap_value = explainer.shap_values(X_test_std[i], characteristic_func=my_characteristic_func)

characteristic_funcには,独自に実装した特性関数を渡すことができる.指定しない場合は,SHAPライブラリのKernel SHAPと同等のものが実行される.

3. skip_features:計算をスキップする特徴量を指定
explainer = shappack.KernelExplainer(model.predict, X_train_std[:100], feature_names=boston.feature_names)
shap_value = explainer.shap_values(X_test_std[i], skip_features=["PTRATIO", "TAX"])

skip_featureには,SHAP値の計算をスキップしたい特徴量を指定する.特徴量は,特徴量名またはインデックス番号で指定できる. なお,skip_featuresを特徴量名で指定する場合は,KernelExplainerクラスのfeature_names引数に特徴量名のリストを渡す必要がある.

可視化

現時点では,ShapPackは独自に可視化の仕組みを持っていないため,SHAP値の可視化のためにはSHAPライブラリを使う必要がある.

import shap
shap.initjs()
shap.force_plot(explainer.base_val[0], shap_value, X_test[i], boston.feature_names)

f:id:hirotsuru314:20210720171941p:plain
図2 Kernel SHAPに計算された予測値への貢献度の可視化

ちなみに,上の結果は,あるデータの住宅価格を平均より高く予測しており,その予測に対してRM(1戸当たりの平均部屋数)やLSTAT(低所得者人口の割合)が強く貢献していることを示している.

ShapPackの詳細

ここでは,新たに追加した三つの機能についての詳細をそれぞれ説明する.

1. マルチプロセスで並列処理できる機能

Kernel SHAPを利用する際の大きな課題の一つとして,計算時間が挙げられる.ShapPackではSHAPライブラリのボトルネックとなる箇所をマルチプロセスで並列処理できるようにすることで,計算速度の向上を目指している.
ここでは,まずSHAPライブラリとShapPackの計算時間の比較評価の結果を示す. 次に,SHAPライブラリのKernel SHAPのボトルネックと,ボトルネックを解消するための並列計算の実装について述べる.

計算時間の評価

下のコードのようにJupyter NotebookでSHAP値の計算時間の10回平均値を測定した.

%%timeit -r 10
shap_value = explainer.shap_values(X_test_std[i])

データやモデルは「ShapPackの使い方」で示した例と同じものを使用し,サーバはCPU8コア,メモリ24GB,バックグラウンドデータセットのサイズは100を採用した. 評価結果を下の表1に示す.

f:id:hirotsuru314:20210721092848p:plain:w500
表1 SHAP値の計算時間(単位は秒)

表1から,SHAPライブラリの計算時間が5.54秒であるのに対し,n_workers=8で設定したShapPackの計算時間は0.684秒であり,本実験条件では約8倍の計算時間の短縮が達成できた.これは,単純にSHAPライブラリだとマルチコアのサーバでも1コアしか計算に使わない一方で,ShapPackでマルチコアをフルに使って計算を実行するためである.
ShapPackでマルチコアを使わないn_workers=1でもSHAPライブラリよりも早い結果を示している.これはSHAPライブラリを参考にShapPackの実装を進める中でいくつかの細かい実装の工夫をしたためだと思われるが,その詳細はここでは割愛する.

SHAPライブラリのボトルネック

Jupyter Notebookを用いる場合,以下のように%%prunをセルの先頭に差し込むだけで比較的簡単にコードのプロファイリングができる.

%%prun
shap_value = explainer.shap_values(X_test_std[i])

SHAPライブラリのSHAP値の計算に対する出力結果を下に示す.

85988 function calls (82014 primitive calls) in 5.572 seconds

   Ordered by: internal time

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
     2    4.719    2.360    4.719    2.360 {sklearn.svm._libsvm.predict}
     1    0.731    0.731    5.455    5.455 _kernel.py:503(run)
  2074    0.022    0.000    0.024    0.000 _kernel.py:477(addsample)
     1    0.012    0.012    5.572    5.572 _kernel.py:204(explain)
    ・・・

この結果から,全体の計算時間(5.572秒)に対してSVMのpredict関数(4.719秒)が支配的であることがわかる.実際に該当する箇所は下のコードの最後のself.model.f(data) である.

# https://github.com/slundberg/shap/blob/master/shap/explainers/_kernel.py
def run(self):
    num_to_run = self.nsamplesAdded * self.N - self.nsamplesRun * self.N
    data = self.synth_data[self.nsamplesRun*self.N:self.nsamplesAdded*self.N,:]
    if self.keep_index:
        index = self.synth_data_index[self.nsamplesRun*self.N:self.nsamplesAdded*self.N]
        index = pd.DataFrame(index, columns=[self.data.index_name])
        data = pd.DataFrame(data, columns=self.data.group_names)
        data = pd.concat([index, data], axis=1).set_index(self.data.index_name)
        if self.keep_index_ordered:
            data = data.sort_index()
    modelOut = self.model.f(data) #Bottleneck!
    ・・・

なぜここの処理が時間がかかるかというと,モデルに渡ってくるdataのサイズが大きいからである.例えば,ボストンの住宅価格のデータセットの場合,特徴量数が13であり,部分集合のサンプリング数のライブラリの推奨値が「2*13+2048=2074」となる.また,特性関数の計算に用いるバックグラウンドデータセット(KernelExplainerクラスの引数の二つ目)のサイズを100とすると,207,400個のデータに対してモデルのpredict関数を実行しなければならない.この計算がKernel SHAPの計算時間のボトルネックになる.

並列計算の実装

サイズの大きなデータに対するモデルの予測を早く計算する方法として,ShapPackではデータを分割してマルチプロセスでモデルの予測を計算するように実装している(図1).もちろん,用いるモデル自体がマルチプロセスでの並列計算をサポートしていれば,そちらの機能を使えばよいが,全てのモデルがマルチプロセス処理に対応しているわけではないので,ShapPack側でマルチプロセス処理を実装し,引数で並列プロセス数を指定できるようにした.

2. 特性関数を独自に実装して組み込める機能

特性関数の実装例

以下にShapPackに組み込む特性関数を独自に実装した例を示す.

def my_characteristic_func(instance, subsets, model, data):
    n_subsets = subsets.shape[0]
    n_data = data.shape[0]
    synth_data = np.tile(data, (n_subsets, 1))
    for i, subset in enumerate(subsets):
        offset = i * n_data
        features_idx = np.where(subset == 1.0)[0]
        synth_data[offset : offset + n_data, features_idx] = instance[:, features_idx][0]
    model_preds = model(synth_data)
    ey = np.zeros(n_subsets)
    for i in range(n_subsets):
        ey[i] = np.min(model_preds[i * n_data : i * n_data + n_data])
    return ey

shap_value = explainer.shap_values(X_test_std[i], characteristic_func=my_characteristic_func)

この例は,SHAPライブラリのKernel SHAPの特性関数 E[f(x_S,X_{\bar{S}})]の期待値 Eを最小値計算 minに置き換えた例である. コードの詳細は述べないが,上のように十数行でKernel SHAPのオリジナルの特性関数を少し拡張したものを実装できる.
特性関数は「特徴量の不在をどのように再現するか」という問題が含まれており,SHAP値の計算の上では重要なポイントである.特性関数を独自に設計するというのは簡単ではないが,用途やデータに合わせてを特性関数を設計したいというケースもある(事例は後述).独自に特性関数を設計して適用したい場合,SHAPライブラリでは既存のソースコードに直接変更を加えなければならないが,ShapPackでは関数を実装して引数として渡すだけで,独自に設計した特性関数を適用できる.

特性関数を独自に設計した事例

特定の用途に適した特性関数を独自に設計している事例として,異常検知の解釈のために特性関数を設計した事例がある. 例えば,解釈したいデータ点(以下,インスタンス)に応じて適応的に参照値を選択するために,インスタンスと学習データとのinfluence weightingを用いている例がある*2.これは,ShapPackに組み込む特性関数の中で実装可能である.

また,インスタンスの近傍での異常スコアの最小化問題を解くことで特徴量の不在を近似する特性関数が提案されている*3. この手法では,インスタンスと部分集合の両方に適応的に参照値を決定することになるが,ShapPackの特性関数では,上の特性関数の実装例のinstancesubsetsでどちらも関数に渡たされる仕様になっているため,こちらもShapPackに組み込む特性関数の中で実装可能である.

3. SHAP値を計算しない特徴量を指定できる機能

SHAPライブラリでは,機械学習モデルの入力となる特徴量すべてに対してSHAP値を計算するが,ShapPackでは指定した特徴量のSHAP値の計算をスキップすることができる.これは例えば,人間が経験的に予測にほとんど影響を与えない特徴量を事前知識として持っている場合や,ある特徴量の影響がないと仮定した上での他の特徴量の貢献度を知りたい場合などを想定している.また,この機能により,指定する特徴量がある場合は計算が必要なSHAP値を減らすことができ,計算速度の向上にも貢献する.

実装としては,指定した特徴量は部分集合に必ず含まれているとし,その存在の有無の効果を測らないようにしている.例えば,以下のように特徴量「PTRATIO」と「TAX」を指定して実行した場合の可視化の結果は,図2の結果から「PTRATIO」と「TAX」が除外されていることがわかる.

explainer = shappack.KernelExplainer(model.predict, X_train_std[:100], feature_names=boston.feature_names)
skip_features=["PTRATIO", "TAX"]
shap_value = explainer.shap_values(X_test_std[i], skip_features=skip_features, n_workers=-1)
feature_names = np.delete(boston.feature_names, explainer.skip_idx)
x_test = np.delete(X_test[i], explainer.skip_idx)
shap.force_plot(explainer.base_val[0], shap_value, x_test, feature_names)

f:id:hirotsuru314:20210721091246p:plain
図3 PTRATIOとTAXの計算を除外した場合の貢献度の可視化

今後の展望

まず,現在のShapPackの大きな制約としてシャープレイ値の推定方法がKernel SHAPに限定されることが挙げられる. 一方,SHAPライブラリはKernel SHAP以外にDeep SHAPやTree SHAPなど豊富なアルゴリズムの実装がある. そのため,ShapPackでも高速性や拡張性を意識しつつKernel SHAP以外のアルゴリズムを実装していきたい.

次に,シャープレイ値を用いた機械学習の解釈の興味深い発展として,特徴量間の因果関係を導入したシャープレイ値がいくつか提案されている. 例えば,Shapley Flow*4や Causal Shapley Values*5などがある.これらの手法について調査してShapPackに実装できないか検討していきたい.

参考文献

協力ゲーム理論のシャープレイ値に基づき機械学習モデルの予測を解釈するKernel SHAPの理論と実装のまとめ

機械学習の幅広い分野への応用が進むにつれ,機械学習がその予測の根拠などを理解できない「ブラックボックス」となることが問題視されており,機械学習の解釈性や説明性が注目されています.今回のテーマである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のクラスは,以下に示す3つの望ましい性質を満たす唯一の解が存在する.

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_trainX_testはそれぞれ学習データ,テストデータが用意されているとする.)

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

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

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

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

以下,それぞれのステップについて詳しく見ていく.(所々コードを引用しているがコードは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のLassoLarsICLassoを用いたラッソ回帰が適用されるようになっている.

デフォルト(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

参考文献

[1] A Unified Approach to Interpreting Model Predictions:SHAPの提案論文

[2] 協力ゲーム理論 - 株式会社 勁草書房:協力ゲーム理論全般についてまとめられた良著

機械学習モデルの局所的な解釈に着目したシステムにおける異常の原因診断手法の構想

著者 鶴田 博文, 坪内 佑樹
所属 さくらインターネット株式会社 さくらインターネット研究所
研究会 第8回WebSystemArchitecture研究会

1. はじめに

インターネットを介して利用するシステムの大規模化に伴い,システムの構成要素数の増大や,構成要素間の関係性の複雑化が進んでいる. そのため,システムの性能に異常が発生したときに,システムの状態を示す指標であるメトリックをシステム管理者が網羅的に目視することや,メトリック間の関係性を把握することができず,システムの異常原因を特定することが難しくなっている.

この問題を解決するために,深層学習などの機械学習モデルを用いて,システムの異常の原因を診断する手法が提案されている[1,2]. これらの手法は,システム管理者が異常の根本原因を絞り込むために活用することが期待できる. しかし,原因診断を行うためには,事前に機械学習モデルの学習が必要であり,それに加えて定期的なモデルの更新も必要であるため,それらに伴う課題が存在する. 課題として,まずモデルの学習と更新に伴う計算コストがかかることが挙げられる. 次に,モデルの入力となる分析対象のメトリックを事前に指定する必要があるため,より原因に近いメトリックが診断結果から除外される可能性がある.

事前にモデルの学習が不要であり,異常発生を起点にその原因を診断できる手法として,統計的因果探索を用いた手法が提案されている[3,4]. これらの手法は,因果グラフにより異常の伝搬経路を特定できることに加え,異常発生時のみしか計算コストを要しないことが利点として挙げられる. しかし,既存手法であるMicroscope[3]やAutoMAP[4]は,分析対象のメトリックを事前に指定しておく必要があることが課題である.

本発表では,事前にモデルの学習や対象メトリックの指定を必要とせず,機械学習モデルの局所的な解釈手法であるSHAP[5]を用いてシステムの異常の原因を診断する手法の検討を行う.

2. 機械学習モデルの局所的な解釈

近年,深層学習をはじめとする機械学習モデルは様々な分野へ応用が進んでいる. 一方,機械学習モデルはその計算過程が複雑であるため,人間が動作を理解することができないブラックボックスとなることが問題視されており,機械学習モデルの解釈性についての研究が注目されている[6].

局所的な解釈は,特定の入力に対するモデルの予測や判断の根拠を提示することであり,代表的な手法としてLIME[7]やSHAP[5]が挙げられる. これらの手法は,予測や判断の根拠となった特徴量を提示する手法である. 例えば,病気になるリスクを予測する機械学習モデルがあったとして,ある人物がリスクが高いと予測されたとする. LIMEやSHAPではその根拠となる特徴量(例えば,体重,血圧など)を予測への寄与の度合いとともに提示する.

局所的な解釈手法は,特に画像認識の分野で数多くの研究が報告されているが,異常の原因診断においてもその有用性が示されつつある[8-10]. 例えば,SHAP等を用いてPCA[8]やオートエンコーダ[9],変分オートエンコーダ[10]などによる異常検知の結果の解釈が,他の手法と比較して,原因の特定精度が高い,もしくは人間の直感に近い解釈を与えるなどの報告がされている.

3. 提案する原因診断手法

「1. はじめに」で述べた通り,既存の深層学習などの機械学習モデルを用いた原因診断手法[1,2]では,原因診断のために用いるモデルを事前に学習し,定期的に更新する必要があるため,そのための計算コストがかかることや,分析対象となるメトリックをシステム管理者が事前に指定しておく必要があることが課題である. 一方,統計的因果探索を用いた原因診断手法[3,4]では,事前にモデルの学習を必要としないが,既存手法では分析対象のメトリックを事前に指定しておく必要がある. この場合,システム管理者が選定したメトリクスの中に異常の根本原因となるメトリックが含まれず,診断結果から原因メトリックが除外される可能性がある. そこで提案手法では,事前にモデルの学習や対象メトリックの指定を必要としない原因診断手法の検討を行う. 提案手法の概要図を以下に示す.

f:id:hirotsuru314:20210530220928p:plain
図1 提案手法の概要図

提案手法は,異常発生後にその原因を診断するための手法であり,異常を検知するためには,Service Level Objective(SLO)やメトリックごとに設定した閾値などを用いることを想定している. また,提案手法は,異常発生時にシステム管理者に原因診断結果を提示し,システム管理者がシステムを異常状態から復旧するための作業を支援することを目指している. そのため,診断結果を提示するまでの時間が短いことが望まれる. この要件を満たすための提案手法の工夫や実際の実行時間の評価は後述する.

提案手法は以下の三つのステップから構成される.

Step 1:メトリックのフィルタリング

提案手法では,事前に分析対象となるメトリックを指定する必要がないため,異常発生後に対象メトリックを選定できる. そのため,異常発生時にほとんど変動がないメトリックなど,その異常への関連の可能性が低いものをフィルタリングできる. メトリックのフィルタリングは,原因診断の精度の向上と後続のステップの実行時間の短縮に有効である.

異常への関連性が低いメトリックをフィルタリングする手法の一つとして,以前の我々の研究成果であるTSifter[11]の活用を検討する. TSifterは,定常性の検定および階層的クラスタリングを用いて迅速に異常に関連するメトリックを絞り込む次元削減手法である.

Step 2:モデルの学習

提案手法では,異常発生後に観測データからモデルを学習するため,高速に学習可能なモデルを用いる必要がある. 本発表では,モデルとして,主成分分析(PCA)[12]の利用を検討する.(今後,非線形のモデル等に拡張
する予定) PCAを用いた異常検知では,観測データに対する次元削減により正常部分空間を求め,テストデータと正常部分空間との距離を異常スコアとする. 提案手法では,PCAで算出される異常スコアを異常検知ではなく,検知後の原因診断に用いることを想定している.

Step 3:異常への貢献度の計算

提案手法では,異常の原因診断を行うために,各メトリックの異常への貢献度を計算する. 貢献度の計算には,協力ゲーム理論のShapley Valueに基づくSHAPを採用した. いくつかのSHAPのアルゴリズムの中でも,あらゆる機械学習モデルに適用可能であるKernel SHAPを採用している.

4. 評価

4.1 実験環境

ベンチマークアプリケーション

実験環境として,Google Kubernetes Engine(GKE)上にマイクロサービスのベンチマークアプリケーションである Sock Shop[13]を構築した. Sock Shopを構成する11コンテナからcAdvisorを用いてCPU使用率などのメトリックを5秒おきに収集した.

異常注入

Sock Shopアプリケーションに対して,擬似的な負荷を生成するために,Locust[14]を利用した. Locustによる負荷を生成している間に,システムの異常を模倣するために,特定のコンテナ(user-dbコンテナ)にstress-ngコマンドを用いてCPU負荷を注入した.

ベースライン手法

原因診断のベースライン手法として,Gaussian Based Thresholding(GBT)を用いた. GBTを用いた原因診断は,学習データの平均値とテストデータの平均値の差分が大きい順番に異常への貢献度が高いとする単純な手法である.

4.2 異常の原因診断

Sock Shopを構成する11のコンテナから取得したメトリック数の合計は601であった. そこからTSifterを用いてメトリックのフィルタリングを行った結果,メトリック数が72まで削減された. 今回の実験条件において,異常の原因となるuser-dbコンテナのメトリックは,フィルタリング後に5つのメトリックが残っており,そのグラフを描いたものが下の図である. 各メトリクスは標準化している.

f:id:hirotsuru314:20210525141228p:plain
図2 user-dbコンテナのメトリック

図2から横軸が190あたりからメトリックが大きく変動しているのがわかるが,これは実際にstress-ngコマンドでCPU負荷をかけた時間帯に相当する. メトリックは,20分間に相当する240点を取得しており,前述の通りフィルタリング後の特徴量数は72であるため,分析対象のデータは,72×240の多次元時系列データとなる. このデータを10分間隔で区切って,前半の120点を学習データ,後半の120点を異常を診断するテストデータとして分割する. この場合,異常を検知した時,直近10分間とそれ以前とのメトリックのずれを分析していることになる. 今回は学習データとテストデータの分割を手動で設定しているが,その分割方法をデータから適応的に決定することは今後の課題とする.

下の図は,一つのタイムステップ(図2の横軸が200の点)における異常への貢献度をSHAPにより求め,その結果を可視化したものである.

f:id:hirotsuru314:20210525141352p:plain
図3 SHAPによる解釈(force plot)

横軸は,PCAにより求めた異常スコアである. 学習データの異常スコアのベース値が0.99程度であるのに対し,テストデータでは異常スコアが1.42まで増大している. 図3は,その異常スコアの増大への特徴量の貢献度を可視化したものである. 赤色で示されているものが異常スコアが増大する方向に押し上げている特徴量を示している. 特徴量名は「c-(コンテナ名)_(メトリック名)」となっている. この結果から,図2の横軸が200の時点において,user-dbコンテナのCPUのメトリックが異常に最も寄与していると解釈しており,user-dbコンテナのネットワーク関連のメトリックが次に貢献度が大きくなっている.

次に,下図はテストデータの120タイムステップのSHAPによる解釈の結果をまとめたsummary plotと呼ばれるものであり,異常への貢献度が大きいものから順に10つのメトリックを上から並べたものである.

f:id:hirotsuru314:20210525141634p:plain
図4 SHAPによる解釈(summary plot)

図3と同様に,図4においても,user-dbコンテナのCPUおよびネットワーク関連のメトリックが大きな貢献度を示している. 本実験における異常は,user-dbコンテナにstress-ngコマンドを用いてCPU負荷を注入したものである. そのため,本実験条件において,SHAPによる解釈は,実際の異常の根本原因と一致した結果を与えている.

一方,ベースライン手法であるGBTベースの原因診断では,提案手法と異なる診断結果が得られた. GBTによる解釈,すなわち学習データとテストデータの平均値の差分が大きい順にメトリックを並べた結果は以下の通りである.

c-carts_fs_usage_bytes
c-orders_fs_usage_bytes
c-carts-db_memory_rss
c-carts-db_memory_working_set_bytes
c-orders-db_memory_rss
c-catalogue_fs_usage_bytes
c-user-db_cpu_user_seconds_total
c-user-db_memory_working_set_bytes
c-front-end_cpu_cfs_throttled_seconds_total
c-user-db_network_transmit_packets_total

本実験における根本原因であるuser-dbのCPUのメトリックは異常への貢献度が7番目となっていた. GBTによる解釈では,正常時でも分散が大きいメトリックが,偶発的に学習データとテストデータの平均値の差分が大きくなった場合,それを異常による変動と見分けることができない. 一方,PCAでは正常時の分散の大きさを考慮しているので,本実験では結果的に正しい解釈を与えられているのではないかと考えている.

本実験の異常パターンにおいては,提案手法で採用しているSHAPの方がベースライン手法よりも良い原因診断の結果を与えることがわかった. 一方,本実験は一つの異常パターンのみを対象としているため,提案手法の有用性を示すためには,より広範な異常パターンに対して評価を行う必要があると考える. しかし,単純な平均値の逸脱により異常を解釈するベースライン手法と協力ゲーム理論に基づくSHAPによる解釈が異なる結果を与えうることは大変興味深い結果である. 今後はより広範な障害パターンに対しての評価を進めていく予定である.

4.3 実行時間

異常の原因診断では,その診断結果をシステム管理者が見たのちに異常状態からの復旧のための対応を行うことを想定しているため,診断結果を出すまでの時間は短いことが望まれる. そこで,提案手法において原因診断の結果を出すまでの時間を評価した.

下の表は提案手法の実行時間の結果を実行ステップごとに示したものである. この結果は,図4のsummary plotを計算する際の実行時間である.

実行時間(s)
Step 1:フィルタリング(TSifter) 1.78
Step 2:PCAの学習 0.02
Step 3:SHAPの計算 62
合計 64

提案手法の実行時間は64 sであり,SHAPの計算が支配的であることがわかった. SHAPの計算にはSHAPの提案者らが開発しているPython製のライブラリ[15]を用いた. また,summary plotでは120のタイムステップにおけるSHAPの計算が実行されるため,今回は8コアサーバで並列処理している.

提案手法が実用的なシステムに有効かどうかを議論するためには,システムの規模(分析対象のメトリック数)を変化させた際の実行時間について評価する必要がある. 提案手法の実行時間は,対象となるメトリック数の増大とともに長くなるため,本実験条件より大規模なシステムへ対応するためには,SHAPの計算の高速化が課題となる. SHAPの計算時間の改善策として,並列数を大きくすることや,shapライブラリ[15]の計算時間を改善することなどを考えている.

5. 今後の展望

まず,提案手法の有用性を示すためには広範な異常パターンに対して評価を行う必要があると考える. そのため,異常の注入やメトリックの取得などを自動化した実験環境を構築し,広範な異常パターンに対して原因診断の精度を定量的に評価する予定である.

次に,対象とするシステムが大規模化した際に,提案手法が実用に耐えうるかを検証する. そのために,メトリック数が増大した場合の原因診断の計算時間の評価を行う.

最後に,今回の採用したPCAは,線形かつ時系列の情報を考慮していない単純なモデルであるため,今後,
非線形のモデルや時系列に対応したモデルを採用し,その有効性を検証する予定である.

スライド

speakerdeck.com

参考文献

[1] A Deep Neural Network for Unsupervised Anomaly Detection and Diagnosis in Multivariate Time Series Data

[2] Detecting and Localizing Anomalies in Container Clusters Using Markov Models

[3] Microscope: Pinpoint Performance Issues with Causal Graphs in Micro-service Environments

[4] AutoMAP: Diagnose Your Microservice-based Web Applications Automatically

[5] A Unified Approach to Interpreting Model Predictions

[6] Peeking Inside the Black-Box: A Survey on Explainable Artificial Intelligence (XAI)

[7] "Why Should I Trust You?": Explaining the Predictions of Any Classifier

[8] Shapley Values of Reconstruction Errors of PCA for Explaining Anomaly Detection

[9] Explaining Anomalies Detected by Autoencoders Using SHAP

[10] On Anomaly Interpretation via Shapley Values

[11] TSifter:マイクロサービスにおける性能異常の迅速な診断に向いた時系列データの次元削減手法

[12] Principal component analysis

[13] Sock Shop A Microservices Demo Application

[14] LOCUST: An open source load testing tool

[15] slundberg/shap

企業研究者の立場からKaggleに取り組む意義を考えた

先日,Kaggleで初めてコンペに挑戦し,その振り返りをブログに書きました. 現在,企業で研究者として働いている私は,Kaggleのコンペに取り組むことは非常に学びが多く,自身の研究活動にも良い貢献をするだろうと確信しました. 私自身Kaggleに取り組むまでKaggleと研究に繋がりを見いだせておらず,実際にコンペに取り組むことでその繋がりが見えてきました. まだ1コンペしか参加経験がないビギナーではありますが,私が考えている研究者としてKaggleに取り組む意義について現時点の考えをまとめたいと思います.

自分の立場について

Kaggleへの取り組みが研究に良い貢献をするかどうかは,研究の分野や内容に依存すると思うので,私の立場をはっきりさせておきます. 私は現在,インターネットインフラサービスを提供するさくらインターネット株式会社の組織内研究所であるさくらインターネット研究所で研究員として働いています. 私は,自分たちのドメインにおける課題の中で人間の手や簡単なルールベースのアルゴリズムでは解決が難しい課題に対して,機械学習を適用して解決する研究に取り組んでいます. したがって,私は機械学習の理論の研究をしているわけではなく,機械学習を応用して研究を行っている立場です.

機械学習を応用する研究

機械学習を課題解決の手段として用いる研究は,多くの分野で増えていると思います. 私が参加した2020年11月開催のIBIS2020の発表の中でも,農業・土木工学・機械設計など様々な分野に機械学習が応用されていました. 私の分野においても例外でなく,2020年に開催されたクラウドコンピューティング技術に関する国際会議であるIEEE CLOUD 2020において,機械学習をはじめとした数理的アプローチで課題解決を行っている研究が多数発表されていました.

我々の研究所ではこういった流れに取り残されず研究を行っていくためにも,研究の問題設定を行うドメイン知識や経験が豊富な研究員と,その課題を機械学習統計学を駆使して解決する研究員で役割分担をしてチームとして研究を行う取り組みを積極的に行っています.以下はその役割分担について,同僚の松本さん発表資料から抜粋したものです.

f:id:hirotsuru314:20210114111528p:plain

Kaggleを知って危機感を覚えた

前述のような自分たちのドメインにおける課題に対して機械学習を適用する研究において,機械学習を用いた問題解決力が研究に直結することは明らかです. 一方,私がKaggleの存在を知ってから,公開されている過去コンペの解法についての発表資料をいくつか読んでみると,その中の多く解法が私が知らない,もしくは名前だけ知っているけど理論や使い方もわからない手法を用いていました.(例を挙げると,LightGBM・BERT・Metric Learningなどなど) 私はこのことに非常に危機感を覚えました.というのも,もし研究における問題設定がうまくできていたとしても,問題解決力が不足していて,実用上有効な精度が出せないと研究が進まないといったことが起きると思います. 例えば,我々の分野においてWebシステムの異常検知という問題設定をした場合,その検知精度が60%であったとしたら,実用性があまりにも低いので,どれだけ素晴らしい問題設定や提案をしたとしても研究として成立しないと思います. 世界中のデータサイエンティストが精度を競い合うKaggleにおいて頻出のモデルやアルゴリズムをキャッチアップできていない今の自分の力ではそういった問題が頻繁に起きるのではないかと思いました.

Kaggleを通じて問題解決力を高める

そこで,Kaggleを通じて機械学習を用いた問題解決力を高めることは,自分の研究活動において有効だと考えました. Kaggleで得られる知識や経験はいろいろありますが,私が特に研究に直結すると思ったのは以下の3点です.

1.最新かつ実用的な手法のキャッチアップ

Kaggleでは枯れた技術のみが使われているわけではなく,最新の研究や論文の成果がどんどん実践で使われています.逆に言うと,state-of-the-artを追わないと勝てないのかもしれません.このため,コンペに参戦することや参戦後に上位解法を学ぶことで,数年以内に提案されたような新しい手法をキャッチアップできると思いました. さらに,新しいだけでなく,Kaggleのコンペで実データに適用されて,精度が検証されていることが非常に有用だと思いました. 手法の提案論文では,(査読の段階で取り除かれることも多いと思いますが)提案に有利な評価であったり,実用的でない擬似データを用いた評価であることが多いです.それがKaggleで実戦投入されて,論文とは異なる実データで精度や実用性が検証されていることには非常に価値があると思いました.

2.実装力

Kaggleのコンペの中でも「Code Competition」と呼ばれるコンペでは,与えられたNotebook環境で成果物をSubmitする必要があるため,メモリの制限やコードの実行時間に制約が課されます.これにより,省メモリかつ高速な実行コードを書く力要求されます. また,前述のような最新の手法の中にはライブラリが公開されていないものもあるため,時には論文を読みながら自分で実装する力も必要になると思います.

3.データの捉え方

Kaggleでは様々な種類のデータ・問題設定でコンペが開催されるため,継続的に参加することで幅広いデータの捉え方を学べると思いました. これは自分たちの手持ちのデータだけを分析していてはなかなか得られない経験なので,非常に有用だと思います.

また,これらの習得できる力を踏まえた問題解決力が,コンペの順位,メダル,称号などを通して客観的な指標として評価できることも魅力の一つだと思います.Kaggleで学べることはまだあると思いますが,上記の3つだけでも,どれも自分の研究にも必要な力であるので,私は研究に必要な機械学習を用いた問題解決力の習得・維持のためにKaggleに継続的に取り組んでいきたいと考えました.
ここまでが企業研究者の立場である私がKaggleに取り組んでいこうと思う主な理由なのですが,Kaggleでコンペに取り組むことが直接研究の業績に繋がるような事例もあるようなのでまとめました.

Kaggleの成果が直接的に研究業績に繋がる事例

ここからは,Kaggleのコンペで上位入賞するような場合に限った話ですのでハードルは高いです.事例を以下の3点にまとめました.

1.Kaggleの解法を論文としてまとめる

  • 「PLAsTiCC Astronomical Classification」というコンペの優勝者が書いた論文がジャーナルに採択

    iopscience.iop.org

  • Kaggleの過去コンペの結果をまとめたようなジャーナル論文もある

    www.sciencedirect.com

2.コンペで入賞して学会で発表する

  • 私が前回参加した「Riiid! Answer Correctness Prediction」では上位入賞者は「AAAI-2021 Workshop on AI Education」で発表できる.

  • 国際学会が主催しているコンペもある.例えば,IEEEが主催した「IEEE-CIS Fraud Detection」で上位入賞者はIEEE CIS Conferenceで発表できる.

3.コンペで入賞してプレスリリースを出す

多くの企業がKaggleのコンペで上位入賞した場合,以下のようにプレスリリースを出しているようです.

www.nssol.nipponsteel.com

www.jri.co.jp

さいごに

私は自分の研究を進めていく上で必要な機械学習を用いた問題解決力の習得・維持のためにKaggleに継続的に取り組んでいきたいと考えています.いずれはコンペで上位入賞できるような力をつけ,Kaggleの成果自体が研究業績になるようにしていきたいと思います.このブログが私と同じような立場や危機感を持っている方々の参考になれば嬉しいです.Kaggleやっていきましょう!

Kaggle初コンペの振り返り〜Riiidコンペで銀メダル獲得〜

Kaggleで開催されていたRiiid! Answer Correctness Predictionに参加しました.結果を簡単にまとめると以下の通りです.

  • 順位:139位(3406チーム中)
  • メダル:銀メダル(上位5%以内)
  • 解法:LightGBMとSAKTのアンサンブル
  • チーム:1人で参加

自身初となるデータ解析コンペでしたが,なんとか銀メダルを獲得することができました.今回はその振り返りを時系列で書いていきたいと思います. 先に感想を述べると,Kaggleは本や論文ではなかなか得られない学びに溢れている上に,世界中の人々と順位を競い合うゲームのような感じでとにかく楽しいので最高すぎました. この記事がこれからKaggleを始める方々の参考になれば嬉しいです.

目次

Kaggleとの出会い

Kaggleの存在を知ったのは4年くらい前ですが,Kaggleを始めるきっかけとなったのは2020年11月に開催された第23回情報論的学習理論ワークショップ (IBIS2020)という学会でした.IBIS2020のチュートリアルセッションにて「Kaggle・実践データ解析入門」という発表がありました.この発表では,Kaggleの称号の最高位であるGrandmasterを持つ登壇者の方が,Kaggleの概要や,「こういうデータのコンペでは初手としてこういう手法が有効である」といった実践的なテクニック,研究や論文の成果をどのようにKaggleに使っているか,など非常に興味深い内容が盛り沢山でした.動画や資料は公開されていませんが,Kaggleの概要の部分に関しては,発表中で紹介されていた以下のスライドがまとまっています.

結局、Kagglerは何を必死にやっているのか? / What is Kaggle? - Speaker Deck

この発表を聞いて,単純に「Kaggle面白そう!」と思いました.さらに,Kaggleでは最新の研究や論文の成果がどんどん実戦で使われていることを知り,Kaggleと研究に繋がりが見えたことが始める大きなきっかけとなりました.実際に登壇者の方は常に最新の論文から新しく提案された機械学習モデルをチェックしていて,精度が良さそうなものがあれば,論文を読みながら実装してKaggleで実戦投入するといったことをやっているようです.現在私は,さくらインターネット研究所で,機械学習をソルバーとして用いて自分たちのドメインクラウドコンピューティングなど)における課題を解決する研究に取り組んでいるため,機械学習を用いた問題解決能力は研究に直結します.そのため,Kaggleを通して最新のモデル,かつ提案論文内での評価とは異なるデータで実用性が検証されているモデルを把握しておくこと,また実践を通じて自身の問題解決能力を高めることは,研究活動に大きく貢献するだろうと確信しました.そんなこんなで興味を持ってすぐにアカウントを作ってKaggleを始めました.

コンペ参戦準備

Kaggleを始めるといっても,コンペへの参加の仕方さえ分からなかったので,まずはKaggleの全体像を掴むところから始めました.具体的には,以下の本を一通り手を動かしながら読みました.

実践Data Scienceシリーズ PythonではじめるKaggleスタートブック (KS情報科学専門書)

この本は,Kaggleをゼロから始める人にとっては非常に良い本で,アカウント作成からコンペへの参加,予測結果のSubmitまでの流れを一通り学ぶことができるため,初めてのコンペ参加へのハードルをかなり下げてくれます.

コンペに登録

Kaggleスタートブックを読み終えた時点ですぐにコンペに参加登録をしようと思い,良さそうなコンペを探しました. Kaggleスタートブックの中で初めてのコンペを選ぶときのコツとして,「すでに開催から一定時間が経過していてコンペ終了まで一ヶ月程度のコンペがおすすめ」と書いてありました.その意図としては,開催から一定時間が経過していると,Discussionや公開Notebookなどの充実した情報が集まっているので,それを参考に自分のモデルを作ることができるためです.これを参考にRiiid! Answer Correctness Predictionに登録しました.コンペの概要は以下の通りです.

  • 開催期間:2020/10/06〜2021/01/08
  • データ:TOEICの学習アプリで各ユーザが過去に取り組んだ問題とその問題に正解したかどうかのデータ
  • 問題:将来あるユーザはある問題に正解できるかを予測する二値分類問題

このような人間の知識の習得をモデリングする問題はKnowledge Trackingというそうで,教育分野では重要な問題なようです.
私がコンペに登録したのが2020/11/25で締切まで約1ヶ月半程度全力で取り組みました.ちなみに,取り組んだ時間としては,毎日寝る前2~3時間程度で,大晦日は格闘技を見るためにやりませんでしたが,それ以外は休まず毎日やりました.

初Submit

まず簡単な特徴量とモデルでSubmitまでやってみようと思ったのですが,このSubmitまでが鬼門でした. 本コンペは,KaggleのNotebook環境を用いてSubmitする必要がある「Code Competition」と呼ばれるコンペでした.このコンペの特性上,以下の2点がSubmit時の大きな制約となりました.

  1. Notebook環境のメモリ上限は16GB

  2. Notebook環境上での予測時間は9時間以内(テストデータは250万行)

この2点の制約により,機械学習やデータサイエンスの力だけでなくエンジニアリング力も試されました. まず,メモリ上限についてですが,学習データのCSVファイルを普段通りpandasでread_csvしたところ,読み込みに数十分かかる上にメモリが上限の16GBを超えてしまいエラーになりました.というのも学習データのCSVは,行数が1億を超えていました.さっそくデータの読み込みさえできずに焦りましたが,公開Notebookに「Tutorial on reading large datasets」という素晴らしいNotebookがあり,データ型をきちんと指定することや,読み書きが早いファイル形式を利用することなど知見が溜められていました.

次に,予測時間の制限についてですが,私は最初,学習データで作った独自の特徴量をPandasのDataFrameで保持し,テストデータに対してPandasのmerge関数でLeft joinすることで特徴量を追加していましたが,どうやらこのDataFrameのjoinの処理が遅いらしく,推論時間を9時間以内おさめることができませんでした.これもまた素晴らしいNotebookに助けられて「LGBM with Loop Feature Engineering」に特徴量を辞書型で保持し,テストデータのDataFrameの行に対してfor文を回し,特徴量を追加していくといった手法でした.これを真似ると推論時間は大幅に改善されました.このように悩んだことがあっても多くのKagglerが有益な情報を公開してくれていて学びながら取り組んで行けることもKaggleの魅力だと実感しました.(補足:Pandasのleft joinを357倍高速化したというNotebookも公開されていてこちらも参考になりそうですがコンペ中は追いきれませんでした)

これらの問題を解決して,やっとSubmitまで行き着きました.モデルには,LightGBMを選択しました.LightGBMは,ここ数年のテーブルデータのコンペの初手として多くのKagglerが用いているようです.なぜ初手にLightGBMを用いるかについては,以下の記事にわかりやすくまとまっていて非常に参考になります.

upura.hatenablog.com

初Submit後にLeaderboard(コンペの順位表)に自分の名前が載っているのを見てとにかくモチベーションが上がりました.

特徴量エンジニアリング

コンペ期間において一番時間を割いたのは特徴量エンジニアリングです.まずは簡単に思いつくようなユーザの正解率や,問題の正解率など,直感的に有用そうな特徴量を追加していきました.あとから知ったのですが,これはTarget encodingと呼ばれる特徴量生成のテクニックでした.Target encodingとは,説明変数に含まれるカテゴリ変数と目的変数を元にして特徴量を作り出すことです.
ここでTarget encodingにより明らかに有効そうな特徴量を追加して手元のスコアが上がっているのに,Submitしたスコアが全然上がらないという問題に悩まされました.これはTarget encodingにおける定番的なミスで「リーク」という問題を起こしていました.例えば,ユーザの正解率という特徴量を作るときに,学習データ全てを使って算出したユーザの正解率を特徴量として使うと,ある時点では将来の情報が漏れていることとなってしまい,精度が過剰に高くなってしまいます.そこで,ある時点のユーザの正解率はそれより過去のユーザの問題の回答から正解率を算出するように変更したところ,手元とSubmit時のスコアに相関が取れるようになりました.私はこれまで数年ほど機械学習の勉強をしてきたのにいざ実戦でやってみると知らないことが多いなと辛い気持ちになりましたが,これが成長というやつです.
自分がコンペに取り組む上で重要な知識が欠落していること認識したので,知識を補うために以下の本を読みました.

Kaggleで勝つデータ分析の技術

この本は,Kaggleのコンペに取り組む上で必要な情報が網羅的に書かれていてめちゃくちゃ参考になりました.前述のTarget encodingのリークの話もしっかり書いてありました.この本を参考に思いつく特徴量をどんどん追加していきました.最終的な特徴量の数は30~40程度でした.

モデルの構築

モデルは,最初のSubmitの時に採用したLightGBMを主に使っていました.特徴量エンジニアリングの成果が出て,早い段階でLightGBMのシングルモデルで銅メダル圏内に届きました.しかし,特徴量エンジニアリングが行き詰まってきてスコアの伸びが小さくなってきたので,モデルのパラメータチューニングに着手しました.パラメータチューニングにはOptunaを用いました.

github.com

Optunaには,LightGBMに特化したLightGBM Tunerという拡張機能があり,Pyhtonコードのimport文を1行変更するだけで簡単にチューニングができて便利でした.しかし,このチューニングはスコアにはあまり貢献しませんでした.
この点については,コンペのDiscussionの「LGBM Parameter Tuning」で議論が展開されており,パラメータチューニングは大きく精度に貢献しないことなどが書かれていました.このようにコンペ開催中にも惜しみなく知識や経験を共有する文化は素晴らしいと思いました.
締切まで残り3週間程度の時期にLightGBMのスコアが停滞してきたので,この時点で他の学習モデルの検討を始めました.

アンサンブル学習

LightGBM以外の学習モデルには,TransformerをベースにしたSelf Attentive Knowledge Tracing(SAKT)を選択しました.その理由は,公開されているNotebookの中で,LightGBMなどの勾配ブースティング系とSAKTのアンサンブル学習で大きくスコアを伸ばしているものがちらほら見つかったためです.SAKTは以下の提案論文と公開Notebookを読みながら実装しました.

arxiv.org

SAKTの実装は完了したものの,LightGBMのスコアには到達しませんでした.しかし,アンサンブル学習を試してみると一気にスコアが伸び,そのときの順位で100位くらいで銀メダル圏内に到達しました.アンサンブルの方法は単純な加重平均で,この時点でのスコア的にはAUCでLightGBMが0.778,SAKTが0.772で5:5のアンサンブルをすると一気に0.786まで上がりました.
このようにモデルのアンサンブル学習を行うと,単体のモデルの最高スコアを大きく上回ることが起きうるというのは「Kaggle Ensembling Guide」を読んで知っていたのですが,どうも直感的には低いスコアの方に引きづられそうな気がしていました.ただ,結果を見ると,これがアンサンブル学習の力か…となりました.また,ここでスコアを大きく向上させられた大きな要因として,アンサンブル学習を行ったのが勾配ブースティング系のモデルと,ニューラルネットワーク系のモデルの異なる種類のモデルをアンサンブルさせたことだと思われます.これもアンサンブル学習にはモデルの多様性が重要であるとよく語られることではあるのですが,本などで学んだこともこうやって実際に体験すると,より知識として定着するので実践から学ぶのは良いなと改めて思いました. あとはコンペ終了まで,ひたすら特徴量エンジニアリングを頑張ってLightGBMとSAKTのアンサンブルでスコアを出すというのを繰り返していました.

コンペ終了:Leaderboardの公開

2021/01/08の午前9時にコンペが終了し,最終的なLeaderboardが公開されました, 結果は3406チーム中139位でなんとかソロで銀メダルを獲得できました.

f:id:hirotsuru314:20210108230540p:plain

コンペ全体を通しての感想と今後

初コンペを終えた感想としては,Kaggleはとにかく学びが多くて,とにかく楽しいので,これからもKaggleを続けていこうと思いました.当初は初コンペでメダルが取れたら嬉しいなーくらいの気持ちでしたが,いざ銀メダルを獲得すると,嬉しさよりも,上位陣との歴然の差を目の当たりにした悔しさの方が勝っていました. 銀メダルというと一見あとちょっと頑張れば金メダルだと思えるかもしれませんが,そこには大きな壁があり,コンペ期間中一度も金メダルにかすりそうな瞬間はありませんでした.さらにコンペ終了後に公開された上位の解法にざっと目を通しましたが,どんな頭しているんだと思うくらいすごくてさらに悔しい気持ちがこみ上げてきました.ただ,こういった悔しい経験がどんどん自分を強くしていくので勝てなくても参加するのが大事だと思いました.これからコンペの上位解法をしっかり読んで勉強します.
今後は,自分もコンペで優勝争いをしたいし,いずれは最高位であるGrandmasterにもなりたいです.そのためにもまずは次のコンペでもメダルをとってKaggle Expert(銅メダル以上2個)を目指します.
長くなりましたが最後まで読んでいただきありがとうございました.Kaggleやっていきましょう!

グラフィカルモデルに基づく因果探索手法の調査

最近,因果推論や因果探索に興味を持ち,勉強している.というのも最近,ゆううきさん と一緒に分散システムの異常の原因を即時に診断するための研究を進めている.原因を診断するためのアプローチとして,サーバやコンテナ等から取得できる様々なメトリック(CPU使用率やメモリ使用率など)を(グラフ理論における)ノードとして,因果グラフを構築することを考えている.メトリック同士の単なる「相関」ではなく,結果と原因の関係である「因果」を捉えようとするアプローチである.例えば,システムの障害が発生した場合,相関だけでは,AとBが関連がありそうというところまでしか言えないが,因果を特定できると理想的には,Aの原因はBであるといった議論ができるため,有用だと考えている.
実際に,前述のような因果グラフを構築して障害の原因を特定しようというアプローチは,以下の例に挙げるようにここ数年で増えている印象がある.

今回の記事は,因果探索手法の中でもグラフィカルモデルに基づくアプローチに絞って,まとめていく.以下のサーベイ論文がよくまとまっていたので全面的に参考にした.

Review of Causal Discovery Methods Based on Graphical Models(2019)

グラフィカルモデル

グラフィカルモデルとは,確率モデルをグラフを用いて記述したものである.ここでいうグラフとは, x軸, y軸があるグラフではなく,グラフ理論のことであり,ノード(頂点)とエッジ(枝)の集合で構成されるものである.グラフィカルモデルを用いることで,どの確率変数同士が直接的に関係しているかを視覚的に表現できる.
上記でエッジと書いた部分は,単なる線と矢印の二つの表現があり,それぞれを用いたグラフを無向グラフおよび有向グラフと呼ぶ.変数間の因果を表現する際に利用されるのは主に後者の有向グラフである.有向グラフの代表的な例としてベイジアンネットワークが挙げられる.

因果探索

因果探索とは,データから因果関係(因果グラフ)を推定する手法である.前述のグラフィカルモデルで考えると,考慮すべき確率変数(ノード)のデータは手元にある状態で,因果の構造,すなわちどのようにエッジが繋がっているかを推定するということになる.多くの分野において,データの根底にある因果関係を見つけ出し,それを利用することは非常に重要であるため,因果探索は,ゲノミクス,臨床医学,疫学など幅広い分野に応用されている. 次から,因果探索の具体的なアプローチについてまとめていく.

制約ベースの因果探索

制約ベースのアプローチでは,まず,観測変数にどのような条件付き独立性が成り立つかを推測する.次に推測された条件付き独立性を制約として,それを満たすような因果グラフを探索する.ここでは,代表的な制約ベースのアルゴリズムとして,PCアルゴリズム[1]を説明する.

PCアルゴリズム

PCアルゴリズムの処理の流れを以下の図を用いて説明する.前提として,PCアルゴリズムでは,未観測共通原因(潜在的な交絡因子)がないことを仮定している.

f:id:hirotsuru314:20201007095925p:plain
PCアルゴリズムによる因果グラフの推定

  1. 完全無向グラフを構築する.
  2. 条件付き独立性が成立するエッジをグラフから除外する.
    ここで,条件付き独立性を検索する際は,条件として与える変数の数を0から1ずつ増やしながら,条件付き独立性が成立しうる全ノードの組合せを対象に行われる.
  3. V-structureに基づきエッジの方向を決定する.
    V-structureでは,三つの変数A,B,Cが「A-B-C」でA-C間にエッジがなく,A,Cを条件付き独立とする条件ノード群にBが含まれていないならば,「A→B←C」の関係が成立するとして,無向グラフ中の一部の因果の方向を推定する.
  4. オリエンテーションルールに基づきエッジの方向を決定する.
    オリエンテーションルールは有向非巡回グラフ(DAG)の定義から導かれるグラフの構造に関するルールである.具体的なルールについては「Estimating High-Dimensional Directed Acyclic Graphs with the PC-Algorithm」にまとめてある.

PCアルゴリズムは,一般にノード数が増加すると上記2の条件付き独立性を検索する際のノードの組合せの数が増大し,計算時間がかかるという問題がある.これに対して,条件付き独立性の検定をGPUで並列化して高速に処理するといった提案[2]もあるようだ.また,PCアルゴリズムは未観測共通原因がある場合のFast causal inferece(FCI)アルゴリズム[3]などへ拡張されている.

スコアベースの因果探索

スコアベースのアプローチでは,同じ条件付き独立性を与える因果グラフの集合であるマルコフ同値類ごとに,モデルのよさを測るスコアをつける.ここでは,代表的なスコアベースのアルゴリズムとして,Greedy Equivalence Search(GES)アルゴリズム[4]を説明する.

GESアルゴリズム

GESアルゴリズムでは,PCアルゴリズムのように完全無向グラフから始めるのではなく,空のグラフから始まり,必要なエッジを追加したり,不要なエッジを削除していきながら因果グラフを構築する.アルゴリズムの各ステップでは,グラフに有向エッジを追加することで,ベイズ情報量規準(BIC)のようなスコアや統計的仮説検定のZスコアで計測される適合性が向上するかどうかが決定され,適合性を最も向上させるエッジが追加される.スコアをそれ以上改善できなくなった場合,GESアルゴリズムは,どのエッジを削除するとスコアが最も改善されるかをエッジごとに検証していく.このような手順でスコア(例えばBIC)を用いて因果グラフを推定していく.

構造方程式モデルに基づく因果探索(非ガウス性・非線形性の活用)

先に挙げた制約ベース・スコアベースの因果探索は,データの構造に対してできるだけ仮定をおかずに,データのみから因果グラフを推測するアプローチである.これは,後述する構造方程式モデルにおける関数形と誤差変数の分布に仮定をおかないノンパラメトリックアプローチに分類される.仮定をおかないため様々なデータに適用できるという利点があるが,因果の構造を正確に推定できないケースも多い.そこで,分析者の事前知識を仮定としてモデルに取り入れるアプローチがある.ここでは,仮定を取り入れる上で重要な概念である構造方程式モデルについて述べ,非ガウス性の仮定を利用したLiNGAMという手法について述べる.

構造方程式モデル

構造方程式モデルは,データの生成過程を記述するツールである.ここで, x yという二つの変数について,以下のような単純な因果グラフを考える.

f:id:hirotsuru314:20201007131606p:plain
因果グラフのサンプル

このような因果グラフにおける構造方程式を以下のように表す.

 \displaystyle
x=e_x \\
y=f_y(x, e_y) \\

ここで, x yは観測変数で, e_x e_yは未観測の変数であり,誤差変数という. e_y yの値を決定するために寄与しうる x以外の全ての変数をまとめて表した変数という意味をもつ.因果グラフをこのような方程式で表現した場合,分析者は以下の2点に仮定を導入することができる.

  1. 関数 f_yの仮定(例えば線形か非線形か)
  2. 誤差変数 e_x e_yの分布の仮定(例えば,ガウス分布か非ガウス分布か)

上記の二つに仮定を導入するものをパラメトリックアプローチ,関数には仮定をおく一方,誤差変数の分布には仮定をおかないものをセミパラメトリックアプローチと呼ぶ.ここでは,代表的なセミパラメトリックアプローチであるLiNGAMという手法を紹介する.

LiNGAM

LiNGAMはLinear Non-Gaussian Acyclic Modelの略であり,関数形に線形性を,誤差変数の分布に非ガウス分布を仮定するアプローチである.一見,二つの仮定をおいているようにも見えるが,非ガウス分布という仮定はガウス分布でさえなければどんな(連続)分布でも良いということからほぼ仮定をおいていないという解釈ができ,セミパラメトリックアプローチに分類されるようだ.LiNGAMについては,機械学習プロフェッショナルシリーズの「統計的因果探索」に詳しく解説されている.ここではさらっとアルゴリズムの概要をまとめる.
まずLiNGAMにおける仮定を整理する.

  • 関数形が線形
  • 誤差変数が非ガウス連続分布
  • 有向非巡回グラフ(DAG)
  • 誤差変数が独立性(未観測共通原因がないことを意味する)

p個の観測変数 x_1,x_2,...x_pに対するLiNGAMモデルは構造方程式モデルを用いて以下のように表される.

 \displaystyle
x_i=\sum_{i \neq j}^p b_{ij}x_i+e_i

これを行列・ベクトルで表現すると以下のようになる.

 \displaystyle
\mathbf{x}=\mathbf{B}\mathbf{x}+\mathbf{e}

ここで, p次元ベクトル \mathbf{x}は観測変数ベクトル, p× p行列 \mathbf{B}は係数行列, p次元ベクトル \mathbf{e}は誤差変数ベクトルである.そして, \mathbf{e}の成分 e_iは独立で,それぞれ非ガウス連続分布に従う. 係数行列 \mathbf{B}は,観測変数の分布 p(x)に基づいて一意に推定可能であり, \mathbf{B}の成分のゼロ・非ゼロパターンから因果グラフを描くことができる. 係数行列の推定には,独立成分分析(ICA)によるアプローチと回帰分析と独立性評価によるアプローチがある[5]. どちらのアプローチもcdt15/lingamで実装されている.OSSのコードも含めてLiNGAMはもう少し詳細に追いかけたい.

さいごに

今回は,サーベイ論文を参考にグラフィカルモデルを用いた因果探索の手法を整理した.ここでいくつかアルゴリズムを紹介したが,最近機械学習技術を活用した多くの手法が提案されている.例えば,深層強化学習[6],敵対的生成ネットワーク(GAN)[7]を活用したものなどがあり,まだまだ奥が深い.また,重要なイシューとして,確率変数が時系列データの場合にどのように時系列の情報を消失せずにモデリングするかというのも大変興味深い.今後も因果探索について着実に勉強していって,ブログにまとめていきたい.

参考

[1] An Algorithm for Fast Recovery of Sparse Causal Graphs
[2] cuPC: CUDA-Based Parallel PC Algorithm for Causal Structure Learning on GPU
[3] Causal inference in the presence of latent variables and selection bias
[4] Optimal structure identification with greedy search
[5] 統計的因果探索
[6] Causal Discovery with Reinforcement Learning
[7] SAM: Structural Agnostic Model, Causal Discovery and Penalized Adversarial Learning

経路積分量子モンテカルロ法を用いた量子アニーリングのシミュレータ

0または1の状態しかとれないビットを用いた古典的な計算機において量子力学の原理に基づいた量子アニーリングをシミュレートできるのだろうか.
量子アニーリングの理論について学んだ(前回の記事)のちに,このような疑問が生まれた.調べてみると,どうやら手元のPCでシミュレータを作るには経路積分量子モンテカルロ法という手法がよく用いられるようだった. そのため今回は,経路積分量子モンテカルロ法について学び,Pythonでの既存実装に触れてみることで,量子アニーリングについての理解をより深めていくことを目的とする.

目次

経路積分量子モンテカルロ法

今回,古典的な計算機を用いて量子アニーリングをシミュレートするために,経路積分量子モンテカルロ法という手法を用いる. 量子系をシミュレートするために,まず量子ゆらぎを含む横磁場イジングモデルの分配関数 Zを,対応する古典イジングモデルの分配関数で表す. そのとき,これから示す計算を追っていくと,横磁場イジングモデルの分配関数が1次元高い古典の分配関数に対応することがわかる.この結果を使って量子系の平衡状態における性質を古典系のマルコフ連鎖モンテカルロ法で調べる,というのが量子モンテカルロの基本的な考え方である.

以下,理論的な背景を少し詳細に見ていく.
スピン数が Nの場合,横磁場イジングモデルハミルトニアンは以下のように表される.

 \displaystyle
\begin{eqnarray}
\hat{H} & = & \hat{H_0}+\hat{H_q} \\
\hat{H_0} & = & -\sum_{i,j=1}^NJ_{ij}\hat{\sigma}_i^z\hat{\sigma}_j^z-\sum_{i=1}^Nh_i\hat{\sigma}_i^z \\
\hat{H_q} & = & -\Gamma\sum_{i=1}^N\hat{\sigma}_i^x \\
\end{eqnarray}


ここで, \hat{\sigma_i^z}, \hat{\sigma_i^x}はパウリ行列の z,x成分であり,以下の行列で表される.

 \displaystyle
\hat{\sigma_i^z} = 
\begin{pmatrix}
1 & 0 \\
0 & -1 \\
\end{pmatrix}
,
\hat{\sigma_i^x} = 
\begin{pmatrix}
0 & 1 \\
1 & 0 \\
\end{pmatrix}


このハミルトニアンを用いて分配関数は以下のように表せる.

 \displaystyle
\begin{eqnarray}
Z & = & Tr[e^{-\beta\hat{H}}] \\
& = & Tr[e^{-\beta(\hat{H_0}+\hat{H_q})}] \\
& = & \sum_{\boldsymbol{\sigma_1}}\langle\boldsymbol{\sigma_1}\mid{e^{-\beta(\hat{H_0}+\hat{H_q})}}\mid\boldsymbol{\sigma_1}\rangle \\
\end{eqnarray}


ここで, \beta は温度の逆数(逆温度)である.
Suzuki・Trotter分解により,

 \displaystyle
e^{-\beta(\hat{H_0}+\hat{H_q})} = \lim_{M\rightarrow\infty}\biggl(e^{-\frac{\beta}{M}\hat{H_0}}e^{-\frac{\beta}{M}\hat{H_q}}\biggl)^M


と表せることを用いると,

 \displaystyle
Z = \lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma_1}}\langle\boldsymbol{\sigma_1}\mid\biggl(e^{-\frac{\beta}{M}\hat{H_0}}e^{-\frac{\beta}{M}\hat{H_q}}\biggl)^M\mid\boldsymbol{\sigma_1}\rangle


ここでM個の積の間に以下の恒等演算子を挟み込むと,

 \displaystyle
\boldsymbol{1}  =  \sum_{\boldsymbol{\sigma_k}}\mid\sigma_k\rangle\langle\sigma_k\mid (k = 1, 2, ..., M)


分配関数は以下のように表される.

 \displaystyle
\begin{eqnarray}
Z & = & \lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma_k},\boldsymbol{\sigma_k'}}\langle\boldsymbol{\sigma_1}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_1'}\rangle\langle\boldsymbol{\sigma_1'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_2}\rangle \\
& \times & \langle\boldsymbol{\sigma_2}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_2'}\rangle\langle\boldsymbol{\sigma_2'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_3}\rangle \\
& \times & ・・・ \\
& \times & \langle\boldsymbol{\sigma_M}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_M'}\rangle\langle\boldsymbol{\sigma_M'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_1}\rangle
\end{eqnarray}


ここで, \mid\boldsymbol{\sigma_k}\rangle N 個のスピン系の直積空間を表す.

 \displaystyle
\mid\boldsymbol{\sigma_k}\rangle=\mid\boldsymbol{\sigma^z_{1,k}}\rangle\otimes\mid\boldsymbol{\sigma^z_{2,k}}\rangle\otimes・・・\otimes\mid\boldsymbol{\sigma^z_{N,k}}\rangle


分配関数の式から, k = 1, 2, ..., M は元の問題の持つ変数が定義された空間とは別の独立な変数のインデックスであり,トロッタ方向または虚時間方向と呼ばれる.
下図はSuzuki・Trotter分解を行なった際にトロッタ方向に次元が追加されている様子を示した概略図である.

f:id:hirotsuru314:20191210223523p:plain:w500
Suzuki-Trotter分解の概略図

これまでの結果から,分配関数を求めるには,以下の2つの値を求めればよいことがわかる.

 \displaystyle
\begin{eqnarray}
\langle\boldsymbol{\sigma_k}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_k'}\rangle \tag{1} \\
\langle\boldsymbol{\sigma_k'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_{k+1}}\rangle \tag{2}
\end{eqnarray}


 (1) \hat{H_0}が対角行列であることから,

 \displaystyle
\begin{eqnarray}
\hat{\sigma}_i^z\mid\sigma_k\rangle & = & \hat{\sigma}_{i, k}\mid\sigma_k\rangle \\
\langle\boldsymbol{\sigma_k}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_k'}\rangle & = & exp\biggl[\frac{\beta}{M}\biggl(\sum_{i,j=1}^NJ_{ij}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}+\sum_{i=1}^Nh_i\hat{\sigma}_{i, k}\biggl)\biggl]\prod_{i=1}^N\delta(\hat{\sigma}_{i, k}, \hat{\sigma}_{i, k}')
\end{eqnarray}


となる.次に (2) は以下の式が成立する.

 \displaystyle
\begin{eqnarray}
\langle\boldsymbol{\sigma_k'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_{k+1}}\rangle = \biggl[\frac{1}{2}\sinh\biggl(\frac{2\beta\gamma}{M}\biggl)\biggl]^{\frac{N}{2}}exp\biggl[\frac{1}{2}\log \coth\biggl(\frac{\beta\Gamma}{M}\sum_{i=1}^N\hat{\sigma}_{i, k}'\hat{\sigma}_{i, k+1}\biggl)\biggl] \\
\end{eqnarray}


よって,分配関数 Zは以下のように表すことができる.

 \displaystyle
Z = \lim_{M\rightarrow\infty}\biggl[\frac{1}{2}\sinh\biggl(\frac{2\beta\Gamma}{M}\biggl)\biggl]^{\frac{NM}{2}}\sum_{\{\hat{\sigma}_{i, k}=\pm1\}}exp\biggl\{\sum_{k=1}^M\biggl[\sum_{i, j=1}^N\frac{\beta J_{ij}}{M}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}-\sum_{i=1}^N\frac{\beta h_{i}}{M}\hat{\sigma}_{i, k} \\
+\frac{1}{2}\ln \coth\biggl(\frac{\beta\Gamma}{M}\sum_{i=1}^N\hat{\sigma}_{i, k}\hat{\sigma}_{i, k+1}\biggl)\biggl] \\


この結果から横磁場イジングモデルの分配関数は対応する磁場なしイジングモデルに1次元トロッタ方向に次元を追加したもので表されることがわかる.対応するハミルトニアンは以下の式で表される.

 \displaystyle
H_{eff} = -\sum_{i, j=1}^N\sum_{k=1}^M\frac{J_{ij}}{M}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}-\sum_{i=1}^N\sum_{k=1}^M\frac{h_{i}}{M}\hat{\sigma}_{i, k}-\frac{1}{2\beta}\ln \coth\biggl(\frac{\beta\Gamma}{M}\biggl)\sum_{i=1}^N\sum_{k=1}^M\hat{\sigma}_{i, k}'\hat{\sigma}_{i, k+1} \\


シミュレーテッドアニーリングとの対比

これまでの計算から,横磁場イジングモデルの分配関数は対応する古典イジングモデルの分配関数に対して虚時間方向に次元を一つ追加したもので表せることがわかった. このことから,量子アニーリングのシミュレータは,一見単にトロッタ数 Mで並列化したシミュレーテッドアニーリングと同じもののように感じるかもしれないが,上記で計算したハミルトニアンを読み解いてみると,そうでないことがわかる.
ハミルトニアン H_{eff}の第1・2項を見ると,同じ層内( kが同一)での相互作用のみであるが,横磁場の効果を表す第三項では,層間の相互作用を持っている.( k k+1の相互作用)このことから時間発展とともに層同士が互いに独立ではなく,相互に作用し合うため,単にトロッタ数分の並列処理をしているわけではないようだ.
どのような相互作用でそれがどう基底状態探索に効いているのかなどについてはまた別途追っていきたい.

計算の流れ

計算の全体的な流れは以下のようになる.

  1. パラメータの初期化
    ・温度  T
    ・横磁場の強さ  \Gamma

  2. スピンの初期化  M層のそれぞれ N個のスピンをランダムに初期化する.

  3. ランダムに層と層中のスピンを選択し,スピンの向きを反転させる.

  4. スピン反転前後のエネルギーの差分を \Delta E とすると,確率 min(1, e^\frac{\Delta E}{T}) でスピンの向きの反転を受け入れる(受け入れない場合は棄却して変更なし)

  5. 任意の回数3,4のステップを繰り返す(モンテカルロステップ)

  6. 横磁場の強さを  \Gamma小さくする.

  7. 任意の回数5,6のステップを繰り返す(アニーリングステップ)

  8. 各層のエネルギーを求め(横磁場の項は除く)一番エネルギーが低い層のスピン配位を解とする

シミュレータを触ってみる

GitHub上ではいくつかの既存実装が存在したが,今回は以下のものを触ってみることにした.

github.com

必要なライブラリ環境をそろえて,リポジトリのトップディレクトリで以下のように実行すると,ヘルプが見れる.ちなみにPython 2系で実装されている.

$ python solve.py qmc -h
usage: Ising Solvers qmc [-h] [-problem PROBLEM] [-steps STEPS] [-dump DUMP]
                         [-spins SPINS] [-P P] [-T T] [-G0 G0] [-Gf GF]
                         [-e0 E0] [-ef EF]

optional arguments:
  -h, --help        show this help message and exit
  -problem PROBLEM  a file containing problem data in three column format, see
                    README for details
  -steps STEPS      number of monte carlo steps (MCS) used, good values depend
                    on problem size
  -dump DUMP        a file to dump move acceptances to, if desired -- slows
                    down performance
  -spins SPINS      initial spin configuration, if desired
  -P P              number of trotter slices (replicas), 20 to 80 is typical
  -T T              ambient temperature in K, 0.015 is typical
  -G0 G0            initial transverse field strength
  -Gf GF            final tranverse field strength, typically near zero
  -e0 E0            initial coupling pre-factor
  -ef EF            final coupling pre-factor

多くの引数があるが,基本的には指定しなくてもデフォルト値が入るようになっていて,必要であれば実行者が上書きできるようにしてある.
マストで渡さないといけない引数は-problem であり,これは,以下のような形式のファイルのパスを渡す必要がある.(sample_dataディレクトリ配下にサンプルが置いてある)

$ head -5 sample_data/ising12.txt
random 12x12 spin ising model with solution energy of -18,972,276
1 2 -28420
2 3 17885
3 4 -146311
4 5 128784

2行以降は以下の形式で

i j J_ij

すなわち,スピン間の相互作用を示す J_{ij}を設定するためのファイルであることがわかる.今回の最適化では実用的な問題が設定してあるのではなく,単にイジングモデルのスピン間にファイルで設定した相互作用を与え,そのイジングモデル基底状態を探索するのみである. また, -dump 引数を指定しておくと,実行の過程がファイル出力される.-steps でアニーリングステップ数が指定できる.
以下のように実行すると,用いたパラメータ等の出力が返ってくる.

$ python solve.py qmc -problem sample_data/ising12.txt -dump output.txt -steps 1000

Ising Spin Glass
+----------------+-------------------------------------------------------------------+
| Field          | Value                                                             |
+----------------+-------------------------------------------------------------------+
| data file      | sample_data/ising12.txt                                           |
| description    | random 12x12 spin ising model with solution energy of -18,972,276 |
| initial config | 0xe696019ff2add3594c219fcbfd9aef9b0dd4                            |
| current config | 0x313fac59f68f16d3ec6a17d6aca964b2025d                            |
| initial energy | 1627456.0000000007                                                |
| current energy | -15899168.000000006                                               |
+----------------+-------------------------------------------------------------------+

Solver: Path-Integral Quantum Monte Carlo
+-----------+------------+
| Parameter | Value      |
+-----------+------------+
| G0        | 3          |
| Gf        | 1e-06      |
| P         | 40         |
| T         | 0.015      |
| dump      | output.txt |
| e0        | 1e-06      |
| ef        | 4          |
| steps     | 1000       |
+-----------+------------+
To replicate starting conditions, run with arguments:
    qmc -dump output.txt -G0 3 -ef 4 -P 40 -steps 1000 -T 0.015 -problem sample_data/ising12.txt -Gf 1e-06 -e0 1e-06 -spins 0xe696019ff2add3594c219fcbfd9aef9b0dd4

出力ファイルを見てみると以下のような形式で出力されている.

$ head output.txt
Path-Integral Quantum Monte Carlo
{'dump': 'output.txt', 'G0': 3, 'ef': 4, 'P': 40, 'steps': 1000, 'T': 0.015, 'Gf': 1e-06, 'e0': 1e-06}
<class 'ising.SpinGlass'>(data_file="sample_data/ising12.txt", spin_configuration="0xe696019ff2add3594c219fcbfd9aef9b0dd4")
3,16.27456,0xe696019ff2add3594c219fcbfd9aef9b0dd4,-1
3.0,15.64962,0xe696019ff2add3594c219fcbfd9aef9b0cd4,7
3.0,13.09112,0xe696019ff2add3594c219fcbfd9aef9b0ddc,28
3.0,13.02468,0xe69e019ff2add3594c219fcbfd9aef9b0dd4,12
3.0,16.31304,0xe696019ff2add3194c219fcbfd9aef9b0dd4,30
3.0,15.9128,0xe696019ff2add3594c219fcbfd9aefbb0dd4,0
3.0,14.19152,0xe697019ff2add3594c219fcbfd9aef9b0dd4,24

4行目以降の結果は,左から横磁場の強さ,全層の中から最も低いエネルギーの値,スピンの状態を16進数で表したもの,何層目のスピンを反転させたのか.
このライブラリの実装では,1アニーリングステップの中で,それぞれの層の中から1スピンずつ反転させている.すなわちモンテカルロステップ数はトロッタ数と同じである.そのモンテカルロステップごとに1行の出力がある.

ステップを重ねるにつれてイジングモデルが低エネルギー状態に向かっていることを確認するために,横軸にアニーリングステップ,縦軸はエネルギーでグラフを書くと以下のようになる.

f:id:hirotsuru314:20191216173856p:plain

アニーリングステップが進むにつれてエネルギーが低い状態に向かっていることから,基底状態を探索している様子が見て取れる.

まとめ

古典的な計算機で量子アニーリングをシミュレートする代表的な手法として,経路積分量子モンテカルロ法がある.経路積分を用いる過程でSuzuki・Trotter展開によって得られるハミルトニアンの分配関数の近似値から,マルコフ連鎖モンテカルロ法を用いることで,量子アニーリングの挙動をシミュレートすることができる. 動作としては,問題とするイジングモデルをトロッタ数分だけコピーして,それらが互いに相互作用しながら基底状態を探索していく. 今回用いたqmcのような既存実装を用いると,実用的な問題ではないものの手元のPCで簡単に量子アニーリングの挙動をシミュレートすることができる.

参考