OmniAnomaly：確率的RNNによる多変量時系列の異常検知を徹底解説

サーバー群のCPU・メモリ・ネットワーク、各種機器の計測値、ロボットのセンサー。こうした「装置(エンティティ)」は数十本の時系列を同時に吐き出します。やっかいなのは、これらの系列が強い時間依存を持つ一方で、本質的に確率的なゆらぎを含むことです。たとえばネットワークの再送回数は、正常でも環境次第で大きく揺れます。この「揺れ」を異常と取り違えれば誤報だらけになり、揺れに紛れた本物の異常は見逃されます。

前回解説したUSADは、敵対的オートエンコーダで「軽量・高速・安定」を実現した手法でした。実はそのUSADが比較対象として何度も持ち出していた当時の最強ベースラインが、本記事の主役 OmniAnomaly です。Su らがKDD 2019で発表しました(Su et al., 2019)。

OmniAnomalyの核心は、GRU(RNN)とVAEを融合した「確率的RNN」で、時間依存と確率的ゆらぎの両方を同時にモデル化することです。正常パターンを頑健な確率的潜在表現として学び、再構成確率(reconstruction probability) が低い観測を異常と判定します。さらに、(1) 確率変数どうしを時間的につなぐ「確率変数接続」、(2) 事後分布を非ガウスにする「planar正規化フロー」、(3) 極値理論にもとづくしきい値の自動決定(POT)、(4) どのセンサーが異常の原因かを示す異常解釈、という実運用に効く4点セットを備えています。

OmniAnomalyを理解すると、次の場面で武器になります。

• IT運用の自動化(AIOps): サーバー群の多変量メトリクスから障害を検出し、原因メトリクスまで特定
• 機器のヘルスモニタリング: 予測困難な計測値でも、再構成ベースで頑健に異常検知
• 産業設備・ロボットの監視: ノイズの大きいセンサー群から微小な異常を拾う

本記事の内容

• なぜ「確率的RNN」なのか — 時間依存 × 確率性の両立
• GRU+VAEのアーキテクチャ(qnet/pnet)と確率変数接続・planar NF
• ELBOによる学習と、再構成確率による異常スコアの定義
• POT(Peaks-Over-Threshold) による自動閾値選択の数理
• 各次元の再構成確率を使った異常解釈
• サーベイ: SMAP/MSL/SMD の仕様と、F1比較・POT有効性・アブレーション結果
• PyTorchによる簡易OmniAnomalyの実装と合成データでの再現実験

前提知識

この記事を読む前に、以下を読んでおくと理解が深まります。

USAD：敵対的オートエンコーダによる多変量時系列の異常検知

OmniAnomalyを当時の最強ベースラインとして比較した後続手法。対比すると理解が深まります。

変分オートエンコーダ(VAE)の理論

ELBO・再構成・KLダイバージェンス。OmniAnomalyの土台です。

OmniAnomalyとは — 「正常パターン」を確率的に覚える

まず全体像をつかみましょう。OmniAnomalyの発想は「異常とは、正常パターンから大きく外れた観測である」という素朴な直感に立ちます。そこで、正常データだけを使って「正常とはどういう状態か」を低次元の潜在表現として学び、その表現でうまく再構成できない観測=異常とみなします。

ここまではオートエンコーダ系と同じですが、OmniAnomalyは2つの工夫を加えます。

第一に、GRUで時間方向のパターンを捉えます。サーバーや機器の挙動は「いつもこの順序で立ち上がって落ち着く」といった文脈を持つため、各時刻を独立に見るのではなく、過去の流れを踏まえて判断します。

第二に、潜在表現を確率分布として持ちます(VAE)。決定的に「1点」で表すのではなく「だいたいこのあたり、ばらつきはこのくらい」と分布で持つことで、時系列の確率的なゆらぎまで表現できます。

そして判定には再構成確率を使います。学習した正常パターンのもとで、ある観測が「どれだけありそうか」を確率で測り、確率が低いほど異常とみなすのです。

論文の全体パイプラインを俯瞰しておきましょう。

出典: Su et al., OmniAnomaly, KDD 2019, Fig.2

この図は2系統に分かれます。上段のオフライン訓練は「データ前処理 → モデル訓練 → しきい値選択」で、正常データだけから確率的RNNを学び、POTでしきい値を一度だけ決めます。下段のオンライン検知は「オンライン検知 → 異常結果 → 異常解釈」で、新しい観測に異常スコアを付け、しきい値で判定し、最後にどの次元が原因かまで提示します。重要なのは、訓練で得た「モデル」と「しきい値」が破線/実線でオンライン側に流れ込む点で、運用時には推論だけが走る——つまり訓練と検知が明確に分離された設計です。本記事は、この各ブロックの中身を上流から順にほどいていきます。

なぜわざわざ確率分布で持つのか。その理由を、決定的なモデルと比べて見ていきましょう。

なぜ確率的なのか — ゆらぎを表現できる強み

