Fire Engine

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

Krylov部分空間を導入して特異スペクトル変換による異常検知の処理を高速化した

1年くらい前に特異スペクトル変換法による異常検知ライブラリを作ったんですが、作ったっきり放置していたので、開発当初からやりたかった計算の高速化処理を書きました。
ずっと放置してた割にはちょいちょいGitHubのスターを押してもらえてて、データサイエンスの流行を感じた。自分ももう一回ちゃんと学び直していこうという気になったので、まずは昔書いたやつの拡張からやっていく。

【目次】

特異スペクトル変換とは?

特異スペクトル変換法の特徴については以前のブログに書いているので、ぜひそちらも読んでください。

特異スペクトル変換法の全体像は以下のようになっています。

f:id:hirotsuru314:20171011214314p:plain
出典:上の図は井手剛氏の著書「入門 機械学習による異常検知―Rによる実践ガイド」のP200 図7.4を元に作成しました。

図のように過去と今のパターンを行列として取り出し、それぞれの特徴抽出をした行列(上図の主部分空間の部分)を生成します。生成した二つの行列の食い違いの大きさを変化度として異常検知を行います。

特異スペクトル変換は、特定の確率分布を仮定していないため、多様な変化に頑強に対応できる点やパラメータチューニングが容易などといった特徴が挙げられます。

一方、特異スペクトル変換の一つの課題として計算コストが高いことが挙げられます。先ほど述べた「行列からの特徴抽出」の部分で、特異値分解と呼ばれる行列計算を用いており、これがかなり計算コストが高いです。したがって、検知を行いたい時間間隔や、パラメータの最適値によっては、リアルタイム計算が厳しくなってきます。

今回は、Krylov部分空間という概念を導入することにより、特異スペクトル変換による異常検知の処理を高速化することが目的です。

Krylov部分空間の導入

理論的な詳細は下記の元論文をご参照ください。この論文は今回実装した異常検知のアルゴリズムについてはもちろんのこと、部分空間同士の距離について丁寧に書かれていて数学的にも大変勉強になりました。

Change-Point Detection using Krylov Subspace Learning

Banpeiの既存実装だと、データから取り出した密行列に特異値分解かけるので、その行列のサイズが大きくなると処理はどんどん重くなりますが、論文のアルゴリズムを用いると、次元数k(通常2~5程度)のKrylov部分空間を導入して、最終的にはk×kの固有値問題に帰着します。

上の論文や他の参考情報(最後に載せている)を元に高速化処理を書きました。

github.com

検証結果

それでは、特異値分解(SVD)を使った既存実装と、今回実装したLanczos法を使った新規実装を比較していきます。(Lanczos法はKrylov部分空間法の計算手法の一つ)

用いるデータは、Banpeiのリポジトリに置いている周波数異常データです。
まず以下のようにしてデータを読み込みます。(ライブラリのインストール方法についてはREADME参照)

import banpei
import pandas as pd

model   = banpei.SST(w=100)
raw_data = pd.read_csv('./tests/test_data/periodic_wave.csv')
data = raw_data['y']

Jupyter Notebookを使うと以下のように%timeというマジックコマンドでプログラムの実行時間を簡単に計測できます。

# 既存実装(SVD)
%time results_svd = model.detect(data, is_lanczos=False)
#出力結果=> 
CPU times: user 4.22 s, sys: 80.4 ms, total: 4.3 s
Wall time: 2.6 s

# 新規実装(Lanczos)
%time results_lanczos = model.detect(data, is_lanczos=True)
#出力結果=> 
CPU times: user 944 ms, sys: 51.2 ms, total: 995 ms
Wall time: 872 ms

用いたサンプルデータ、既存実装により計算した変化度スコア、新規実装により計算した変化度スコアを描くと以下のようになりました。

f:id:hirotsuru314:20190203093607p:plain

実行時間に関しては、4.3sから995msとなり、約4.3倍早くなりました。
この結果は完全に取り出した行列のサイズに依存しています。今回の例で用いた行列のサイズは50×100ですが、このサイズが大きくなればなるほど、新規実装の方が顕著に処理が早くなります。

次に変化度スコアのグラフについてですが、一番下の新規実装(Lanczos)は既存実装(SVD)に比べてグラフがノイジーで粗くなっているように見えます。この理由は詳細には調べていませんが、異常検知の観点でいくと、しっかり変化度スコアが立ち上がってくれているなら、多少のグラフの粗さは問題ないかと思っています。

さいごに

難解な数学的な処理をプログラムに落とすと謎の達成感があるんだけど、書いたものは使ってなんぼなので、その辺りは今後の課題としてやっていきます。
データサイエンス周りは1年間くらい追ってなかったので、基礎的なインプットも再開してやっていく!

参考