はじめに
今回は、特異スペクトル変換法というアルゴリズムをPythonで実装します。このアルゴリズムは時系列データの異常検知に対して非常に強い力を発揮します。また、ハイパーパラメータ（人が調整する必要のあるパラメータ）が少なく、比較的チューニングが容易であることも特徴の一つです。数学の理論については深追いはせず、アルゴリズムの概要と実装まで書いていきたいと思います。
【目次】
時系列データについて
時系列データとは、時間の推移ととともに観測されるデータのことです。昨今、様々な企業がデータ活用を推進していますが、世の中の実務の現場に貯まっていく多くのデータは時系列のデータです。
データ分析において、時系列データの取り扱いは、他のデータ（クロスセクションデータなどと言われる）に比べて、理論的にも実銭的にも難しいです。例えば、分析をする上で、ヒストグラムを描くということをよくやると思いますが、ヒストグラムを描くということは、暗にデータの相関を無視し、互いに独立であることを仮定しています。 時系列データでは、隣り合う時刻の観測値同士には明らかに相関関係があり、「各観測値が互いに独立である」という仮定は成り立ちません。このようなデータをきちんと取り扱えることは非常に有用です。
【参考】 実務の現場に多い時系列データ分析の際に注意しておきたい点を列挙してみる - 六本木で働くデータサイエンティストのブログ
時系列データの異常と変化点検知
一口に異常と言ってもいろいろありますが、明らかに他の大多数のデータから外れた値を持っている場合は非常に検知が容易です。このような異常を検知することは外れ値検知と呼ばれ、正規分布を用いた単純な外れ値検知の話を以前記事に書きました。
一方、時系列データの異常パターンの特徴としては、各観測データの時刻（横軸）を適当にシャッフルすると検知できなくなるということです。実際に例を見てみましょう。
データの真ん中あたりで波が激しくなって、また戻っているのがわかると思います。変化点という観点で言うと、波の周波数が激しくなっている点（横軸が30のあたり）と、またゆるやかに戻っている点（横軸が60のあたり）が変化点です。このように人間の目でデータを見ると、どこが異常かだいたいわかりますが、コンピュータにこの異常を検知させるのはけっこう難しいです。これを自動で検知してやろうというのが今回の目的です。今回は、1変数の時系列データにおける異常検知の話です。
特異スペクトル変換法の概要
時系列データの異常検知の手法はいろいろありますが、今回は特異スペクトル変換というアルゴリズムを実装してみます。特異スペクトル変換法の概略は下の図に集約されます。
（出典：上の図は井手剛氏の著書「入門 機械学習による異常検知―Rによる実践ガイド」のP200 図7.4を元に作成しました。）
履歴行列とテスト行列
一次元時系列データの分析をする際に、ある任意の幅を設定して、
個の隣接した観測データをまとめて
次元のベクトルとして取り出すということをよくやります。さらに
を時間の経過とともにずらしながら取り出すことで、時系列片をつくっていきます。この取り出す領域のことをスライド窓と呼び、スライド窓により取り出したベクトルを部分時系列と呼びます。
上の概略図からわかるように、特異スペクトル変換法ではこの部分時系列をさらにまとめて行列としてデータを取り出します。この行列自体の特徴が、ある時刻の時系列データの特徴となるわけです。ある時刻のまわりデータを用いて生成した行列をテスト行列と呼び、時刻
より以前のデータを用いて生成した行列を履歴行列と呼びます。特異スペクトル変換法では、この二つの行列の食い違いの大きさを変化度として異常検知を行います。
特異値分解
特異スペクトル変換法の理論の理解には線形代数の知識がある程度必要なので、場合によっては「特異スペクトル変換法は現在と過去の行列を取り出して、その食い違いの大きさをもとに変化度を計算しているのだな」くらいの理解でよいかもしれません。（後述のPythonによる実装まで読み飛ばしていただいてかまいません。）
特異値分解は、固有値分解を長方形行列に拡張したものと見なすことができます。詳細な説明については下記の参考リンクがわかりやすいと思います。固有値分解ってなに？という方にはプログラミングのための線形代数という本が超絶おすすめです。
http://www.misojiro.t.u-tokyo.ac.jp/~murota/lect-kisosuri/singularvalue050129.pdf
数式を追うことも大事ですが、今回やっていることをざっと書くと、前述の履歴行列とテスト行列に特異値分解を行い、特異値上位の左特異ベクトルを取り出し行列としてまとめてます。（行列のランクより少ない数を取り出す） これはマイナーな方向を無視し、ノイズ除去または次元削除をしたことに相当すると考えることができます。また、取り出したそれぞれの列ベクトルで張られる空間を主部分空間などと呼びます。
後に出てくるコードを見ればわかりますが、いかにも難しそうな特異値分解もNumPyの力を借りれば、関数を1個呼び出すだけで計算結果が返ってきます。
numpy.linalg.svd — NumPy v1.13 Manual
変化度の定義
時刻における過去と現在の行列の主部分空間から、両者の食い違いを定量化することで、時刻
の変化度を定義します。これは部分空間同士の距離を求める問題になるのですが、空間同士の距離を求めるというのはすごく理論的にも難しいです。今回は、行列2ノルムが行列の最大特異値に等しいということを利用して、異常度を算出しています。
まとめると、特異スペクトル変換法は、まずスライド窓を使って、現在と過去のデータを行列として取り出し、特異値分解に基づき行列の特徴パターンを求め、そこから両者の変化度を算出するアルゴリズムです。
Pythonによる実装
作りたいものは、入力として1次元データのリスト（Pythonの組み込みリストやNumPyのarray）を渡すと、出力として特異スペクトル変換法に基づき変化度に変換した1次元データリストを返すというものです。
※今回処理の流れを掴みやすくするため、あえて2つの関数に詰め込んで実装しています。
import numpy as np def extract_matrix(data, start, end, w): row = w column = end - start + 1 matrix = np.empty((row, column)) i = 0 for t in range(start, end+1): matrix[:, i] = data[t-1:t-1+row] i += 1 return matrix def sst(data, w, m=2, k=None, L=None): """ Parameters ---------- data : array_like Input array or object that can be converted to an array. w : int Window size m : int Number of basis vectors k : int Number of columns for the trajectory and test matrices L : int Lag time Returns ------- Numpy array contains the degree of change. """ # Set variables if not isinstance(data, np.ndarray): data = np.array(data) if k is None: k = w // 2 if L is None: L = k // 2 T = len(data) # Calculation range start_cal = k + w end_cal = T - L + 1 # Calculate the degree of change change_scores = np.zeros(len(data)) for t in range(start_cal, end_cal + 1): # Trajectory matrix start_tra = t - w - k + 1 end_tra = t - w tra_matrix = extract_matrix(data, start_tra, end_tra, w) # Test matrix start_test = start_tra + L end_test = end_tra + L test_matrix = extract_matrix(data, start_test, end_test, w) # Singular value decomposition(SVD) U_tra, _, _ = np.linalg.svd(tra_matrix, full_matrices=False) U_test, _, _ = np.linalg.svd(test_matrix, full_matrices=False) U_tra_m = U_tra[:, :m] U_test_m = U_test[:, :m] s = np.linalg.svd(np.dot(U_tra_m.T, U_test_m), full_matrices=False, compute_uv=False) change_scores[t] = 1 - s[0] return change_scores
実際にこのコードを使って、先に挙げた周波数異常のデータを変化度に変換し、描画した結果が下記になります。（縦軸は見やすいようにスケールを調整しています。）
change_scores = sst(data, 50)
変化度を示している赤のラインを見ると周波数変化が起こっている2点で、きわめて明瞭にピークをもつことがわかります。このように、波の周波数が激しくなっている点も、ゆるやかに戻っている点のどちらも検知できています。若干遅れて検知しているのは、パラメータの調整で改善できる部分でもありますが、ある程度限界もあります。
今回の例は周波数の異常パターンでしたが、同アルゴリズムを使って、データの傾きの変化点なども検知することができます。
特異スペクトル変換法の課題
特異スペクトル変換法の課題のひとつとして、計算コストが高いことが挙げられます。一連の処理の中でも特に特異値分解の部分の計算量がネックになります。計算が重く時間を要すると言うことは、データのリアルタイム処理には使えないということであり、時系列解析の幅を制限する要因にもなり得ます。
この課題を解決する方法の一つとして、クリロフ部分空間法という手法と併用することで計算の高速化することができるようです。また、そもそもリアルタイム処理には別のアルゴリズムを使うという選択肢もあるでしょう。例えば、ChangeFinderは逐次的な計算アルゴリズムになっているため、リアルタイム処理には向いています。（一方、ハイパーパラメータの闇は深そうなので、要は適材適所を見極める必要があるということですね。）
おわりに
現在、今回記事に書いた特異スペクトル変換法による異常検知も実装した『Banpei』という異常検知パッケージを絶賛開発中です。
github.com
今後、様々なアルゴリズムを実装したり、ストリーミングデータの取り扱いができるように機能拡張したり、などいろいろ計画中です。
私自身エンジニア歴がまだ1年にも満たない未熟者なので、書き方がおかしいところなどあればぜひご指導をよろしくお願いいたしますorz