左図は決定的なモデル(予測ベースのLSTM-NDTなど)のイメージです。各観測を潜在空間の「1点」に押し込めます。これは「未来は1つに決まる」という前提で、予測誤差で異常を測ります。しかし現実の時系列には本質的に予測できないゆらぎがあり、正常なゆらぎを異常と誤判定しやすくなります。

右図がOmniAnomalyです。潜在を「分布」として持つため、「正常な範囲のばらつき」を分布の広がりとして吸収できます。ノイズの大きい区間では分布を広く、安定した区間では狭く——というように、ゆらぎの大きさそのものを学習できるわけです。結果として、正常なゆらぎに振り回されず、真に分布から外れた異常だけを拾えます。

これがOmniAnomalyの「Robust(頑健)」の正体です。論文も「再構成ベースのモデルは予測ベースより頑健。なぜなら実データはしばしば予測不能だから」と述べています。

この頑健さは潜在空間の見え方にも現れます。論文は学習後の3次元潜在 $z$ を可視化しています。

出典: Su et al., OmniAnomaly, KDD 2019, Fig.6

左が全体、右がその拡大です。点は時刻ごとの潜在表現で、正常も異常もひとかたまりの滑らかな帯を作っています。注目は、異常時刻の点(色の違う点)が潜在空間で大きく飛ばず、正常域の近くに留まっている点です。これは「異常を $z$ 空間で無理に遠ざける」のではなく、「正常パターンを滑らかな低次元多様体として学び、再構成確率の低さで異常をあぶり出す」という設計を裏づけます。だからこそ正常のゆらぎに $z$ が過敏に反応せず、頑健になるのです。

では、この確率的RNNを具体的にどう組むのか。ネットワーク構造を見ましょう。

アーキテクチャ — qnetとpnet

OmniAnomalyは、VAEと同じく2つのネットワークからなります。観測 $x$ から潜在 $z$ を推論する qnet($q_\phi(z|x)$)と、潜在 $z$ から観測を生成(再構成)する pnet($p_\theta(x|z)$)です。

論文のグラフィカルモデルとネットワーク詳細図を見ると、設計思想が一目でわかります。

出典: Su et al., OmniAnomaly, KDD 2019, Fig.3

上段(a)が確率変数の依存関係を表すグラフィカルモデルです。(a1)qnet では観測 $x_t$ と決定的なメモリ $e_t$ から潜在 $z_t$ を推論し、(a2)pnet では潜在 $z_t$ から再構成 $x’_t$ を生成します。両者に共通する肝は、潜在ノード $z_{t-1} \to z_t$ の横向きの矢印——つまり潜在変数そのものが時間方向に連鎖している点です(白丸が確率変数 $z, x$、黒丸が決定的変数 $e, d$)。下段(b)は実装で、左の(b1)qnetは「GRU → Concat($z_{t-1}$と$e_t$) → Dense層 → Linear/Softplusで$\mu,\sigma$ → Planar Normalizing Flow」、右の(b2)pnetは対称に「GRU → Dense層 → Linear/Softplusで再構成の$\mu,\sigma$」という流れです。確率変数接続(Concat)と正規化フローが、ありふれたVAEとの差分だと視覚的に確認できます。

この図のデータフローを、式で1本ずつ追っていきましょう。

qnet(推論ネットワーク)

時刻 $t$ で、観測 $x_t$ と前ステップのGRU隠れ状態 $e_{t-1}$ をGRUセルに入れて $e_t$ を得ます。

$$ \begin{equation} e_t = (1 – c^e_t)\circ\tanh(w_e x_t + u_e(r^e_t\circ e_{t-1}) + b_e) + c^e_t\circ e_{t-1} \end{equation} $$

ここで $r^e_t$ はリセットゲート、$c^e_t$ は更新ゲートで、$\circ$ は要素積です。この決定的な $e_t$ が「長期の時間情報を運ぶ内部メモリ」の役割を果たします。

次が重要な工夫、確率変数接続(stochastic variable connection) です。$e_t$ に前ステップの潜在変数 $z_{t-1}$ を連結して全結合層 $h_\phi$ に通し、潜在 $z_t$ の平均と標準偏差を出します。

$$ \begin{equation} \mu_{z_t} = w_{\mu_z}\,h_\phi([z_{t-1}, e_t]) + b_{\mu_z} \end{equation} $$

$$ \begin{equation} \sigma_{z_t} = \mathrm{softplus}\big(w_{\sigma_z}\,h_\phi([z_{t-1}, e_t]) + b_{\sigma_z}\big) + \epsilon \end{equation} $$

$z_{t-1}$ を入力に使うことで、潜在変数どうしが時間的に依存します。これがLSTM-VAE(単にVAEのMLPをLSTMに置き換えただけ)との決定的な違いです。LSTM-VAEは観測の時間依存は扱えても、潜在変数の時間依存を無視していました。OmniAnomalyは潜在の系列にも依存構造を入れることで、より良い表現を学びます。

planar正規化フロー

