Generalized Mahalanobis depth in the reproducing kernel Hilbert spaceをPythonで実装した話

はじめに

2011年に以下の論文が出版されている．

Yonggang Hu, Yong Wang, Yi Wu, Qiang Li & Chenping Hou, "Generalized Mahalanobis depth in the reproducing kernel Hilbert space," Statistical Papers volume 52, pages 511-522 (2011). link.springer.com

この論文において，再生核ヒルベルト空間（Reproducing Kernel Hilbert Space; RKHS）上でMahalanobis depthを計算する方法が提案された．それを具体的にPythonで実装してみたので紹介するのが本記事の主旨である．

統計的depthとは

統計的depthは多次元空間のmedianを考える際に出てくる統計量である．多数のサンプルデータ点（一般には多次元ベクトルの集合）からなるdata cloudが与えられたとき， depthにより各点がdata cloudの「中心」からどれくらい近いか・離れているか (= centrality)が定量化される．値は0から1の範囲を取り，最大値を与える点をmedianと呼ぶ（ことにする）．たとえ多次元空間であってもdepthの計算手段/アルゴリズムを与えさえすれば，medianが定義できるというわけである．

depthにより，多次元空間における各点ごとに一次元の順序 (order) が導入される．このときの順序は「中心から外向き (center-outward)」である．当然ながらその「中心」にはmedianがいる．順序を用いることで多次元データ群上で分位数（quantile）と階数（rank）に相当するものが計算可能となり，統計的depthは統計学分野に幅広い応用を持つ．

depth を計算するための手段（depth関数）も様々なものが提案されており，具体的には

Halfspace depth
Mahalanobis depth
Simplicial depth
Simplicial volume depth
Zonoid depth
Lens depth
Spatial depth
Projection depth

などがある（各depthの出典は省略）．

Mahalanobis depthとは

ある多次元空間のデータ点 $\mathbf{x}$ に対して， data cloud $X$ から計算されるマハラノビス距離を $d^{2} (\mathbf{x}, X)$ と書いたとき，Mahalanobis depth (MHD) は次式で計算できる．

$$ MHD (\mathbf{x} \mid X) = \frac{1}{1 + d^{2}(\mathbf{x}, X)} $$

論文とnotationが全然違うが許してほしい（ただのサボり！）．

$x \in [0, \infty) $ なる $x$に対して

$$ y = \frac{1}{1+x} $$

のグラフを思い浮かべれば，depthの値の範囲が $(0, 1]$ となることはすぐに分かる．

MHDはマハラノビス距離の逆数もどきなので，その等高線は楕円（一般には超楕円体）に限られる．データを生成する（真の）分布が楕円対称ならばうまくdepthもフィットするが，一般のサンプルデータ群は複雑な形状をしているため，必ずしも最適とは限らない．さらにマハラノビス距離を計算する際の分散共分散行列の正定値性を仮定している点も制約になる．一般には必ずしも正定値にならず特異 (singular)な行列となり，逆行列が計算できない，つまりマハラノビス距離が計算できなくなる．

Generalized Mahalanobis depth およびその実装

論文ではまず上記MHDの欠点を克服するための Generalized Mahalanobis depth (GMHD) が導入される．さらにGMHDの具体的な実現方法として，再生核カーネル関数に対応するRKHS上でGMHDを計算する式が定理として与えられる．具体的にはkernel PCAに相当する固有値問題を解き，固有値と固有ベクトルから計算できる．論文中ではこのdepthは特にkernel mapped GMHD (kmGMHD)と呼ばれている．

kmGMHDの実装は簡単である．scikit-learnベースで話をすると，

学習データに対してkernel PCAをfitし，データ全体をtransformする．手元には変換された特徴ベクトルが手に入る（RKHSにおける座標成分に相当）．
非ゼロの固有値に対応する特徴ベクトル群を取り，固有値の平方根の逆数によりベクトルの各成分をスケーリングする．
特徴ベクトルの各成分の二乗和を計算する．この量を仮にrと置く．
depth = 1 / (1 + r) により計算される．

以下のリンクより，実際に実装したクラスを確認することができる．

上記の説明に対応しているのは以下のコードたちである（抜粋）．

x_transformed = self.kpca.fit_transform(X)  # [n_samples, n_components]
non_zeros = np.flatnonzero(self.kpca.eigenvalues_)
x_transformed = x_transformed[:, non_zeros]
x_transformed = x_transformed / np.sqrt(self.kpca.eigenvalues_[non_zeros])
self.depth_values = 1 / (1 + np.sum(x_transformed**2, axis=1))

実験結果

depth のデモは下記のnotebookを参考に．

これらのnotebookにおいて，各種のtoy datasetに対してdepthを計算した．Figure 1にMHDの結果を，Figure 2 にkmGMHDの結果をそれぞれ示す． Figure 1より，MHDの等高線は楕円型であり，data cloudに必ずしもきれいにフィットしていないことが分かる．