2022-05-23

論文解説 Group Equivariant Convolutional Networks

データサイエンス論文

前回の記事では、Equivariant Neural Networksというデータの対称性に着目した深層学習の設計の新しいパラダイムについて概観した。

blog.tsurubee.tech

今回は、2016年に登場したEquivariant Neural Networksの先駆け的な存在であるGroup Equivariant Convolution Networksの論文を解説する。

論文紹介
- アブストラクト和訳
群論の予備知識
従来のCNNの同変性
Group equivariant Convolutional Neural Networks（G-CNNs）
評価結果
- Rotated MNIST
- CIFAR-10
まとめ

論文紹介

タイトル：Group Equivariant Convolutional Networks（ICML 2016）

著者：Taco S. Cohen, Max Welling

アブスト ラクト和訳

我々は、対称性を利用することでサンプルの複雑さを軽減する畳み込みニューラルネットワークの自然な一般化であるGroup Equivariant Convolutional Neural Networks（G-CNNs）を紹介する。G-CNNsはG-convolutionsという新しいタイプの層を用い、通常の畳み込み層よりも実質的に高度な重みの共有を享受できる。G-convolutionsは、パラメータ数を増やすことなく、ネットワークの表現力を向上させる。G-convolutionsは使いやすく、並進、反射、回転によって生成される離散的な群に対して、無視できるほどの計算量で実装することが可能である。G-CNNはCIFAR10と回転したMNISTでSOTAの性能を達成した。

群論の予備知識

ここでは、論文を理解するために必要な群論の知識を簡単に説明する。

群の定義

空でない集合 $G$ 上に二項演算 $G \times G \rightarrow G$ が定義されていて、次の性質を満たすとき、 $G$ を群という。

(1) 結合法則

　任意の $a, b, c \in G$ に対して、 $(ab)c=a(bc)$ が成り立つ。

(2) 単位元の存在

　ある $e \in G$ が存在し、任意の $a \in G$ に対して、 $ae=ea=a$ を満たす。
　元 $e$ を単位元という。

(3) 逆元の存在

　任意の $a \in G$ に対して、 $ab=ba=e$ を満たす $b \in G$ が存在する。
　元 $b$ を逆元という。

群の具体例

1. 並進群（Translation group）

群の簡単な例として、2次元整数の並進の集合である $\mathbb{Z}^2$ がある。これは、元として2次元の座標 $(n, m)$ 、演算として単純な足し算を考える。そうすると、結合法則を満たすことは自明であり、単位元は $(0, 0)$ 、逆元は $(-n, -m)$ となることから、前述の群の定義を満たすことが確認できる。

2. p4m群

任意の並進移動と、任意の回転中心における90度の回転操作を元とした集合を考える場合、集合は、並進移動・回転操作の合成に対して群となる。また、p4群に鏡映操作を加えた群をp4m群と呼び、この二つの群が、論文中のG-CNNsで対象としている群である。

p4m群の演算は、4つの整数 $m, r, u, v$ を用いて以下のように行列として表現できる。

ここで、 $m \in \left\{0, 1\right\}$ 、 $0 \leq r \lt 4$ 、 $(u, v) \in \mathbb{Z}^2$ である。

準同型写像

$f$ を群 $G$ から群 $G'$ への写像とする。任意の $a, b \in G$ に対して、

$f(ab) = f(a)f(b) \tag{1}$

を満たす写像 $f$ を準同型写像という。ここで、 $ab$ は群 $G$ 上の演算であるが、 $f(a)f(b)$ は群 $G'$ 上の演算であることに注意されたい。準同型写像は、群の演算の性質が引き継がれていることを示しており、このことを「群の構造を保つ」という。

群の作用

群の本質は、それがある対象に「作用する」ことである。

（線形代数と群の表現Ⅰの「はじめに」より）

とあるように、作用というのは群論において重要な概念である。
$G$ を群、 $X$ を集合とすると、 $G$ が $X$ へ（左から）作用するとは、写像 $L: G \times X \ni (g, x) \longmapsto L(g, x) = L_g(x) \in X$ であり、次の二つの性質を満たすものである。

(1) $L(e)=I_X$

　群 $G$ の単位元に対応する $X$ 上の置換は $X$ 上の恒等変換である