$\mu_{z_t}, \sigma_{z_t}$ から得た $z^0_t \sim \mathcal{N}(\mu_{z_t}, \sigma_{z_t}^2 I)$ は対角ガウスです。しかし真の事後分布はガウスとは限りません。そこでplanar正規化フローで、可逆変換を $K$ 回繰り返して非ガウスに変形します。

$$ \begin{equation} z^k_t = f^k(z^{k-1}_t) = z^{k-1}_t + u\,\tanh(w^\top z^{k-1}_t + b) \end{equation} $$

最終出力 $z_t = z^K_t$ を潜在変数として使います。

なぜこの形なら確率密度を正しく追えるのか。変数変換の公式より、変換後の密度は元の密度にヤコビアンの行列式を掛けて補正します。対数で書くと、フロー全体の対数密度は次のように初期ガウスの対数密度から各層の対数ヤコビアン行列式を引いた形になります。

$$ \begin{equation} \log q_K(z^K_t) = \log q_0(z^0_t) – \sum_{k=1}^{K}\log\left|\det\frac{\partial f^k}{\partial z^{k-1}_t}\right| \end{equation} $$

planar変換のありがたさは、この行列式が $D\times D$ 行列を組まずにスカラー1個で計算できる点にあります。$\psi(z) = (1-\tanh^2(w^\top z + b))\,w$ とおくと、

$$ \begin{equation} \left|\det\frac{\partial f^k}{\partial z}\right| = \left|1 + u^\top \psi(z)\right| \end{equation} $$

と、行列式が内積1回で済みます($\partial f/\partial z = I + u\,\psi(z)^\top$ がランク1更新なので行列式定理が使える)。これが潜在次元 $D$ に対して $O(D)$ で回り、フロー長 $K=20$ を重ねても軽い理由です。なお変換が可逆であるには $w^\top u \geq -1$ が必要で、実装ではパラメータ $u$ をこの条件を満たすよう再パラメータ化します(後述コードの u_hat)。

上図のように、単純なガウス(左)に可逆変換を重ねると、曲がった複雑な分布(右)へと変形できます。各変換は可逆なので確率密度を正しく追跡でき、ガウスでは表せない事後分布を柔軟に近似できます。これが入力データの複雑なパターンを捉える助けになります。式の上では、上の対数ヤコビアン項がそのままELBOの $\log q_\phi(z|x)$ に効いてきます。後段の学習でこの項がどう現れるかを見ましょう。

pnet(生成ネットワーク)

pnetはqnetと対称的な構造で、潜在 $z_t$ から観測を再構成します。まず線形ガウス状態空間モデル(SSM)で $z$ 空間の変数を時間的につなぎ($z_t = O_\theta(T_\theta z_{t-1} + v_t) + \epsilon_t$)、GRUセル $d_t$ を経て、全結合層 $h_\theta$ から再構成の平均 $\mu_{x_t}$ と標準偏差 $\sigma_{x_t}$ を出します。再構成 $x’_t$ は $\mathcal{N}(\mu_{x_t}, \sigma_{x_t}^2 I)$ からサンプリングされます。

ここがポイントです。再構成を「1点 $\mu_x$」ではなく「分布 $\mathcal{N}(\mu_x, \sigma_x^2)$」として出すので、後で確率として異常を測れます。学習で $\sigma_x$ も推定するため、ゆらぎの大きい次元は自動的に「許容幅」を広げられます。

ネットワークができたら、これを正常データだけで学習します。

学習：ELBOの最大化

OmniAnomalyはVAEと同じく、ELBO(変分下限) を最大化して学習します。長さ $T+1$ の系列 $x_{t-T:t}$ に対し、損失は次の形です。

$$ \begin{equation} \mathcal{L}(x_{t-T:t}) \approx \frac{1}{L}\sum_{l=1}^{L}\Big[\log p_\theta(x_{t-T:t}\,|\,z^{(l)}_{t-T:t}) + \log p_\theta(z^{(l)}_{t-T:t}) – \log q_\phi(z^{(l)}_{t-T:t}\,|\,x_{t-T:t})\Big] \end{equation} $$

3つの項の意味を順に見ましょう。第1項 $\log p_\theta(x|z)$ は負の再構成誤差(再構成確率の対数)で、「正常をどれだけ忠実に再現できたか」。第2項と第3項の和 $\log p_\theta(z) – \log q_\phi(z|x)$ はKL正則化で、潜在分布を事前分布に近づける役割です。第2項の事前 $\log p_\theta(z)$ は線形ガウスSSMで時間的につながれ、第3項の事後 $q_\phi(z)$ はplanar NFで変形されています。

ここで先ほどの正規化フローが効いてきます。第3項の事後密度はフローを通した後の密度なので、

$$ \begin{equation} \log q_\phi(z_t|x_{t-T:t}) = \log q_0(z^0_t) – \sum_{k=1}^{K}\log\left|1 + u_k^\top \psi_k(z^{k-1}_t)\right| \end{equation} $$

