[TL;DR] 以下の論文紹介とpython実装例: Lu, Yue, et al. "Matrix Profile XXX: MADRID: A Hyper-Anytime and Parameter-Free Algorithm to Find Time Series Anomalies of all Lengths." ICDM2023.
§1. はじめに
「ついに出た」と言いつつ論文の発表日は9ヵ月以上前ですが。
前回、前々回と時系列ディスコード(Time Series Discords)探索による単変量異常検知を紹介してきました。これらの手法は深層学習を初めとする複雑な機械学習異常検知手法と比べてパラメータが少ない点で扱いやすく、特にDAMPに至っては部分時系列長パラメータのみで異常検知が可能でした。DAMPに限らず大多数の時系列異常検知手法は部分時系列長パラメータを有していますが、このパラメータはきちんと調整しないと十分な性能が発揮できない極めて重要なパラメータです。
今回紹介するMADRID [1]では、その部分時系列長パラメータの調整さえも自動化し、パラメータ調整の手間なしに異常を求めるものです。
まずは実装例＆出力例をご紹介します。以下のような単変量時系列データがあります。
開始から10000点目（緑の部分）を学習データとして定義したとき、10000点目以後の時系列データの中から異常個所を検出するというような問題設定です。本データの異常は以下の通りです。
11500点目～12500点目を拡大しています。赤の領域部分が異常な波形となっていることが分かります。この赤の領域ですが、筆者がマニュアルで区切った訳ではなく、MADRIDによって自動で検出された異常です。
以下、は、上図を出力するためのPythonプログラムです。MADRIDのコア関数部分のプログラムは私のgithubの方に公開しましたのでそちらを参照ください。
import numpy as np
import matplotlib.pyplot as plt
from madrid import MADRID
#1次元時系列データの読み込み
TS = np.loadtxt("Data/UCR_Anomaly_FullData/119_UCR_Anomaly_ECG1_10000_11800_12100.txt")
train_test_split = 10000 #学習データとテストデータの分割点
best_discord_loc, _, _, _, _ = MADRID(TS, train_test_split)
#可視化
plt.figure(figsize=(10,2))
plt.grid()
plt.plot(TS,lw=2., color="k")
plt.axvspan(best_discord_loc[0], best_discord_loc[1],
facecolor='r', alpha=0.1,label="anomaly")
plt.xlim(best_discord_loc[0]-300,best_discord_loc[1]+300)
plt.show()
MADRID関数に入力されるのは時系列データTSと、その学習データの分割点のみで、ハイパーパラメータの調整なしに異常な部分時系列の(開始点,終了点)がタプル形式で返されます。上記プログラムは、受け取ったタプル形式の出力情報用いてmatplotlibでハイライトしてます。
§2. MADRID
MADRIDは、あらゆる長さの時系列ディスコード(Time Series Discords)を高速に探索する方法です。時系列ディスコードは前々回に説明しましたが、要は異常な部分時系列です。機械学習的な文脈で言うと、部分時系列に関する1-Nearest-Neighbor法を用いた異常検知を行ってます。
テストデータのあらゆる長さの各部分時系列に関して、(上記の例では)長さ10000の学習データで1-Nearest-Neighborなんてやったら日が暮れてしまうのではと思われるかもしれません。何の工夫もなしに実装した場合、その直感は正解でしょう。しかし、時系列ディスコードが最初に提案 [2] されて以後、高速化への努力が積みに積み重なりまして、MADIRIDにおいては、上記の「UCR_Anomaly_ECG1」データでの出力に1分40秒程度しかかかっていません。（参考までに、CPUはi5-13400です。）
§2.1 アルゴリズム概要
長年研究されてきたの時系列波形分析技術の集大成であるMADRIDのアルゴリズムを詳しく説明するのは困難ですが、ここでは前回紹介したDAMPアルゴリズムからの差分を用いて説明してみます。
DAMPは、単変量時系列データにおける時系列ディスコード検索の最先端手法ですが、パラメータに部分時系列長を必要とします。MADIRIDでは様々な部分時系列長パラメータで何度もDAMPを実行するというのが基本のメカニズムです。
（ここで§1での誇張を謝罪しなければなりません。パラメータフリーと銘打ちしましたが、実際には調べる部分時系列長パラメータの最小値と最大値を設定する必要があります。私の実装では学習データのサイズに応じてそこを自動決定していますが、これは完全なやり方ではないでしょう。）
§2.2 高速化への工夫
色々な部分時系列長でDAMPを実行するに際し、計算高速化のための重要なアイデアがあります。部分時系列長が16のとき、DAMPによってディスコードを調べて、その位置が12123だったとします。その前後の部分時系列長、例えば部分時系列長が17や15のときも、ディスコードの位置は12123かその付近ではないかと予想出来ます。そこがディスコードでなかったとしても、比較的高いディスコードスコアが得られるのではないかと考えられます。
前回の記事で述べた通り、DAMPはアルゴリズム内のループの中で、良いBest_So_Farが得られている状況であれば枝刈りを効率的に行うことが出来ます。よって、部分時系列長=16でのディスコードについて調べ終わった後、部分時系列長が17や15のときは、最初に12123の位置の部分時系列からディスコードスコアを調べ始めればよいのです。
イメージとしては下図になります。下図では、より一般化して、部分時系列長=16の検索でディスコード位置=12123を発見した後（赤■）は、部分時系列長=8~24の全てに対し、その位置のディスコードスコアを最初に計算する（黄■）ということをやっています。同様に、部分時系列長=8又は24を検索した後も同じ手続きを行います。
出展： [1] の図9
§3. その他の出力例
昆虫のEPGデータ
UCR Anomaly Dataset [3]より、Asian Citrus Psyllidという虫にEPG装置を取り付けて得られた時系列データです。
TS = np.loadtxt("UCR_Anomaly_FullData/147_UCR_Anomaly_Lab2Cmac011215EPG3_5000_16390_16420.txt")
train_test_split = 5000
best_discord_loc,_,_,_,_ = MADRID(TS, train_test_split)
plt.figure(figsize=(10,2))
plt.grid()
plt.plot(TS,lw=2., color="k")
plt.axvspan(best_discord_loc[0], best_discord_loc[1], facecolor='r', alpha=0.1,label="anomaly")
plt.xlim(best_discord_loc[0]-300,best_discord_loc[1]+300)
イタリアの電力需要データ
同じくUCR Anomaly Archive [3] より、Asian Citrus Psyllidという虫にEPG装置を取り付けて得られた時系列データです。異常個所では、3波形が連続して異なるパターンとなってます。人のチェックでも結構注意して見ないと気付かないんじゃないかなと思いますが、MADRIDでは自動で検出してくれました。
TS = np.loadtxt("Data/UCR_Anomaly_FullData/155_UCR_Anomaly_PowerDemand4_18000_24005_24077.txt")
train_test_split = 10000
best_discord_loc,_,_,_,_ = MADRID(TS, train_test_split)
plt.figure(figsize=(10,2))
plt.grid()
plt.plot(TS,lw=2., color="k")
plt.axvspan(best_discord_loc[0], best_discord_loc[1], facecolor='r', alpha=0.1,label="anomaly")
plt.xlim(best_discord_loc[0]-150,best_discord_loc[1]+150)
産業用コンベヤシステムのセンサデータ
最後に、本家論文でも実験に使われていた産業用コンベヤシステムのセンサデータ（通称HRSS）です。非常に些細な波形異常を上手く拾えてます。
以上、様々なデータでの異常検知結果の紹介でした。公平性のため述べておかないといけないこととして、上記のデータでは検出が上手くいったものの、どうみても異常波形っぽいのに上手く拾ってくれないデータもちらほらありました。こちらは技術的な深堀の余地がありそうです。
§4. おわりに
本記事では、パラメータフリー単変量時系列異常検知技術MADRIDを紹介しました。自身の専門性バイアスがかかってるのかもしれませんが、現状単変量時系列異常検知を行うのであれば真っ先に検討すべき技術なんじゃないかなと思ってます。
MADRIDのPython版プログラムの方は私のGithubの方で公開しました。私の知る限り、MADRIDのPython実装はこれが初です。
とは言え、機能面では最低限しか備えておらず、著者のオリジナルMATLABプログラムはいくつかの改良もされてるみたいなので、MATLABが使える環境にある方はそちらを参照する方がいいかもしれません。
参考文献
[1] Lu, Yue, et al. "Matrix Profile XXX: MADRID: A Hyper-Anytime and Parameter-Free Algorithm to Find Time Series Anomalies of all Lengths." 2023 IEEE International Conference on Data Mining (ICDM). IEEE, 2023.
[2] Keogh, Eamonn, Jessica Lin, and Ada Fu. "Hot sax: Efficiently finding the most unusual time series subsequence." Fifth IEEE International Conference on Data Mining (ICDM). IEEE, 2005.
[3] Wu, Renjie, and Eamonn J. Keogh. "Current time series anomaly detection benchmarks are flawed and are creating the illusion of progress." IEEE transactions on knowledge and data engineering 35.3 (2021): 2421-2429.