(2) $L(g_1g_2)=L(g_1)L(g_2)$ 　 $(g_1,g_2\in G)$

　群 $G$ の二つの元の積 $g_1 g_2$ に対応する $X$ 上の置換は、 $g_1$ および $g_2$ にそれぞれ対応する置換の合成（合成写像）である。

(2)の条件は、式(1)の準同型写像と同じ形をしている。すなわち、群 $G$ が集合 $X$ に作用するとは， $G$ の各元 $g$ に対して、準同型写像の定理を満たす $X$ 上の変換 $L(g)$ が与えられているときをいう。

同変性

入力 $x$ について以下の式を満たすとき、「変換 $\Phi$ は操作 $T_g$ に対して同変である」という。同変性のより詳細な説明は前回の記事を参照。

$\displaystyle \Phi(T_gx)=T'_g\Phi(x) \tag{2}$

従来のCNNの同変性

CNNの畳み込み層から得られる特徴マップ $f$ は、各ピクセルの座標 $(p,q) \in \mathbb{Z}^2$ において、 $K$ 次元のベクトルを返す関数 $f: \mathbb{Z}^2 \rightarrow \mathbb{R}^K$ と表現できる。 $K$ はチャネル数を示している。
ここで、特徴マップ $f$ にある群 $G$ の元 $g$ を作用させることを考えると、下式のように表せる。

$\lbrack L_g f\rbrack(x) = f(g^{-1}x) \tag{3}$

ここでなぜ $g$ の逆元が登場するかなど関数への作用については、線形代数と群の表現Ⅰの「8.2 群の関数への作用」の解説が直感的でわかりやすい。ここでは、操作を $g$ として抽象化して書いているが、例えば、 $g$ が並進 $t = (u, v) \in \mathbb{Z}^2$ を表現するのであれば、 $g^{-1}x$ は単純に $x - t$ を意味する。

ここから、従来のCNNが並進移動に対して同変性を持つことを示す。フィルタ $\psi$ を用いた畳み込み処理（正確には相互相関であるが、以降慣例に従い畳み込みと呼ぶ）は以下の式で表される。

$\displaystyle \lbrack f \star \psi \rbrack(x)=\sum_{y\in \mathbb{Z}^2} \sum_k f_k(y)\psi_k(y-x) \tag{4}$

$k$ はチャネルのインデックスを示している。ここで、一つのチャネルに着目すると、以下の式変形から並進移動の後に畳み込みをしたものは、畳み込みをした後に並進移動をしたものと同じであることがわかる。

$\displaystyle \begin{eqnarray} \lbrack \lbrack L_tf \rbrack \star \psi \rbrack(x) & = & \sum_{y\in \mathbb{Z}^2} f(y-t)\psi(y-x) \\ & = & \sum_{y\in \mathbb{Z}^2} f(y)\psi(y+t-x) \\ & = & \sum_{y\in \mathbb{Z}^2} f(y)\psi(y-(x-t)) \\ & = & \lbrack L_t \lbrack f \star \psi \rbrack \rbrack(x) \\ \end{eqnarray}$

2行目は、 $y \rightarrow y+t$ を代入している。この結果は、畳み込み処理は並進群に対して同変であることを意味している。また、従来のCNNにおける特徴マップ $f$ は、群 $\mathbb{Z}^2$ 上の関数と捉えることができる。 G-CNNsではこれを任意の群に対して一般化する。

Group equivariant Convolutional Neural Networks（G-CNNs）

G-CNNsとは

前述のように従来のCNNの畳み込み層はデータの並進移動に対する同変性を持つが、下図のように回転操作に対しては同変性を持たない。つまり、入力画像を畳み込み処理して得られた特徴マップを左上に並進移動させたものは、入力画像を同様に左上に並進移動した画像を畳み込み処理して得られた特徴マップと一致するが、同様のことは回転操作に対しては成り立たない。

これは、例に挙げた手書き数字のように画像内の認識対象に上下左右の明確な向きがある場合には問題になりにくいが、同じ画像でも例えば衛星画像や顕微鏡画像など明確な向きが存在しないものも多く、この場合は回転しても画像の意味は変わらないため、回転に対する同変性が望まれる。 Group equivariant Convolutional Neural Networks（G-CNNs）では並進移動に加えて、鏡映および90度の倍数の回転操作に対する同変性を持つように一般化した畳み込み層（G-convolutions）を実現した。