と、初期ガウスの対数密度から対数ヤコビアン項を引いた形で計算されます。$\log|\det|$ が引かれることで「フローが分布を広げた/曲げたぶん」が密度に正しく反映され、$z$ がガウスから離れるほど第3項は小さくなります。これは「事後をガウスに縛らず、必要なら非ガウスに動いてよい」という自由を与えつつ、その自由のコストを密度で会計する仕組みです。SSM事前(第2項)とフロー事後(第3項)が噛み合って、潜在系列に時間的な滑らかさを保たせます。

平たく言えば「正常をよく再構成しつつ、潜在空間を素直な形に保つ」よう学習します。再構成確率の期待値はモンテカルロで $L$ 個の $z$ をサンプリングして近似し(再パラメータ化トリック)、勾配を流します。最適化はAdam、勾配爆発を防ぐため勾配クリッピング(ノルム10)を使います。

それでは、学習済みモデルでどう異常を判定するのでしょうか。

推論：再構成確率による異常スコア

検出時、観測 $x_t$ の異常スコアは再構成確率そのものです。直前の $T$ 観測 $x_{t-T:t}$ から $x_t$ を再構成し、その条件付き確率の対数を取ります。

$$ \begin{equation} S_t = \log p_\theta(x_t\,|\,z_{t-T:t}) \end{equation} $$

スコアの読み方に注意です。$S_t$ が高いほど「正常パターンに沿っていて、ありそう」=正常。低いほど「ありそうにない」=異常です。USADの再構成誤差(高いほど異常)とは符号が逆になります。上図のように、正常な観測は再構成分布のピーク付近にあり高確率、異常な観測は裾にあり低確率です。実装では扱いやすさのため、スコアを反転した負の対数尤度(NLL) を使い「高いほど異常」に揃えることが多いです。

ここがUSADのような決定的モデルとの本質的な差です。ガウス出力のNLLを書き下すと、

$$ \begin{equation} -\log p_\theta(x^i_t|z) = \frac{(x^i_t – \mu^i_{x_t})^2}{2\,(\sigma^i_{x_t})^2} + \log \sigma^i_{x_t} + \frac{1}{2}\log 2\pi \end{equation} $$

となります。USADなどの再構成誤差(MSE)が $(x – \mu_x)^2$ そのものなのに対し、OmniAnomalyは誤差を学習した分散 $(\sigma^i_{x_t})^2$ で割って正規化します。だから「正常でも大きく揺れる次元」は $\sigma_x$ を大きく学んで誤差を割り引き、揺れの小さい次元はわずかなずれも見逃さない——という次元ごと・時刻ごとの自動スケーリングが効きます。決定的モデルはこの分母を持たないので、ノイズの大きい正常区間を異常と取り違えやすいのです。後の実装でこの差をはっきり数値で確認します。

スコアが出たら、最後に「どこからを異常とするか」のしきい値が必要です。OmniAnomalyはこれも自動化します。

POT：極値理論による自動しきい値

異常検知の地味だが重要な難所がしきい値の決定です。手で調整すると運用が回りません。OmniAnomalyは極値理論(EVT) の POT(Peaks-Over-Threshold) でこれを自動化します。

POTの発想はこうです。データ全体の分布を仮定するのは難しいが、「裾(極端な値)」だけなら一般化パレート分布(GPD)でよく近似できる、という極値理論の定理に基づきます。これはPickands–Balkema–de Haanの定理で、十分高いしきい値 $th$ を超える超過量 $X – th$ の条件付き分布が、しきい値を上げるにつれてGPDに収束する、という主張です。

$$ \begin{equation} \bar{F}_{th}(s) = P(X – th > s \mid X > th) \simeq \left(1 + \frac{\hat{\gamma}\,s}{\hat{\beta}}\right)^{-1/\hat{\gamma}} \end{equation} $$

つまり「異常スコアの大きい側がどう減衰するか」だけをGPDで素直にモデル化できるわけです。手順は、まず初期しきい値(高い分位点)$th$ を置き、それを超える超過量だけを取り出してGPDを最尤推定で当てはめ、目標の超過確率 $q$ に対応する最終しきい値を外挿します。

$$ \begin{equation} th_F \simeq th – \frac{\hat{\beta}}{\hat{\gamma}}\left(\left(\frac{q N’}{N’_{th}}\right)^{-\hat{\gamma}} – 1\right) \end{equation} $$

ここで $\hat{\gamma}, \hat{\beta}$ はGPDの形状・尺度パラメータ(最尤推定)、$N’$ は観測数、$N’_{th}$ は初期しきい値を超えた数です。調整するのは「初期分位点」と「超過確率 $q$」の2つだけで、データ分布に仮定を置かない点が強みです。後の実験で、POTが全探索で求めた最良しきい値にほぼ匹敵することを確認します。

異常を検出したら、運用者は次に「なぜ異常なのか」を知りたくなります。

異常解釈：原因センサーの特定

OmniAnomalyのもう一つの貢献が異常解釈です。再構成確率は次元ごとに分解できます。

$$ \begin{equation} S_t = \sum_{i=1}^{M} S^i_t, \qquad S^i_t = \log p_\theta(x^i_t\,|\,z_{t-T:t}) \end{equation} $$

