異常検知と変化検知: 9章メモ(部分空間法による変化検知)

以下の赤い本の9章を読みます。キャラクターは適当です。誤りがありましたらご指摘いただけますと幸いです。

異常検知と変化検知 (機械学習プロフェッショナルシリーズ)異常検知と変化検知 (機械学習プロフェッショナルシリーズ)
井手 剛 杉山 将

講談社 2015-08-08
売り上げランキング : 56236

Amazonで詳しく見る
by G-Tools
他の章のメモ:まだ
f:id:cookie-box:20180513082851p:plain:w60

「変化を検知したい」ってときありますよね。この本の9章を読みたいです。

f:id:cookie-box:20180513082908p:plain:w60

導入が雑…。

f:id:cookie-box:20180513082851p:plain:w60

最初の累積和法は、「正常状態の分布」「異常状態の分布」を所与として、逐次的に観測される変数が異常状態の方に近い場合にだけ異常さのポイントを積み上げる手法のようですね。…「正常状態の分布」「異常状態の分布」がわかってたら苦労しなくないですか?

f:id:cookie-box:20180513082908p:plain:w60

いやまあそこがわかってたらかなりありがたいけど、じゃあどちらの分布に近いのかというのを対数尤度比という形で定量化して、さらにそれを時系列データにおける継続的な異常状態の検知に適用するために累積和という指標にまとめる、というのも大事な仕事なんじゃないかな。

f:id:cookie-box:20180513082851p:plain:w60

そう、その「継続的な異常状態の検知」という点についてもなんですが、この累積和って異常の継続性を考慮してますか? 「たまに異常になる分には目をつぶって、継続的に異常になったら検知したい」ということで累積和をとったということですが、これって異常が連続していようと数回おきであろうと足し上げていくだけですよね。異常が継続していなくても累積和はどんどん降り積もっていきません?

f:id:cookie-box:20180513082908p:plain:w60

そういう場合は定期的にリセットするんでしょ。結局ある一定の期間内の異常さの和に閾値を置いているのと等しいとは思うよ。ある一定の期間内に多くの異常が発生することがここでの継続的な異常なんだと思う。

f:id:cookie-box:20180513082851p:plain:w60

なるほど。…次に紹介されているのは近傍法ですね。この手法で少なくとも所与とするのは「正常な部分時系列のセット」ですね。それで、検査対象の部分時系列を入力すると、最も近い正常な部分時系列までの距離(の対数の定数倍)を返すといった具合ですね。

f:id:cookie-box:20180513082908p:plain:w60

k=1 の k 近傍法で、明示的に異常標本を与えずにとにかく異常度のみ算出する場合はその手続きになるね。

f:id:cookie-box:20180513082851p:plain:w60

115ページから、ある時点(まで)のデータが異常かどうかを判定するというよりは、「検査データセット」が「参照データセット」と異質か否かという話になっていますね。その1つの抽象的な定式化が、まず「参照データの分布 p」と「検査データの分布 p'」を求めた上で、2つの分布の相違度を出すということですね。「参照データの分布」と「検査データの分布」の相違度は例えば (9.8) 式のように定義できると…対数尤度比の期待値ということですが、どうも既視感がある数式です…これは何かの予兆なんでしょうか…。

f:id:cookie-box:20180513082908p:plain:w60

いや、カルバック・ライブラー情報量ですね。

f:id:cookie-box:20180513082851p:plain:w60

はい。

f:id:cookie-box:20180513082908p:plain:w60