群上の関数

G-CNNsにおける特徴マップは、群 $\mathbb{Z}^2$ 上の関数ではなく、任意の群 $G$ 上の関数であり、特に論文中では群p4およびp4m上の関数として定義されている。この群上の関数は非常にイメージしづらいが、論文のFigure 1では以下のようにp4群上の特徴マップのイメージが描かれている。

出典：『Group Equivariant Convolutional Networks』のFigure 1

図の左のように、90度の倍数の回転に関連する4つの特徴マップを円上に描く。そうすると、この図の各ピクセルは回転座標（ピクセルが現れる特徴マップ）と2つの並進座標（特徴マップ内のピクセル位置）を持っている。p4上の関数に90度回転 $r$ をかけると、各特徴マップはその赤い矢印に従って、同時に90度回転を行う。この操作の結果が図の右に示されている。

G-convolution：群同変な畳み込み層

G-CNNsでは独自に3つの層（G-convolution、G-pooling、nonlinearity）を定義しているが、ここではG-convolutionについてのみ解説する。 G-CNNsのp4群上の畳み込み層の処理の全体像が下のスライドにわかりやすくまとまっている。

出典：AMMI Course "Geometric Deep Learning" - Lecture 8 (Groups & Homogeneous spaces)

ここで重要なポイントとして、G-CNNsの畳み込み層は第1層目と第2層目以降ではその性質が大きく異なる。なぜなら、第1層目のフィルタは $\mathbb{Z}^2$ 上の関数であるが、第2層目以降では $G$ 上の関数になるからである。

第1層目

式(4)の従来のCNNの畳み込み処理は、画像や特徴マップに対してフィルタをシフトさせ、内積を計算する。このフィルタのシフトをある群 $G$ からの一般的な変換に置き換えることで、G-CNNの第1層で用いられるG-correlationを得ることができる。

$\displaystyle \lbrack f \star \psi \rbrack(g)=\sum_{y\in \mathbb{Z}^2} \sum_k f_k(y)\psi_k(g^{-1}y) \tag{5}$

上式の入力画像 $f$ とフィルタ $\psi$ はともに $\mathbb{Z}^2$ 上の関数であるが、特徴マップ $f \star \psi$ は群 $G$ 上の関数である。

第2層目以降

第2層目以降では、入力の特徴マップが群 $G$ 上の関数であるため、フィルタ $\psi$ も群 $G$ 上の関数である必要がある。

$\displaystyle \lbrack f \star \psi \rbrack(g)=\sum_{h\in G} \sum_k f_k(y)\psi_k(g^{-1}h) \tag{6}$

この式が群の作用に対して同変になることの証明は論文の式(12)にある。従来のCNNやG-CNNsの畳み込み処理の式は、並べて比較してみると何が変わっているのかわかりやすい。

「CNN」から「G-CNNs第1層目」では、 $-x$ という並進移動が群の演算 $g^{-1}$ （例えば回転操作）に一般化されている。また「G-CNNs1層目」から「G-CNNs第2層目以降」では、特徴マップとフィルタがともに群 $G$ 上の関数となる。数式上では $\mathbb{Z}^2$ 上の $y$ が $G$ 上の $h$ に置き換えられているのみである。

これらの式を用いて、実際にG-CNNsの畳み込み層がどのように特徴マップを処理しているのかは、下の記事に視覚的に解説されているので非常に参考になる。

medium.com

評価結果

論文の評価では、Rotated MNISTとCIFAR-10を用いた実験を行い、G-CNNsは当時のSOTA性能を達成した。

Rotated MNIST

データセットとして、ランダムに回転させたMNISTを用いた。ベースラインとなる通常のCNNモデルとして、3×3畳み込みの7層（最終層は4×4）、各層20チャネルのCNNを用いた（Z2CNN）。次に、各畳み込みをp4-Convolutionで置き換え、Z2CNNとパラメーター数をほぼ同一にするためにフィルタ数を $\sqrt{4}=2$ で割った。最後の畳み込み層の後に回転に対するMax Poolingを追加した（P4CNN）。また、ネットワークアーキテクチャにおいて、中間層で同変性を持たせることが重要であり、早期に不変性を持つことは望ましくないという仮説の検証のために、P4CNNの各畳み込み層の後に、回転に対するMax Poolingを行ったモデルも比較した（P4CNNRotationPooling）。評価結果は、以下の表の通りである。