各次元(センサー) $i$ の再構成確率 $S^i_t$ を計算し、確率が低い(=再構成しにくい)次元ほど異常への寄与が大きいとみなします。これらを昇順に並べた上位数本を「原因候補」として運用者に提示します。

上図は、後述する合成データで検出した異常窓について、各次元の異常寄与(平均NLL)を示したものです。寄与が突出した次元(赤)が「原因センサー」候補になります。実運用では、ネットワーク障害ならネットワーク系メトリクスが上位に来る、といった形でトラブルシュートを助けます。論文ではSMDで解釈精度 HitRate@100%=0.80, @150%=0.89 を達成し、LSTM-VAE(0.50/0.62)を大きく上回りました。

理論が揃ったので、論文の実験をサーベイとしてまとめます。

サーベイ①：データセット

OmniAnomalyは3つの多変量時系列データセットで評価されました。うちSMDはこの論文が新たに公開したものです。

いずれも複数の機体・サーバーの集合体で、各エンティティを個別に学習・評価して平均します。観測はすべて1分間隔です。USAD・TranADなど後続研究もこの3つを標準ベンチマークとして使い続けており、多変量時系列異常検知の「定番セット」になっています。

主要なハイパーパラメータも公開されています(Appendix B)。

潜在次元がわずか3で十分なのが目を引きます。論文の感度分析(Fig.7)では、z次元を3〜32で変えてもF1はほぼ変わらず、「小さすぎるとunderfit、ある程度以上は安定」という挙動でした。

出典: Su et al., OmniAnomaly, KDD 2019, Fig.7

3データセット(SMAP/MSL/SMD)とも $z=3$ あたりで立ち上がり、それ以降は横ばいです。とくにSMAP(青)は $z=1,2$ で急落していて、潜在が小さすぎると正常パターンを表しきれずunderfitすることがわかります。一方 $z\geq 3$ では32次元まで上げてもほとんど変わらないため、論文は計算コストの軽い $z=3$ を全データで採用しています。この「ある程度以上は鈍感」という性質は、後ほど合成データでも自分で再現します。

サーベイ②：検出性能(F1)の比較

LSTM-NDT・DAGMM・LSTM-VAEなど代表的な教師なし手法との比較が次表です(point-adjust適用のF1)。太字が各列の最良値です。

point-adjustを使わない「素のF1best」(全しきい値を探索した上限)でも傾向は同じで、論文Fig.4が手法横断で示しています。

出典: Su et al., OmniAnomaly, KDD 2019, Fig.4

紫(OmniAnomaly)は MSL・SMD・Total で頭一つ抜けており、SMAP のみ赤(LSTM-NDT)にわずかに譲ります。注目すべきは負け方のばらつきが小さいことで、Donut・EncDec-AD・DAGMM・LSTM-VAE はデータセットごとに大きく上下するのに対し、OmniAnomalyはどの棒も高い位置で安定しています。これが「3データセットすべてで適合率・再現率がともに高い」という頑健性の見える化です。

OmniAnomalyはMSL・SMDで最良、SMAPでもLSTM-NDTにわずかに及ばない2位です。全体ではF1 0.86を達成し、当時の最強ベースライン(LSTM-NDT)を0.09上回りました。特筆すべきは頑健性で、OmniAnomalyは3データセットすべてで適合率・再現率がともに0.74超(SMAP P0.74/R0.98、MSL P0.89/R0.91、SMD P0.83/R0.94)。どのベースラインもこれを達成できませんでした。

各手法の「負け方」も示唆に富みます。DAGMMは時間情報を使わない(1観測のみ入力)ため、時系列では不利。LSTM-NDTは予測ベースで、MSLのような予測困難な系列で崩れます。LSTM-VAEは観測の時間依存は扱えても潜在変数の時間依存を無視するため、OmniAnomalyの確率変数接続に一歩譲ります。「時間依存」と「確率性」の両立がいかに効くか、という論文の主張を裏づける結果です。

サーベイ③：POTの有効性とアブレーション

POTで自動選択したしきい値が、どれだけ「理論上の最良」に近いかを見ます。F1(POT)と、全しきい値を探索した F1best の比較です。

差はわずか0.003〜0.077で、POTは手探りなしに最良に近いしきい値を選べることが分かります(右は後述の合成データでの実測)。実運用では「全探索」は正解ラベルが要るので不可能ですから、ラベル不要で良い閾値を出せるPOTの価値は大きいです。

さらに論文は4つの構成要素(GRU/確率変数接続/planar NF/POT)のアブレーションを行い、いずれも性能に寄与することを確認しました。GRUは単純RNNやLSTMより良く、確率変数接続はqnet・pnet両方に入れるのが最良、planar NFは事後分布の表現力を上げて性能を改善します(MSLでは効果が小さい——疎な入力は単純な潜在でも表せるため)。

ちなみに前回のUSADは、まさにこのOmniAnomalyを主な比較対象に据え、「同等の精度を平均547倍速く」出すことを売りにしました。OmniAnomalyは精度面で強力ですが、確率的RNN+正規化フローのぶん計算は重い、というトレードオフです。