KL情報量は非対称だけど、これは「参照データの分布が厚いところで検査データの分布が薄いことを問題視する」って向きだね。参照データにはよく見られるようなデータが検査データに全然見られなかったらおかしいだろうと。狭い参照データ分布をカバーするように広い検査データ分布よりは、広い参照データ分布をカバーし切らない狭い検査データ分布の方が相違度が大きくなるね。例えば、
   D_{\rm KL}(p=N(0, 1^2),\,  p'=N(0, 2^2)) \approx 0.312 (参照データ分布の方が狭い場合)
   D_{\rm KL}(p=N(0, 2^2), \, p'=N(0, 1^2)) \approx 0.807 (参照データ分布の方が広い場合)

f:id:cookie-box:20180513082851p:plain:w60

参照データ分布としてそんなにだらだら広がった分布を仮定しない方がよさそうですね。検査データ分布側の密度の欠落に厳しくなりそうです。

f:id:cookie-box:20180513082908p:plain:w60

116ページからは、(9.8) 式の具体的な計算手法が1つ紹介されていて、つまり、特異スペクトル変換法だね。

f:id:cookie-box:20180513082851p:plain:w60

フォンミーゼス・フィッシャー分布…扱ったことがないですね…7章を参照するとこれは、M次元超球面上に定義される分布ということですが…なぜこんな分布を導入するのでしょう…。

f:id:cookie-box:20180513082908p:plain:w60

参照データ分布と検査データ分布にフォンミーゼス・フィッシャー分布(長いので以下vMF分布)を仮定するということだけど、(少なくともこの本の)特異スペクトル法の手順の中で「集中度」は推定せず参照データ分布も検査データ分布も適当な同一の集中度を取ると仮定しているから、結局推定対象パラメータは「平均方向」のみで、だから「参照データセットと検査データセットをそれぞれ表現するvMF分布って何だろう」というよりは、実質的には「参照データセットと検査データセットをそれぞれ代表する単位ベクトルって何だろう」ってことだね。そういう単位ベクトルを求めたらそれらを平均方向とするvMF分布どうしのKL情報量を相違度とするわけだけど。ともあれ単位ベクトルの求め方が肝心そうだね。いや…単位ベクトル自体を求めるわけじゃなくて、そのような単位ベクトルどうしの内積に相当する値を直接求めにいってるね。

f:id:cookie-box:20180513082851p:plain:w60

その手順は…参照データを横に並べた行列(M×n)と、検査データを横に並べた行列(M×k)をそれぞれ特異値分解して、それぞれからM次元の左特異ベクトルがM本ずつ出てきますが、特異値上位の何本かずつを抜き出してきて(参照データと検査データから違う本数ずつでもよい)横に並べた行列を  U, \, Q として、 U^{\top}Q の行列2ノルムこそが求める内積に相当する値であると…はて…。

f:id:cookie-box:20180513082908p:plain:w60

参照データの第i成分を全て束ねたM本のベクトルをどのようなM個の係数の組み合わせで線形結合するとそのL2ノルムが極値になるかを考える(ただし、M個の係数の2乗和は1という制約の下で)。極値を与える係数のセット(=左特異ベクトル)がMセット求まる。最も大きい極値を与える左特異ベクトルは参照データを特徴付けるベクトルだろう、って感じだね。だから参照データのそれと検査データのそれの内積を取れば相違度がわかるだろうと。ただ最も大きい極値以外無視するというのもあれだっていうなら特異値上位の何本かずつ左特異ベクトルを選ぶわけだけど、例えば2本ずつ選んだら内積が4パターン計算できちゃう。どうやってトータルの相違度とすべきかわからない。そのときは、ここでまたしても各行ベクトル(列ベクトルでもいいけど)の各成分をどのように線形結合するとL2ノルムが極値になるかという問題を解いて、最大の極値(最大特異値)をトータルの相違度とするんだね。この操作は  U^{\top}Q の行列2ノルムを取ることと等しい。

f:id:cookie-box:20180513082851p:plain:w60

うーん…そうやって求めた左特異ベクトル(たち)がその元データセットならではのベクトル(たち)というのはいいんですが、「元データセットをよく代表するベクトル(たち)です」っていうのはよくわかんないですね…。いやでも…元データセットが仮に2次元のデータたちで、どれも成分どうしが (3x, 4x) のような比率をもつデータだったとしたら、特異値最大の左特異ベクトルとしてはおそらく (3/5, 4/5) のようなものが選ばれるはずで、そういう方向との内積を取った和が最も大きくなるデータたちですよ、というのはデータセットをよく形容するベクトルのようには思えます。

f:id:cookie-box:20180513082908p:plain:w60

それが118ページでいう「一番人気のある方向」だね。

f:id:cookie-box:20180513082851p:plain:w60

ちょっと図示してみたいですね。ランダムな2次元データセットを生成して、左特異ベクトルをプロットしてみましょう。こんな感じでいいんでしょうか。灰色のベクトルたちを代表するのが赤いベクトルになると期待しているんですが…。

n <- 7
x <- c()
for (i in 1:n) {
  t <- runif(1, min=0, max=2)
  x <- c(x, 2*cos(pi*t), 2*sin(pi*t))
}
x <- matrix(x, nrow=2)

plot(c(-2,2), c(-2,2), xlab="", ylab="", type="n")
for (i in 1:n) {
  lines(c(0, x[[1,i]]), c(0, x[[2,i]]), lwd=2, col="dimgray")
}
res <- svd(x)
lines(c(0, res$u[[1,1]]), c(0, res$u[[2,1]]), lwd=4, col="red")
lines(c(0, res$u[[1,2]]), c(0, res$u[[2,2]]), lwd=4, col="darkorange")
print(res)

データセットを全て長さ2のベクトルにしているのは見やすさのためだけです。長さ1だと左特異ベクトルと重なったときに見えなくなっちゃうし、長さばらばらなのも見づらかったんで。
f:id:cookie-box:20180613230237p:plain:w240
122ページ以降は特異スペクトル変換法の高速化の話ですね。ランチョス法という手法を用いるそうです。123ページの定義9.2は…計算せよ? ああ、こういう問題を定義するというわけですか…。なるほど、この問題を効率よく計算する手段があれば、効率よく特異スペクトル変換が実現でき…ないですよね? 「参照データの左特異ベクトルとの内積を取ろう」という時点でその相手である検査データの左特異ベクトルが既にあればいいですが、それを求めるのにやっぱり特異値分解が必要ですよ?

f:id:cookie-box:20180513082908p:plain:w60

だから、特異スペクトル変換の定義を modify するみたいだね。もう検査データ側の特異値分解はしないことにして、代表的なデータを選んで(全て選んで後で相違度を何らかの方法で集約してもいいけど)左特異ベクトルの変わりにするみたいだね。

f:id:cookie-box:20180513082851p:plain:w60

「一番人気のある方向」とは何だったのか…。でもランチョス法は近似的な3重対角化なんですね。ただネットで調べるとループ計算中に a_j が数値誤差で直交しなくなっていくことに注意とかありますね。でも数値計算の授業でランチョス法ってやりましたっけ。QR法を習ったのは覚えているんですがその前後にやってたんですかね。もう10年くらい前ですから記憶がありませんね。

f:id:cookie-box:20180513082908p:plain:w60

いや部長は高校1年生だよね!? 幼稚園の年長に数値計算の授業はないよね!?

(他の章のメモがあれば)つづく