『Group Equivariant Convolutional Networks』のTable 1

提案するP4CNNはZ2CNNの誤差を半減させることがわかった。また、各層が回転に対する不変性を持つP4CNNRotationPoolingは、Z2CNNよりも優れているが、P4CNNよりは誤差が大きいことから、中間層で同変性を持たせることが性能向上に寄与していることが示された。

CIFAR-10

モデルとして、p4、p4m、および標準的な $\mathbb{Z}^2$ の畳み込みを、2種類のベースラインアーキテクチャ（All-CNN、ResNet44）で比較した。またData Augmentationの影響を評価するために、CIFAR10のデータを水平反転と小さな平行移動により増強したCIFAR10+も合わせて評価した。評価結果を以下の表に示す。

『Group Equivariant Convolutional Networks』のTable 2

CIFAR10、CIFAR10+のいずれにおいても、ResNet44のp4m-CNNが最も高い性能を示した。また、All-CNNをベースラインアーキテクチャとして用いた場合においても、パラメータ数は同等であるにもかかわらずp4m-CNNが高い性能を示した。

まとめ

今回の記事では、従来の畳み込み層をp4群やp4m群に対して同変性を持つように一般化したG-CNNsの論文について解説した。G-CNNsは、パラメータ数を増やすことなくネットワークの表現力を向上させることができ、回転したMNISTとCIFAR-10で従来のCNNを大きく上回る性能を達成した。G-CNNsの大きな制約としては、対象とする群が離散群に限られることである。また、計算量が群のサイズに比例して大きくなるため、サイズの大きな群に対する適用も難しい。しかし、G-CNNsの提案後に、計算量が群のサイズに依存しない同変なネットワークや連続群に適用可能なネットワークが提案されている。今後はそれらの論文についても見ていきたい。

2022-05-12

Equivariant Neural Networksの概論：群論を用いた深層学習の設計の進展

データサイエンス論文

最近、Equivariant Neural Networksというデータの対称性に着目した深層学習の設計の新しいパラダイムの存在を知り、非常に興味を持っている。Equivariant Neural Networksは、深層学習の更なる汎化性能の向上や学習の効率化、適用分野の拡大などを実現する可能性を秘めている。今回は、Equivariant Neural Networksとは何なのか、これまでの研究事例、どういう分野で応用できるのかなど、最近私が学んだことを概論的にまとめていく。

Equivariant Neural Networksの概要
不変性と同変性
群論との関係
Geometric Deep Learningとの関係
CNNと同変性
様々なEquivariant Neural Networks
応用分野
さいごに
- 更新：2022年5月23日公開
参考文献

Equivariant Neural Networksの概要

ニューラルネットワーク（以降、NN）を用いた深層学習は、対象とする問題やデータに応じて最適な設計を行う必要がある。問題やデータが持つ固有の構造や制約は帰納バイアス（inductive bias）と呼ばれ、学習データ以外の有効な事前知識として深層学習の設計に積極的に活用される傾向にある。Equivariant Neural Networksとは、データに潜む対称性を強力な帰納バイアスとして捉え、それをネットワークの設計に取り入れたNNである。

対称性とは、ある対象に特定の操作を適用しても結果が変わらない性質のことであり、身の回りに多く存在する。例えば、手書き数字のデータセットであるMNISTを用いた画像分類問題を例に挙げて考えると、下のように「7」という数字を上下左右に並進移動させても、回転させても、それが「7」であるという分類の結果には変わりがない。これをデータは並進移動や回転に対して対称であるという。

Equivariant Neural Networksでは、これらの対称性を不変性（Invariance）を一般化した同変性（Equivariance）として表す。不変性と同変性について詳細な定義は後述する。Equivariant Neural Networksはそれ自体が新しいネットワークの設計を生み出すものであるが、既存のNNの設計を対称性を通して統一的に解釈できるという側面もある。例えば、畳み込みニューラルネットワーク（以降、CNN）の畳み込み層は並進移動に対する同変性を持つし、TransformerのSelf-Attention層は置換に対する同変性を持っている。したがって、Equivariant Neural NetworksはNNの設計の新しい指針を与えるだけでなく、これまでのNNのアーキテクチャに対して俯瞰的な視点を与えてくれる。