数式とサーベイを踏まえて、実際に手を動かしましょう。

Pythonでのスクラッチ実装

ここからは簡易版OmniAnomaly(GRU-VAE + planar NF)をPyTorchで実装し、合成多変量時系列で再現します。確かめたいのは、(1) 確率的モデルがノイズの大きい正常に惑わされず異常を捉えるか、(2) POTが妥当な閾値を自動で出すか、(3) 潜在次元にどれだけ頑健か、です。

不等分散な合成データ

OmniAnomalyの強みが出るのは「正常でもゆらぎが大きい」データです。そこで一部チャネルに時間変動するノイズ振幅(不等分散) を持たせ、テスト系列には正常に近い異常(振幅変化・相関破壊・固着・点異常)を埋め込みます。

import numpy as np
import torch
import torch.nn as nn
from scipy.stats import genpareto
from sklearn.metrics import roc_auc_score, average_precision_score, precision_recall_curve, f1_score

np.random.seed(3); torch.manual_seed(3)
M, T, EPOCHS = 5, 30, 45   # 変量数・系列長(窓長)・エポック

def gen_series(N, anomalies=False, seed=0):
    rng = np.random.default_rng(seed); t = np.arange(N)
    f1 = np.sin(0.05*t); f2 = np.sin(0.017*t + 1.0)
    base = np.vstack([f1+0.3*f2, 0.8*f1-0.2*f2, 0.5*f1+0.6*f2,
                      np.sin(0.11*t), 0.4*f2+0.3*np.sin(0.07*t)])
    # 不等分散ノイズ:一部chは時間変動するノイズ振幅(正常でもゆらぐ=確率モデル向き)
    noise_amp = np.vstack([0.05+0.30*np.abs(np.sin(0.02*t)),
                           0.05+0.25*np.abs(np.cos(0.015*t)),
                           0.06*np.ones(N), 0.06*np.ones(N), 0.06*np.ones(N)])
    X = base + noise_amp*rng.standard_normal((M, N)); lab = np.zeros(N, int)
    if anomalies:
        for s0 in [300, 1100]:                       # 振幅変化
            ln = 15; X[:, s0:s0+ln] *= 1.30; lab[s0:s0+ln] = 1
        for s0, ch in [(560, 2), (1500, 3)]:         # 相関破壊
            ln = 14; X[ch, s0:s0+ln] = -X[ch, s0:s0+ln]; lab[s0:s0+ln] = 1
        s0 = 1800; ln = 16; X[:, s0:s0+ln] = X[:, s0:s0+1]; lab[s0:s0+ln] = 1  # 固着
        X[2, 820:822] += 1.1; lab[820:822] = 1       # 点異常
    return X.T, lab

Xtr, _ = gen_series(2600, False, 11); Xte, lab = gen_series(2200, True, 29)
mu, sd = Xtr.mean(0), Xtr.std(0)+1e-8; Xtr, Xte = (Xtr-mu)/sd, (Xte-mu)/sd

def seqs(X, lab=None):
    N = len(X); W = np.stack([X[i:i+T] for i in range(N-T+1)])
    if lab is None: return W.astype(np.float32)
    y = np.array([lab[i:i+T].max() for i in range(N-T+1)])
    return W.astype(np.float32), y

Wtr = seqs(Xtr); Wte, yte = seqs(Xte, lab)
Wtr_t, Wte_t = torch.tensor(Wtr), torch.tensor(Wte)
print("系列数 train=%d test=%d 異常率=%.3f" % (len(Wtr), len(Wte), yte.mean()))

チャネル1・2は時間とともにノイズ振幅が0.05〜0.35で変動します。これは「正常だがゆらぐ」状況で、純粋な再構成誤差(MSE)を使うモデルはこのゆらぎを誤検出しがちです。異常率は約0.12になります。

planar正規化フロー層

可逆性を保つため、パラメータ $u$ を再パラメータ化(u_hat)し、対数ヤコビアン行列式も返します。

class Planar(nn.Module):
    def __init__(self, d):
        super().__init__()
        self.w = nn.Parameter(torch.randn(d)*0.1)
        self.u = nn.Parameter(torch.randn(d)*0.1)
        self.b = nn.Parameter(torch.zeros(1))
    def forward(self, z):
        wu = (self.w*self.u).sum()
        u_hat = self.u + (-1 + torch.log1p(torch.exp(wu)) - wu)*self.w/(self.w.pow(2).sum()+1e-8)
        lin = z@self.w + self.b
        f = z + torch.outer(torch.tanh(lin), u_hat)
        psi = (1 - torch.tanh(lin)**2).unsqueeze(-1)*self.w
        logdet = torch.log(torch.abs(1 + psi@u_hat) + 1e-8)
        return f, logdet

u_hatへの変換は、planar変換が可逆になるための条件 $w^\top u \geq -1$ を満たすための標準的なテクニックです。

GRU-VAE本体

GRUエンコーダ→確率的潜在(+planar NF)→GRUデコーダ→ガウス出力、という構成です。prob/stochastic/use_flowフラグで、OmniAnomaly・確率なしの純MSE-AE・フローなしVAEを切り替えられます。

