はじめに
こんにちは、Misocaの洋食（@yoshoku）です。仕事はデータ分析です。
この記事は、Misoca Advent Calendar 2017の23日目の記事です。
今回は、データ分析なエンジニアらしく、機械学習アルゴリズムの実装に挑戦します。また、MisocaはRubyな会社なので、Rubyの線形代数ライブラリであるNumo::NArrayを使います。Numo::NArrayは、PythonでいうところのNumpyです。だいたい同じことができます。
実装する機械学習アルゴリズムは、ICML'17で発表されたCoherence Pursuit [1] です。シンプルなアルゴリズムで、外れ値に強い主成分分析を実現します。深層学習が流行っていますが、ベーシックな機械学習アルゴリズムも、忘れずに押さえていきたいところです。
Coherence Pursuit
まずは、本稿で実装するCoherence Pursuit（CoP, 一貫性追跡法）について、ざっくりと簡単に解説します。CoPは、ロバスト主成分分析の新しいアルゴリズムとして提案されています。主成分分析（Principal Component Analysis, PCA）は、（複数の解釈がありますが）もとのデータと、部分空間に射影してそれを復元したデータとの誤差が、小さくなるような直交射影を求めます。一般に、PCAは、外れ値の影響を受けやすいと言われています。この問題を解決するために、様々な外れ値に対して頑健（Robust）なPCAが提案されています。Robust PCAの多くが、反復アルゴリズムを採用しています。これに対して、CoPでは、反復アルゴリズムを使用しないため、高速でシンプルだと主張されています。
問題設定
$m$次元、$n$点のデータ$D\in\mathbb{R}^{m\times n}$が与えられたとします。CoPでは、column-wise modelというものを採用し、データ$D$を、外れ値ではないデータ$A\in\mathbb{R}^{m\times n_{1}}$と、外れ値データ$B\in\mathbb{R}^{m\times n_{2}}$、置換行列$T$を用いて次の様に表します。
D = [A\hspace{1em}B] T
個々のデータ（$D$の列ベクトル）には、外れ値でないもの（$A$の列ベクトル）と、外れ値であるもの（$B$の列ベクトル）があり、それらが順不同に混在して（$T$で並び替えられて）、データ$D$を構成していると考えます。CoPは、この$A$にある$r$次元の部分空間を求めることを考えます。
アルゴリズム
CoPでは、まず、相互一貫性（mutual coherence）というものを定義します。相互一貫性は、コサイン類似度の絶対値で表されます。二つの非ゼロベクトル$\mathbf{v} _{1}\in\mathbb{R}^{m}$と$\mathbf{v} _{2}\in\mathbb{R}^{m}$の相互一貫性は、次の様に定義されます。
\frac{|\mathbf{v}_{1}^{\top}\mathbf{v}_{2}|}{||\mathbf{v}_{1}||\hspace{0.1em}||\mathbf{v}_{2}||}
CoPは、外れ値データは他のデータと強い類似性を持たないと考えます。たしかに、外れ値データは、特徴空間上で、多くのデータが集まるところから離れて、孤立して存在していると想像できます。この場合、外れ値データは、他のデータとは距離が大きく離れる（強い類似性を持たない）ことになります。
CoPのアルゴリズムは、全てのデータ間で相互一貫性を計算し、他のデータとの相互一貫性の値が大きくなるようなデータで張られる空間を部分空間として得ます。論文では、この相互一貫性にもとづくデータの選択方法により、２種類のアルゴリズムが提案されています。本稿では、シンプルな方法を解説します。とても簡単で、相互一貫性が大きくなるデータを$k$個選択し、その$k$個のデータを$A$と考え、通常のPCA（もしくは特異値分解）を適用することで、部分空間への射影行列$U$を得るというものです。以下にまとめます。
- 各データを$\ell_{2}$ノルムにより正規化したもの$\mathbf{x}_ {i}=\mathbf{d}_ {i}/||\mathbf{d}_ {i}||_ {2}$からなる行列$X\in\mathbb{R}^{m\times n}$を用意する。
- 行列$G=X^{\top}X-I$を用意し、その列ベクトルのノルム$p_ {i}=||\mathbf{g}_ {i}||$を要素とするベクトル$\mathbf{p}\in\mathbb{R}^{n\times 1}$を用意する。列ベクトルのノルムは、$\ell_ {1}$もしくは$\ell_ {2}$を選択する。
- ベクトル$\mathbf{p}$の要素で大きいものを$n_{1}$個を選択し、対応する行列$X$の列ベクトルから行列$A$を構成する。
- 行列$A$に対して通常のPCAを適用することで、射影行列$U\in\mathbb{R}^{m\times r}$を得る。
部分空間の次元数$r$のほかに、外れ値ではないサンプルの数$n_{1}$が、パラメータとなります。これらパラメータは、一般にはデータが与えられた時点では未知であるため、調整が必要です。論文では、$n_{1}$を与えずに適応的にサンプリングする方法も提案されています。しかし、適応サンプリングに関して調整すべきパラメータはあるので、どちらが良いとは言い切れないところがありあす。
実装
Numo::NArrayを使って、CoPを実装していきましょう。ここで、機械学習ライブラリの文化？にあわせて、論文では列ベクトルがデータとなっていましたが、実装では行ベクトルをデータとします（ちょうど転置した状態）。
まずは、Numo::NArrayをインストールします。特異値分解したいので、Numo::Linalgもあわせてインストールします。
$ gem install numo-narray numo-linalg
コードは次の様になります。Numo::NArrayのお陰でシンプルに書けます。
require 'numo/narray'
require 'numo/linalg'
# 余談ですが、私は自前でビルドしたOpenBLASを/usr/local/libに配置しているので、
# 次の様にパスを指定する形で、Numo::Linalgをrequireしてます。
#
# require 'numo/linalg/linalg'
# Numo::Linalg::Loader.load_openblas('/usr/local/lib')
# 一貫性追跡法による射影行列を計算する
# @param samples [Numo::DFloat] (shape: [n_samples, n_dimensions]) もとのデータ
# @param n_inliers [Integer] 外れ値でないデータの数
# @param n_subspace_dimensions [Integer] 部分空間の次元数
# @return [Numo::DFloat] (shape: [n_dimensions, n_subspace_dimensions]) 部分空間への射影行列
def coherence_pursuit(samples, n_inliers, n_subspace_dimensions)
# もとのデータのサンプル数と次元数を得る。
n_samples, n_dimensions = samples.shape
# L2ノルムで正規化する。
norms = Numo::NMath.sqrt((samples**2).sum(1))
normalized_samples = samples / norms.tile(n_dimensions, 1).transpose
# 相互一貫性によるベクトル（論文中のベクトルp）を計算する。
gram_mat = normalized_samples.dot(normalized_samples.transpose) - Numo::DFloat.eye(n_samples)
coherences = gram_mat.abs.sum(1)
# 外れ値ではないデータ（inlier）を取り出す。
inlier_ids = coherences.sort_index.reverse[0...n_inliers]
inliers = samples[inlier_ids, true]
# 射影行列を求める。
_s, _u, vt = Numo::Linalg.svd(inliers)
vt[true, 0...n_subspace_dimensions]
end
実験
LIBSVM Dataにある手書き数字データのUSPSを使用して、PCAとCoPの復元誤差を比較したいと思います。手書き数字の「1」のデータに、外れ値として「8」のデータを加えて、「1」の部分空間を得られるかを試みます。LIBSVM形式のデータをNumo::NArrayの行列で読み込むにはSVMKitを使用します。
実験データの準備は次のようになります。
require 'svmkit'
# USPSデータを読み込む。
# （ラベルが1〜10で割り当てられているので数字画像と一致させるために1を引きます）
base_samples, labels = SVMKit::Dataset.load_libsvm_file('usps.t')
labels -= 1
# 数字「1」に対応するデータを取り出す。
one_ids = labels.to_a.map.with_index { |val, idx| idx if val == 1 }.compact
one_samples = base_samples[one_ids, true]
# 同様にして数字「8」に対応するデータを取り出す。
eight_ids = labels.to_a.map.with_index { |val, idx| idx if val == 8 }.compact
eight_samples = base_samples[eight_ids, true]
# 数字「1」のデータに、外れ値として「8」のデータを少し足して、実験データを作る。
samples = Numo::NArray.concatenate([one_samples, eight_samples[0...10, true]], axis: 0)
実験データが用意できたので、PCAとCoPそれぞれで射影行列を求め、復元誤差を比較してみましょう。ちなみに、部分空間の次元数は適当に16としました。特別な意図はありません。CoPのパラメータである、外れ値でないデータの数も適当に100としてみました。後で調べたところ、数字「1」のデータは264個でした。
def pca(samples, n_subspace_dimensions)
_s, _u, vt = Numo::Linalg.svd(samples)
vt[true,0...n_subspace_dimensions]
end
# PCAとCoPそれぞれで射影行列を求める。
pca_proj_mat = pca(samples, 16)
cop_proj_mat = coherence_pursuit(samples, 100, 16)
# PCAとCoPそれぞれで、数字「1」のデータでの復元誤差を計算する。
pca_diff = one_samples - (one_samples.dot(pca_proj_mat)).dot(pca_proj_mat.transpose)
puts "Reconstruction Error (PCA): %f" % Numo::Linalg.norm(pca_diff)
cop_diff = one_samples - (one_samples.dot(cop_proj_mat)).dot(cop_proj_mat.transpose)
puts "Reconstruction Error (CoP): %f" % Numo::Linalg.norm(cop_diff)
これを実行すると次の様になりました。CoPの方が復元誤差が小さく、PCAと比較して、外れ値でないデータである数字「1」のデータの部分空間を推定できていることがわかります。
Reconstruction Error (PCA): 223.949304
Reconstruction Error (CoP): 174.558816
おわりに
Numo::NArrayを使うと、Rubyでも機械学習アルゴリズムが簡単に実装できますぞ〜!!お試しあれ〜!!