不変性と同変性

不変性と同変性について数式を用いて説明する前に、それらがどういった性質を持っているかを図を使って簡単に説明する。関数 $f$ が入力として数字の配列を受け取り、その配列の要素の置換に対する不変性と同変性の例を下に示す。

不変性は、入力配列の順序に関わらず同じ値を返す、つまり置換という操作に対して結果が変わらない。これは例えば、最大値の計算や総和の計算が該当する。次に、同変性を見てみると、入力配列の順序に対応して、出力配列の順序が変わっている。このように入力に応じて「同じように変わる性質」が同変性であるというイメージを持っておくとよい。

ここまでで不変性と同変性のイメージは共有できたと思うので数式を用いた説明に移る。下の図のように入力 $x$ に対してある操作 $g$ （例えば、並進移動や回転など）と変換 $\phi$ （例えば、NNの畳み込み層など）を適用することを考える。

以下の式を満たすような $g'$ が存在するとき、「変換 $\phi$ は操作 $g$ に対して同変である」という。

$\displaystyle \phi(g(x))=g'(\phi(x))$

つまり、入力 $x$ に操作 $g$ を適用した結果 $g(x)$ に対して変換 $\phi$ を適用した結果 $\phi(g(x))$ が、入力 $x$ に変換 $\phi$ を適用した結果 $\phi(x)$ に対して $g$ に対応する操作 $g'$ を適用した結果 $g'(\phi(x))$ が一致することを示している。これは、入力に対して操作した後に変換をしても、変換をした後に操作をしても結果は変わらないことを意味している。

また、以下の式を満たすとき、「変換 $\phi$ は操作 $g$ に対して不変である」という。

$\displaystyle \phi(x)=\phi(g(x))$

これは、操作 $g$ を適用してもしなくても結果が変わらないことを意味している。同変の定義における $g'$ を恒等変換とした場合に不変となるため、不変は同変の特殊例であることがわかる。

Equivariant Neural Networksは、ここまでで説明した同変性をNNの設計に活用するわけであるが、なぜ不変性ではなく、同変性を用いるのかについては、Taco Cohen氏の学位論文[1]で以下のように述べられている。

We argue and show empirically that in the context of deep learning it is better to learn equivariant rather than invariant representations, because invariant ones lose information too early on in the network.
(DeepLによる翻訳) 我々は、深層学習の文脈では、不変な表現よりも同変な表現を学習する方が良いことを主張し、経験的に示している。なぜなら、不変な表現はネットワークの早い段階で情報を失いすぎてしまうからである。

つまり、NNの各層の処理が特定の操作に対して不変である場合、例えば画像Aが画像Bに対して、上に並進移動しているとか、時計回りに90度回転しているといった相対的な位置や向きなどの情報が失われてしまう。一方、同変の場合、このような相対的な位置や向きの情報をNNの中間層で保持することができる。

群論との関係

群論は、「対称性を記述する数学」や「対称性をはかる数学」と言われることからデータの対称性に着目したEquivariant Neural Networksとの繋がりがあることがわかる。ここでは、群の正確な定義は割愛するが、群とは特定の条件を満たす集合である（実際には演算とセットで考えるがここでは深く立ち入らない）。前述の不変性と同変性の説明のところで、ある変換が「操作」に対して不変/同変であることの定義を述べたが、この「操作」を元（集合の要素）として、それら元のなす集合を群として考える。そして、Equivariant Neural Networksでは、特定の群の全ての元に対して、NNの各層の処理が同変となるように設計する。

Geometric Deep Learningとの関係

Equivariant Neural Networksについて調査しているときに、対象領域が非常に近く、よく出てくるキーワードとして「Geometric Deep Learning」という言葉がある。2021年に公開されたGeometric Deep Learningの文献[2]では、Geometric Deep Learningを以下のように説明しており、対称性に着目していることがわかる。

We call this geometrisation attempt "Geometric Deep Learning", and true to the spirit of Felix Klein, propose to derive different inductive biases and network architectures implementing them from first principles of symmetry and invariance.
(DeepLによる翻訳) 我々はこの幾何学的な試みを「Geometric Deep Learning」と呼び、Felix Kleinの精神に忠実に、対称性と不変性の第一原理から様々な帰納的バイアスとそれを実装したネットワークアーキテクチャを導出することを提案する。