class OmniVAE(nn.Module):
    def __init__(self, M, H=32, Z=3, use_flow=True, stochastic=True, prob=True):
        super().__init__()
        self.stochastic, self.use_flow, self.Z, self.prob = stochastic, use_flow, Z, prob
        self.enc = nn.GRU(M, H, batch_first=True)
        self.mu = nn.Linear(H, Z); self.ls = nn.Linear(H, Z)
        self.flows = nn.ModuleList([Planar(Z) for _ in range(2)]) if use_flow else None
        self.dec = nn.GRU(Z, H, batch_first=True)
        self.xmu = nn.Linear(H, M); self.xls = nn.Linear(H, M)
    def forward(self, x, sample=None):
        if sample is None: sample = self.stochastic
        B, L, _ = x.shape
        h, _ = self.enc(x); mu = self.mu(h); ls = torch.clamp(self.ls(h), -6, 3)
        z = mu + torch.randn_like(mu)*torch.exp(0.5*ls) if sample else mu
        logdet = torch.zeros(B, L)
        if self.use_flow and self.stochastic:
            zf = z.reshape(B*L, self.Z)
            for fl in self.flows:
                zf, ld = fl(zf); logdet = logdet + ld.reshape(B, L)
            z = zf.reshape(B, L, self.Z)
        d, _ = self.dec(z); xmu = self.xmu(d); xls = torch.clamp(self.xls(d), -6, 3)
        return xmu, xls, mu, ls, logdet
    def loss(self, x):
        xmu, xls, mu, ls, logdet = self(x)
        if not self.prob:                                   # 純MSE(決定的AE)
            return ((x-xmu)**2).sum(-1).mean()
        nll = 0.5*(((x-xmu)**2)*torch.exp(-xls) + xls + np.log(2*np.pi)).sum(-1)  # 再構成NLL
        if self.stochastic:
            kl = -0.5*(1 + ls - mu**2 - torch.exp(ls)).sum(-1) - logdet           # flow補正
            return (nll + kl).mean()
        return nll.mean()
    @torch.no_grad()
    def score(self, x):                                     # 異常スコア(高い=異常)。決定的(z=μ)で再現性確保
        xmu, xls, _, _, _ = self(x, sample=False)
        if not self.prob:
            return ((x-xmu)**2).mean(-1).mean(1).cpu().numpy()
        nll = 0.5*(((x-xmu)**2)*torch.exp(-xls) + xls + np.log(2*np.pi)).sum(-1)
        return nll.mean(1).cpu().numpy()

prob=Trueでは再構成を分布として扱い、ゆらぎの大きい次元は $\sigma_x$ を大きく学んで「許容」します。prob=Falseの純MSE-AEはこの調整ができません。

学習と比較

def train(model, epochs=EPOCHS, seed=3, lr=3e-3):
    torch.manual_seed(seed)
    opt = torch.optim.Adam(model.parameters(), lr=lr); hist = []
    for ep in range(epochs):
        perm = torch.randperm(len(Wtr_t)); tot = []
        for i in range(0, len(Wtr_t), 128):
            L = model.loss(Wtr_t[perm[i:i+128]])
            opt.zero_grad(); L.backward()
            torch.nn.utils.clip_grad_norm_(model.parameters(), 10.0)  # 勾配クリッピング
            opt.step(); tot.append(float(L))
        hist.append(np.mean(tot))
    return hist

def best_f1(y, sc):
    p, r, _ = precision_recall_curve(y, sc); f = 2*p*r/(p+r+1e-12); return float(np.nanmax(f))

omni = OmniVAE(M, use_flow=True, stochastic=True, prob=True);  h_o = train(omni)
ae   = OmniVAE(M, use_flow=False, stochastic=False, prob=False); h_a = train(ae)   # 純MSE GRU-AE
sc_o = omni.score(Wte_t); sc_a = ae.score(Wte_t)
print("OmniAnomaly AUC=%.3f AP=%.3f F1=%.3f" % (
    roc_auc_score(yte, sc_o), average_precision_score(yte, sc_o), best_f1(yte, sc_o)))
print("GRU-AE(MSE) AUC=%.3f AP=%.3f F1=%.3f" % (
    roc_auc_score(yte, sc_a), average_precision_score(yte, sc_a), best_f1(yte, sc_a)))

実行結果は次のとおりです。

訓練損失は安定して収束します。確率的モデルでも、勾配クリッピングと適切な学習率で問題なく学習できます。

OmniAnomalyはGRU-AEを大きく上回りました(AUC 0.887 vs 0.711、AP 0.716 vs 0.363)。理由は明快です。純MSE-AEは「ノイズの大きい正常区間」も誤差大とみなして誤検出しますが、OmniAnomalyは $\sigma_x$ を学んでそのゆらぎを許容するため、真に分布から外れた異常だけを拾えます。確率的モデリングの効果が、不等分散データではっきり出ています。

POTで閾値を自動決定