また、文献の5.2節「Group-equivariant CNNs」では、群同変なCNNについて解説されている。このように両者は重なる部分が多い研究領域であるが、Geometric Deep Learningは同変性のみにフォーカスしていないことなどからより対象領域が広いと考えられる。

CNNと同変性

同変性をNNの設計に組み込むことで大きな成功を遂げたモデルの代表的な例として、CNNが挙げられる。以下の図で示したようにCNNの畳み込み層は、並進移動に対する同変性を持っている。図は不変性と同変性の定義の説明の際に用いた図と対応している。

図のように、入力画像「7」を畳み込み処理して得られた特徴マップを左上に並進移動させたものは、入力画像を同様に左上に並進移動した画像を畳み込み処理して得られた特徴マップと一致する。つまり、入力画像を並進移動させても、対応する特徴マップも同様に並進移動する。これは、畳み込み処理が、画像全体で同じ重みを共有したフィルタを用いて、画像に対してフィルタを並進移動させながら内積をとるといった処理を行っているためである。畳み込み層はこの重みの共有により並進移動に対する同変性と全結合層と比べたパラメータ数の削減を実現している。実際にCNNでは、出力層の前にGlobal Average Pooling（GAP）を適用することが多く、これによりモデル全体としては並進に対して不変性となる。

それでは、画像に対して回転操作を行った場合はどうだろうか。結論は、以下の図のようにCNNの畳み込み層は回転に対して同変とならない。

つまり、画像の畳み込み処理により得られた特徴マップを反時計回りに90度回転したものは、画像を反時計回りに90度回転したものに畳み込み処理をして得られた特徴マップと一致しない。これは、例に挙げた手書き数字のように画像内の認識対象に上下左右の明確な向きがある場合には問題になりにくい。しかし、同じ画像でも例えば衛星画像や顕微鏡画像など明確な向きが存在しないものも多く、この場合は回転しても画像の意味は変わらないため、NNの処理が回転に対する同変性や不変性を持つことが望まれる。

この問題を従来のCNNで解決する方法の一つとして、Data Augmentationが挙げられる。つまり、様々な回転操作を施した画像を用意してCNNで学習するという方法である。しかし、この方法では単純にデータ量が膨大になる。また、モデル自体はある画像とそれを回転させた画像の関係性を捉えることはなく、全く異なる画像を同じラベルに紐付けるように学習する。このことからも学習効率が良くないことが想像できる。

様々なEquivariant Neural Networks

Equivariant Neural Networksは2016年頃から非常に活発に研究が行われている。以下のGitHubページには、Equivariant Neural Networksに関連する論文やチュートリアルなどの情報がまとめられている。

github.com

また、E. J. Bekkers氏は、「Group Equivariant Deep Learning」というタイトルの講義をYouTubeで公開しており、既存のEquivariant Neural Networksを数学的な背景から解説してくれている。その講義スライドの中で同変なCNNに関するこれまでの歴史がまとめてある。

出典: Group Equivariant Deep Learning - Lecture 1.5: A Brief History of G-CNNs

ここでは、既存のEquivariant Neural Networksの中からいくつかピックアップして簡単に紹介する。

Group Equivariant Convolution Networks (2016)

Group Equivariant Convolution Networks（G-CNNs）[3]は、Equivariant Neural Networksの先駆け的な存在である。G-CNNsは、並進移動、鏡映、90度の倍数の回転操作を元とするp4m群の作用に対して同変性を持つCNNを実現している。以下の記事は図が豊富で、G-CNNsの処理の内容がわかりやすく描画されている。

medium.com

本ブログ公開後の2022年5月23日にG-CNNsの論文解説をブログにまとめた。

blog.tsurubee.tech

Steerable CNNs (2017)

前述のG-CNNの一つの課題として、計算量が対象とする群のサイズ（元の個数）に比例して大きくなることが挙げられる。Steerable CNNs [4]は、計算量が群のサイズに依存しない同変なCNNを実現している。これにより、同変なCNNをよりサイズが大きな群にスケールすることが可能となる。論文中では、G-CNNと同様にp4mの群の作用に同変なCNNについて説明や評価が行われている。

Spherical CNNs (2018)

従来のCNNは、2次元平面上のデータを取り扱う問題に対して広く適用されているが、Spherical CNNs [5]では、球面上のデータを取り扱う。これは例えば、ドローンなどから取得される全方位ビジョンや地球気象のモデリングのための球面画像などが挙げられる。Spherical CNNsは、3次の特殊直交群 SO(3)の作用に同変な球面上の畳み込み層を提案している。

towardsdatascience.com

B-Spline CNNs (2020)

B-Spline CNNs [6]は、リー代数上で定義されるB-spline基底関数を用いて畳み込みカーネルを拡張することで、同変な畳み込み層を任意のリー群に一般化するフレームワークを提案している。論文中の評価では、病理画像における癌検出および顔のランドマーク推定の二つのタスクにおいて、それぞれ同変な畳み込み層を設計し、どちらも従来のCNNを大きく上回る性能を達成した。

E(n) Equivariant Graph Neural Networks (2021)

E(n) Equivariant Graph Neural Networks (EGNNs) [7]は、E(n)変換に対して同変なグラフニューラルネットワークである。E(n)とは、n次元ユークリッド空間における等長変換群であり、回転、並進、鏡映、置換から構成される。これまでに紹介したモデルと異なり、グラフ構造に着目しているので、創薬などの分子を対象データとする様々なタスクへの応用が期待される。

応用分野

Equivariant Neural Networksはデータに何らかの対称性を見いだせる場面で活用することで、NNの性能向上が期待できる。ここでは、現在の主な応用分野について紹介する。

医療画像

画像内の認識対象が明確な上下左右の向きを持っていない場合、向きが変更される回転に対して同変性・不変性が求められる。このような場面は、医療画像によく現れる。E. J. Bekkersらは、下図のような病理組織画像における有糸分裂の検出（Mitosis detection）、網膜画像における血管のセグメンテーション（Vessel segmentation）、電子顕微鏡画像における細胞境界のセグメンテーション（Cell boundary segmentation）の三つのタスクにおいて、2次元の特殊ユークリッド群SE(2)に対して同変なCNNを設計し、従来のCNNを上回る性能を達成した [8]。

また、P. Müllerらは、MRI画像に3次元の特殊ユークリッド群SE(3)に対して同変なCNNを適用している[9]。下図は、MRI画像における多発性硬化症病変のセグメンテーションの結果を示している。

分子設計・材料設計

分子は3次元構造を持つため、3次元の回転などの操作に対して同変性・不変性が求められる。また、材料科学の分野において，所望の物性を有する材料を開発することは重要な目標の一つであるため、分子やそれをマクロに捉えた材料の物性を深層学習を用いて予測する研究が盛んに行われている。この材料物性の予測に対しても、Equivariant Neural Networksが活用され始めている。 T. Leらは、同変なグラフNNであるEQGATを設計し、低分子の量子力学的特性の予測においてSOTA性能を達成した[10]。

また、生成モデルを用いて、所望の分子を生成するアプローチも研究されている。V. G. Satorrasらは、E(n)変換に対して同変なNormalizing FlowsであるE-NFsを提案し、下図のような3次元の分子構造を生成することを可能にした[11]。

創薬

医薬品として使われる低分子化合物やタンパク質（抗体）も分子であるため、ここで述べる創薬は分子設計に含めても良いが、背景として創薬にフォーカスしているものを切り出して紹介する。 H. Stärkらは、低分子化合物が特定の標的タンパク質にどのように結合するかを予測するために、SE(3)に対して同変なNNであるEQUIBINDを提案した[12]。EQUIBINDは、受容体の結合位置とリガンドの結合姿勢・向きを予測することができる。また、O. E. Ganeaらは、タンパク質複合体の立体構造を予測するために、SE(3)に対して同変なNNであるEQUIDOCKを提案した[13]。EQUIDOCKは、結合ポケットの近似やドッキングポーズの予測を行うことができる。

さいごに

この記事では、Equivariant Neural Networksとは何かから始まり、その応用事例までを概論的に紹介した。Equivariant Neural Networksは、本文中でも紹介した既存のモデルを深掘りして理解しようとしたときに、群論や微分幾何学などの高度な数学の知識が要求される。現在私は、これらの数学を学びながら理解を進めているため、今後はEquivariant Neural Networksのいくつかの既存のモデルについて、学んだことを各論的にまとめた記事を書いていく予定である。