def pot_threshold(train_scores, q=0.05, init_q=0.90):
    init = np.quantile(train_scores, init_q)
    peaks = train_scores[train_scores > init] - init
    c, loc, scale = genpareto.fit(peaks, floc=0)        # GPDを裾に当てはめ
    n, Nt = len(train_scores), len(peaks)
    if abs(c) < 1e-6:
        return init - scale*np.log(q*n/Nt)
    return init + (scale/c)*((q*n/Nt)**(-c) - 1)         # 式(7)

thF = pot_threshold(omni.score(Wtr_t), q=0.05, init_q=0.90)
pred = (sc_o >= thF).astype(int)
print("POT閾値=%.3f  F1(POT)=%.3f  F1(best)=%.3f" % (thF, f1_score(yte, pred), best_f1(yte, sc_o)))

結果は F1(POT)=0.579, F1(best)=0.660。正解ラベルを一切使わずに、全探索の最良(0.660)に近い閾値を自動で選べました。差は論文のSMDでのギャップ(0.076)と同程度で、「ラベルなしでも実用十分なしきい値を自動決定できる」というPOTの主張が合成データでも確認できます。

スコアの時系列を見ると、赤帯(真の異常)でスコアがしっかり跳ね、POT閾値(橙)がそれらを拾える高さに自動で引かれています。

潜在次元への頑健性

最後に、潜在 $z$ の次元を変えてF1がどう動くかを見ます。

for Z in [1, 2, 3, 5, 8, 16]:
    torch.manual_seed(3)
    m = OmniVAE(M, Z=Z, use_flow=True, stochastic=True, prob=True); train(m, epochs=35)
    print("z=%2d  F1=%.3f" % (Z, best_f1(yte, m.score(Wte_t))))

$z=1,2$ ではunderfit気味でF1が0.4台に落ちますが、$z\geq 3$ では0.65前後で安定します(z=3:0.647、z=5:0.660、z=8:0.757、z=16:0.675)。論文のFig.7(z=3〜32で安定)と同じ傾向で、「小さすぎると情報が落ちるが、ある程度以上は鈍感」という性質が確認できました。だから論文も全データで z=3 という小さな値を採用しています。

OmniAnomaly と USAD の使い分け

同じ多変量時系列異常検知でも、両者は設計思想が異なります。

精度と解釈性を重視し、計算資源に余裕があるならOmniAnomaly。とにかく軽量・高速で回したいならUSAD、という住み分けです。どちらも今なお標準的なベースラインとして使われています。

まとめ

本記事では、OmniAnomalyを理論から実装、論文サーベイまで通して解説しました。

• 核心: GRUとVAEを融合した確率的RNN。時間依存と確率的ゆらぎを同時にモデル化し、再構成確率の低さで異常を検出
• 4つの工夫: 確率変数接続(潜在の時間依存)、planar NF(非ガウス事後分布)、POT(自動閾値)、次元別再構成確率(異常解釈)
• 学習: ELBO最大化。再構成項+KL正則化
• 実験: SMAP/MSL/SMDでF1 0.86、当時のSOTAを0.09上回り、適合率・再現率ともに全データで0.74超の頑健性。POTは全探索best-F1に肉薄
• 再現: 不等分散な合成データで、OmniAnomalyがGRU-AEをAUC・AP・F1で大きく上回り、POTが妥当な閾値を自動選択、z次元に頑健であることを確認

次のステップとして、以下も参考にしてください。

USAD：敵対的オートエンコーダによる多変量時系列の異常検知

OmniAnomalyを軽量・高速に置き換えた後続手法。対比で理解が深まります。

時系列の異常検知手法を体系的に解説

統計的手法からLSTM-AEまで。OmniAnomalyの位置づけが俯瞰できます。

VAEによる異常検知

再構成確率にもとづく異常検知の基礎。

参考文献

• Y. Su, R. Liu, Y. Zhao, W. Sun, C. Niu, D. Pei. “Robust Anomaly Detection for Multivariate Time Series through Stochastic Recurrent Neural Network.” Proc. 25th ACM SIGKDD (KDD ’19), 2019, pp. 2828–2837.
• J. Audibert et al. “USAD: UnSupervised Anomaly Detection on Multivariate Time Series.” KDD ’20, 2020.
• D. J. Rezende, S. Mohamed. “Variational Inference with Normalizing Flows.” ICML, 2015.
• A. Siffer et al. “Anomaly Detection in Streams with Extreme Value Theory.” KDD ’17, 2017.

さらに深く ― 確率モデルの系譜とベンチマーク

OmniAnomalyと同じく「再構成の分散（不確実性）」を扱う確率モデルや、コマンドを外生入力にする系譜は以下で深掘りしています。

• コマンド条件付け×関係ベースで「真の異常」を見分ける（2軸サーベイ）
• DC-VAE：拡張畳み込みVAEと分散推定で過渡を誤検知しない異常検知
• Telemanom：LSTM予測残差と動的閾値で異常を検知する
• ESA-ADB：実運用を映す宇宙機テレメトリ異常検知ベンチマーク