更新：2022年5月23日公開

blog.tsurubee.tech

参考文献

[1] T. Cohen, "Equivariant Convolutional Networks" (2021).
[2] M. M. Bronstein et al., "Geometric Deep Learning: Grids, Groups, Graphs, Geodesics, and Gauges" (2021).
[3] T. Cohen et al., "Group Equivariant Convolutional Networks" (2016).
[4] T. Cohen et al., "Steerable CNNs" (2017).
[5] T. Cohen et al., "Spherical CNNs" (2018).
[6] E. J. Bekkers, "B-Spline CNNs on Lie Groups" (2020).
[7] V. G. Satorras et al., "E(n) Equivariant Graph Neural Networks" (2021).
[8] E. J. Bekkers et al., "Roto-Translation Covariant Convolutional Networks for Medical Image Analysis" (2018).
[9] P. Müller et al., "Rotation-Equivariant Deep Learning for Diffusion MRI" (2021).
[10] T. Le et al., "Equivariant Graph Attention Networks for Molecular Property Prediction" (2022).
[11] V. G. Satorras et al., "E(n) Equivariant Normalizing Flows" (2021).
[12] H. Stärk et al., "EQUIBIND- Geometric Deep Learning for Drug Binding Structure Prediction" (2022).
[13] O. E. Ganea et al., "Independent SE(3)-Equivariant Models for End-to-End Rigid Protein Docking" (2022).

2021-07-21

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

データサイエンス

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

blog.tsurubee.tech

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

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

ShapPackの概要
ShapPackの使い方
ShapPackの詳細
今後の展望
参考文献

ShapPackの概要

github.com

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

マルチプロセスで並列処理できる機能
特性関数を独自に実装して組み込める機能
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_workers，characteristic_func，skip_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の特性関数では，上の特性関数の実装例のinstanceとsubsetsでどちらも関数に渡たされる仕様になっているため，こちらも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に実装できないか検討していきたい．

参考文献

*1:A Unified Approach to Interpreting Model Predictions

*2:Additive Explanations for Anomalies Detected from Multivariate Temporal Data

*3:On Anomaly Interpretation via Shapley Values

*4:Shapley Flow: A Graph-based Approach to Interpreting Model Predictions

*5:Causal Shapley Values: Exploiting Causal Knowledge to Explain Individual Predictions of Complex Models

2021-07-19

協力ゲーム理論のシャープレイ値に基づき機械学習モデルの予測を解釈する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の概要
シャープレイ値
SHAP (SHapley Additive exPlanations)
Kernel 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

参考文献

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

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

2021-06-07

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

著者	鶴田博文, 坪内佑樹
所属	さくらインターネット株式会社さくらインターネット研究所
研究会	第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

2021-01-14

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

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点です．

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

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

２．実装力

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

３．データの捉え方

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

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

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

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

１．Kaggleの解法を論文としてまとめる

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

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

www.sciencedirect.com

２．コンペで入賞して学会で発表する

私が前回参加した「Riiid! Answer Correctness Prediction」では上位入賞者は「AAAI-2021 Workshop on AI Education」で発表できる．
国際学会が主催しているコンペもある．例えば，IEEEが主催した「IEEE-CIS Fraud Detection」で上位入賞者はIEEE CIS Conferenceで発表できる．

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

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

www.nssol.nipponsteel.com

www.jri.co.jp

さいごに

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

2021-01-09

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

データサイエンス Kaggle

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

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

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

Kaggleとの出会い
コンペ参戦準備
コンペに登録
初Submit
特徴量エンジニアリング
モデルの構築
アンサンブル学習
コンペ終了：Leaderboardの公開
コンペ全体を通しての感想と今後

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の全体像を掴むところから始めました．具体的には，以下の本を一通り手を動かしながら読みました．

実践Ｄａｔａ　Ｓｃｉｅｎｃｅシリーズ　ＰｙｔｈｏｎではじめるＫａｇｇｌｅスタートブック (ＫＳ情報科学専門書)

この本は，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時の大きな制約となりました．

Notebook環境のメモリ上限は16GB
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やっていきましょう！