時系列データ解析システム、方法及びプログラム
【課題】 時系列データの異常度の検知精度を高めることを目的とする。
【解決手段】 テスト用時系列データと、参照用正常時系列データの各々について、先ず相関係数行列を作成する。次に、graphic lassoのアルゴリズムにより、その各々の相関係数行列から、逆行列である、疎の精度行列を作成する。精度行列が得られると、テスト用時系列データと参照用正常時系列データの各々について、好適には多変量ガウスモデルによって近傍性の確率分布が記述できる。すると、負のエントロピーの計算により、異常度が、カルバック・ライブラー距離として計算できる。この技法は、近傍性保存原理に基づく時系列データの異常検出手法に、各センサの周りの近接グラフ構造をデータから自動的に推測する機能が組み込まれ、さらには、局所的な統計モデルの直接比較に基づいて、理論的に首尾一貫した異常度を与える、という点で有利である。
【解決手段】 テスト用時系列データと、参照用正常時系列データの各々について、先ず相関係数行列を作成する。次に、graphic lassoのアルゴリズムにより、その各々の相関係数行列から、逆行列である、疎の精度行列を作成する。精度行列が得られると、テスト用時系列データと参照用正常時系列データの各々について、好適には多変量ガウスモデルによって近傍性の確率分布が記述できる。すると、負のエントロピーの計算により、異常度が、カルバック・ライブラー距離として計算できる。この技法は、近傍性保存原理に基づく時系列データの異常検出手法に、各センサの周りの近接グラフ構造をデータから自動的に推測する機能が組み込まれ、さらには、局所的な統計モデルの直接比較に基づいて、理論的に首尾一貫した異常度を与える、という点で有利である。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、センサなどから取得された時系列データを解析するためのシステム、方法及びプログラムに関する。
【背景技術】
【０００２】
従来より、製造装置や監視装置などに取り付けたセンサから取得した時系列データを監視することにより、異常の有無を監視するという処理が行われている。
【０００３】
一方、近年、センサを取り付けた自動車を走行させ、そのようなセンサから得られた時系列データを解析する試みも行われている。自動車は、エンジンの冷却水の温度センサ、エンジンの吸入空気の温度センサ、オイル温度センサ、燃料噴射装置用の吸気管内圧力センサ、ターボチャージャ用の過給圧センサ、スロットルポジョンセンサ、ステアリング舵角センサ、車高センサ、液面センサ、電磁による回転速度センサ、ノックセンサ、加速度センサ、流量センサ、酸素センサ、希薄空燃比センサなど多くのセンサを有し、最新の自動車では、センサの数は、数百にも達する。
【０００４】
このような複数のセンサの時系列出力データを用いて行いたいことの１つに、異常診断がある。すなわち、図１(b)に示すのが、そのような異常診断を行いたいテスト・データとする。このようなテストデータが、
Ｄ = {x(t)|x(t)∈Ｒp, t = 1,2,...,N}で表されるとする。ここで、Ｒpは、p次元実数空間をあらわし、よって、x(t)は、実数成分をもつp次元ベクトルである。その１つの成分が、１つのセンサに対応する。tは、時間のインデックスで、ある増分毎の時間に対応する。よって、x(1),x(2),...,x(N)は、時系列に沿った、N個のベクトルである。
【０００５】
一方、参照用正常データが、テストデータと同一フォーマットで用意されているとする。これは、例えば、自動車の正常走行中に測定したセンサの値とする。参照用正常データは、
Ｄ〜 = {x~(t)|x~(t)∈Ｒp, t = 1,2,...,N}で表されるとする。
【０００６】
すると、異常診断を行うということは、Ｄ〜に対して、Ｄを比較して、どのセンサからの信号が異常を示しているかを判別する、ということである。
【０００７】
このような比較を行うためのアプローチとして、従来より知られている方法には、下記のようなものがある。
(1) 時間相関が強い極限からのアプローチとして、自己回帰モデルのような時系列モデリングを利用する技法、例えば、特開平６−１０９４９８号公報に開示されているような技法である。
(2) 独立同一分布の仮定に基づく統計解析の手法を利用する技法。
(3) システムの運動法則に基づくカオス時系列解析などを用いる技法であり、例えば、特開平８−２７８８１５号公報に開示されているような技法である。
【０００８】
しかし、特に自動車のセンサのデータの場合、(a) 非常に動的であり、(b) 時系列信号同士にしばしば強い相関がある、(c) センサ毎に事前知識は与えられない、という特徴がある。すなわち、信号が非常に動的で複雑であるため、データが異なれば、時系列の様相ががらりと変わる。従って、同一センサからの信号であっても、「重ねて比べてみる」ことが意味をなさない。つまり、時系列からの定常性が仮定できない。このような理由から、伝統的な時系列予測の手法を用いて時系列予測モデルを作り、予測からのずれの検出という技法で問題を解くことが難しい。
【０００９】
とはいえ、センサのデータは、独立同一分布に従うとはみなせず、しかも、その確率分布は一般に、複雑なものになる。さらに、センサ間には一般に依存関係があり、各センサを多変量系として扱わざるを得ない。ここで依存関係というのは、例えば、「アクセルを踏むと加速する」という関係で、これは、スロットルポジョンセンサからの出力と、回転速度センサや加速度センサからの出力とが関連している、というようなことである。ところが、多変量系では一般に、統計的漸近論の成立に必要なデータ数が、センサの個数であるpに指数的に依存する。従って、N→∞の漸近論が、多くの現実的な状況では使えない。
【００１０】
一方、自動車、発電設備などの場合、センサを設置される対象は非常に複雑な物理系となっており、システム全体の運動方程式を立てることは、ほぼ絶望である。
【００１１】
このような事情から、上記従来技術を、自動車のセンサのデータのような多変量系と見なされるデータに適用することはできない。
【００１２】
そこで本願発明者の一人を含む発明者らは、本出願に譲渡された、２００７年１月３０日出願の米国特許出願第１１／６６８７４５号明細書において、センサ信号の相関異常の解明を行う発明として、相関行列を元に各センサを、多次元尺度構成法でユークリッド空間に埋め込み、埋め込まれた空間での斥力ポテンシャルを異常度と結びつける方法を提案した。しかし、この技法は、ノイズ除去機能が十分でなく、埋め込み過程が不安定になりがちで、また、算出される異常度が規格化されていない、という点で、上述の問題に当てはめるに十分でない。
【００１３】
さらに、本願発明者の一人を含む発明者らは、本出願に譲渡された、２００７年１２月１８日出願の米国特許出願第１１／９５９０７３号明細書において、近傍性保存原理というアイデアに基づいて、上記の困難を解決することを意図した異常度検出手法を提案した。この手法の特徴は、センサ信号同士の類似度の破れから異常を見出すこと、各センサについて、近傍にあるセンサ以外は考慮しないこと、及び、近傍グラフのタイトさの変化から、異常度を計算することである。
【００１４】
この技法は、自動車データの解析や中規模の発電設備の異常検知に利用され、実用的な成功を収めている。しかし、この技法の限界は、センサの近傍をうまく見出す手段をもたないことである。例えば、他のどれとも独立であるような変数では、近傍と呼べる仲間はいない。一方、測定系が多重化されているようなセンサでは、その多重度に応じた仲間を見出せるはずである。しかし、事前にそのようなセンサ個別の相違を解析する手段は未知であり、変数毎の個性を無視して、一律に近傍数ｋを決めるしかなかった。そのため、とりわけ他と独立性が強いようなセンサで擬陽性（問題でないのに、問題とみなさる）が多発するという問題があった。
【００１５】
さらに、異常度の定義がヒューリスティックスに基づいていたため、単一のセンサの陽性反転に基づくエラーが検出できないとか、ノイズが多い状況では検知能力が十分でない、という問題があった。
【００１６】
【特許文献１】特開平６−１０９４９８号公報
【特許文献２】特開平８−２７８８１５号公報
【特許文献３】２００７年１月３０日出願の米国特許出願第１１／６６８７４５号明細書
【特許文献４】２００７年１２月１８日出願の米国特許出願第１１／９５９０７３号明細書
【非特許文献１】J. Friedman, T. Hastie, and R. Tibahirani, Sparse inverse covariance estimation with the graphical lasso. Biostatics, 9(3): 432-441, 2008
【非特許文献２】O. Banarjee, L.E.Ghaoui, and G. Natsoulis, Convex optimization techniques for fitting sparse gaussian graphical models. In in Proceedings of ICML, pages 89-96. Press, 2006.
【非特許文献３】J. Friedman, T. Hastie, H. Hofling and Tibshinari, Pathwise coordinate optimization. Annals of Applied Statistics. I(2):302-332
【発明の開示】
【発明が解決しようとする課題】
【００１７】
この発明の目的は、各センサの周りの近接グラフ構造をデータから自動的に推測する機能を提供することによって、上記従来技術の問題点を解決することにある。
【課題を解決するための手段】
【００１８】
本発明において、入力として与えられるのは、下記のものである。
(1) テスト・データＤ = {x(t)|x(t)∈Ｒp, t = 1,2,...,N}
(2) 正規化定数ρ ここで、0 < ρ < 1であり、その値は、ユーザが適宜選択する。
(3) 参照用正常データＤ〜 = {x~(t)|x~(t)∈Ｒp, t = 1,2,...,N}の相関係数行列Ｓ〜
【００１９】
これによって出力されるのは、各センサの異常度e1, e2,...,epである。
【００２０】
本発明に従うアルゴリズムの概要は、次のとおりである。以下の処理はすべて、コンピュータのプログラムによって、自動的に実行されることに留意されたい。まず、corr()という関数表示で、相関係数行列を求める手続きをあらわす。すると、最初のステップは、下記のとおり、テスト・データＤの相関係数行列Ｓを求めることである。
１．Ｓ ← corr(Ｄ)
２．次のステップは、Ｄ、及びＤ〜に対する精度行列の計算である。精度行列とは基本的に、相関係数行列の逆行列であるが、特に本発明においては、L1ノルムの正規化項を付加した式を最適化するように解くことによって、精度行列が疎行列になるような解を得る。このような最適化の式を解くために好適なアルゴリズムの１つに、graphical lassoがある。それを、glasso()という関数表示であらわすと、
Λ ← glasso(ρ,Ｓ)、Λ〜 ← glasso(ρ,Ｓ〜)で、テスト・データＤの精度行列Λと、参照用正常データＤ〜の精度行列Λ〜を得る。
３，次のステップは、i = 1,2,...,pに対する特徴量行列の計算である。特徴量行列は、次のようにあらわされる。
【数１】
但し、li(A|B)というのは、negentropy(A,B)のことである。
４．次のステップは、下記の式によって、i = 1,2,...,pに対する異常度を計算することである。
ei = max{li(Λ|Ｓ) - li(Λ〜|Ｓ),li(Λ|Ｓ〜) - li(Λ〜|Ｓ〜)}
これは、相対エントロピー、またはカルバック・ライブラー距離として知られているものである。
【発明の効果】
【００２１】
以上のように、この発明によれば、精度行列が疎行列になるような処理の工夫を施して、好適にはgraphic lassoを用いてこの最適化の式を解き、各センサについて、カルバック・ライブラー距離が異常度をあらわすようにしたことによって、近傍性保存原理に基づく時系列データの異常検出手法に、各センサの周りの近接グラフ構造をデータから自動的に推測する機能が組み込まれ、さらには、局所的な統計モデルの直接比較に基づいて、理論的に首尾一貫した異常度が与えられる。
【００２２】
これによって、擬陽性の少ない、より適切なセンサ時系列データ異常判別技法が提供される。
【発明を実施するための最良の形態】
【００２３】
以下、図面に従って、本発明の実施例を説明する。これらの実施例は、本発明の好適な態様を説明するためのものであり、発明の範囲をここで示すものに限定する意図はないことを理解されたい。また、以下の図を通して、特に断わらない限り、同一符号は、同一の対象を指すものとする。
【００２４】
図２を参照すると、本発明の一実施例に係るシステム構成及び処理を実現するためのコンピュータ・ハードウェアのブロック図が示されている。図２において、システム・バス２０２には、ＣＰＵ２０４と、主記憶（ＲＡＭ）２０６と、ハードディスク・ドライブ（ＨＤＤ）２０８と、キーボード２１０と、マウス２１２と、ディスプレイ２１４が接続されている。ＣＰＵ１０４は、好適には、３２ビットまたは６４ビットのアーキテクチャに基づくものであり、例えば、インテル社のＰｅｎｔｉｕｍ（商標）４、インテル社のＣｏｒｅ（商標） ２ ＤＵＯ、ＡＭＤ社のＡｔｈｌｏｎ（商標）などを使用することができる。主記憶２０４は、好適には、５１２ＫＢ以上の容量、より好ましくは、１ＧＢ以上の容量をもつものである。
【００２５】
ハードディスク・ドライブ２０８には、個々に図示しないが、オペレーティング・システム及び本発明に係る処理プログラムなどが、予め格納されている。オペレーティング・システムは、Ｌｉｎｕｘ（商標）、マイクロソフト社のＷｉｎｄｏｗｓ Ｖｉｓｔａ、Ｗｉｎｄｏｗｓ ＸＰ（商標）、Ｗｉｎｄｏｗｓ（商標）２０００、アップルコンピュータのＭａｃ ＯＳ（商標）などの、ＣＰＵ２０４に適合する任意のものでよい。
【００２６】
ハードディスク・ドライブ２０８にはさらに、センサの基準用の正常の時系列データと、センサのテスト用時系列データとが個別のファイルとして格納されている。異なる状況で記録された複数のデータが、個別のファイルとして保存されていてよい。すると、ユーザは、テスト時に、保存されているファイルのうちの１つを選択することによって、所望の時系列データの異常検知を行うことができる。
【００２７】
同様に、基準用の正常の時系列データのファイルも、複数用意してハードディスク・ドライブ２０８に保存していてもよい。
【００２８】
キーボード２１０及びマウス２１２は、オペレーティング・システムが提供するグラフィック・ユーザ・インターフェースに従い、ディスプレイ２１４に表示されたアイコン、タスクバー、ウインドウなどのグラフィック・オブジェクトを操作するために使用される。キーボード２１０及びマウス２１２はまた、異常検知のためのテスト用センサ時系列データが格納されたファイルを指定するためにも使用される。
【００２９】
ディスプレイ２１４は、これには限定されないが、好適には、１０２４×７６８以上の解像度をもち、３２ビットtrue colorのＬＣＤモニタである。ディスプレイ２１４は、時系列データの波形や、異常検知の結果を表示したりするために使用される。
【００３０】
ハードディスク・ドライブ２０８にはさらに、本発明に係る異常検知を行うためのプログラムが格納されている。このプログラムは、Ｃ＋＋、Ｃ＃、Ｊａｖａ（商標）、Ｐｅｒｌ、Ｒｕｂｙなどの既存の任意のプログラム言語で書くことができる。オペレーティング・システムとして、Ｗｉｎｄｏｗｓ Ｖｉｓｔａ、Ｗｉｎｄｏｗｓ ＸＰ（商標）、Ｗｉｎｄｏｗｓ（商標）２０００などを使用する場合には、Ｗｉｎ３２ ＡＰＩの機能を利用して、ＧＵＩも含むアプリケーション・プログラムとして実装することができる。しかし、本発明に係る異常検知を行うためのプログラムは、ＣＵＩとしても実装することが可能である。
【００３１】
図３は、センサからのデータの間の関連性を示す図である。図示されているように、本発明の異常検知の仕組みでは、iを1からpまでの間の任意の整数としたとき、センサからのデータxiにつき、データxi以外のデータx1, ..., xj, ...,xi-1,xi+1,...,xpとの間の関係をより際立たせるようなアルゴリズムが適用される。すなわち、弱い関係同士のデータの関連性は無視してしまい、かなり関係の大きいデータ同士の関係性のみが残されるようになされる。このような関係性を行列表示すると、ほとんどの成分が0の疎行列になる。換言すると、このような関係性をあらわす行列である、精度行列が、疎行列であるような解を求めることが、本発明の１つの目標である。図３では、データxiに対して、xjとはある程度関係があり、xi+1とは強い関係があり、それ以外の点線であらわす関係は、無関係とみなすことを示している。
【００３２】
なお、精度行列については、C.M. Bishop, "Pattern Matching and Machine Learning". Springer Verlagの2.3.1節などを参照されたい。
【００３３】
図４は、本発明に係る処理の概要を示すフローチャートである。このフローチャートを実行するプログラムは、ハードディスク・ドライブ２０８に保存され、ユーザの操作により、オペレーティング・システムの働きで、主記憶２０６にロードされて実行可能となる。
【００３４】
ステップ４０２では、手続きcorr()によって、テスト・データＤの相関係数行列が計算される。
テスト・データＤの定義を再掲すると、 {x(t)|x(t)∈Ｒp, t = 1,2,...,N}である。ここで、時系列データの次元数pは、例えば50程度で、数百になることもある。Nは、時分割の数で、一般的には100から1000程度である。時分割の刻みの適切な間隔は、適用例によって異なるが、ここでは0.1秒である。但し、後の計算量に効いてくるのは、Nよりも寧ろ、次元数pである。というのは、pが直接、相関係数行列の行数及び列数になるからである。このデータは、ハードディスク・ドライブ２０８に保存されており、処理のため、本発明の処理プログラムにより、主記憶２０６に読み出される。
【００３５】
テスト・データＤの相関係数行列は、次のようにして計算される。先ず、x(t)のi番目の成分を、xi(t)とあらわすことにする。ここで、i∈{1,..,p}である。
【００３６】
そこで、処理プログラムは、xi(t)の、t = 1,2,...,Nについての平均miと標準偏差σiとを計算する。そして、処理プログラムは、平均miと標準偏差σiを用いて、次の式により、xi(t)を標準化する。
【数２】
【００３７】
こうしておいて、処理プログラムは、次の式により、相関係数行列Ｓを求める。
【数３】
【００３８】
なお、ここでは参照用正常データＤ〜 = {x~(t)|x~(t)∈Ｒp, t = 1,2,...,N}の相関係数行列Ｓ〜は、所与であると想定しているが、もしまだ計算されていないなら、ここで同様に計算する。
【００３９】
ステップ４０４で、処理プログラムは、下記のように、関数glasso()を用いて、テスト・データＤの精度行列Λと、参照用正常データＤ〜の精度行列Λ〜とを計算する。glassoとは、graphical lassoアルゴリズムを意味する。
Λ = glasso(ρ,Ｄ)
Λ〜 = glasso(ρ,Ｄ〜)
【００４０】
精度行列Λとは、定義によれば、相関係数行列のSの逆行列である。しかし、Ｓの逆行列を計算することは一般に不可能である。なぜならＳは一般に正則ではないからである。これは例えば、ある程度強い相関がある変数の対が存在する場合にそうなる。センサー系の時系列データではそのような状況は常に起こっていると考えなければならない。仮に、幸運にもＳが正則だったとしても、単純に逆行列を計算して精度行列Λを求めればよいかというと否である。それでは、精度行列Λが疎行列にならず、近傍を自動で選択するという本発明の狙いを満たさないためである。そのための工夫は後で説明するとして、疎行列である精度行列Λが求まったとしよう。すると、グラフィカルモデリングの理論によれば、精度行列Λの(i,j)成分(i≠j)が0であれば、xiとxjは、これら以外の全ての変数を与えたときに条件付独立である。すなわち、
【数４】
【００４１】
本発明は、この条件付き独立性を、近傍探索に利用する。すなわち、統計的に独立でないと判断されれば近傍とみなし、そうでなければ近傍でないとする。
例えば、x2について、精度行列Λ2,j ≠ 0となるjが4と9だけであれば、x2の近傍がx4とx9だけ、ということが分かる。従って、多くの行列要素がゼロであるような疎行列として精度行列Λを最適に決定することは、近傍の自動選択と等価である。
【００４２】
さて、前にも述べたように、センサ・データ自体が多変量正規分布に従うという保証はない。そこで寧ろ、近傍か否かを決める基準として、多変量ガウスモデルを援用することにする。そこで今は、平均ゼロを前提としているので、データに当てはめようとしているモデルは、平均ゼロ、精度行列Λのp次元ガウス分布である。それは、下記の式のようにあらわされる。
【数５】
ここで、detは行列式をあらわす。
【００４３】
ここで、特定の(i,j)、例えば、(i,j) = (1,2)に着目する。そして、x1とx2以外に適当な数字として、例えば、0をいれてみる。すると、指数関数の肩は、次のようになる。
【数６】
【００４４】
従って、Λ1,2 = 0であれば、x1とx2を結ぶ交差項はなく、すると、それらは、他の全てのデータを止めた時に独立となることが分かる。
【００４５】
そこで、精度行列Λの疎性（sparsity）の議論に移る。疎性を一旦離れて、精度行列をデータから素朴に最尤推定することを試みる。すると、テスト・データＤの対数尤度は、次のようになる。
【数７】
【００４６】
ここで、trは行列の対角和である。従って、尤度を最大にする精度行列を求めるためには、上記の式の{ }内を最大化すればよいことになる。
【００４７】
しかし生憎、このようにして求めた精度行列は、一般に疎にならない。そこで、本発明は、L1ノルムに基づく正規化項を付与して最尤推定すると疎行列が得られるという性質を利用する。このようにL1正規化項を付して疎な解を得る手法が統計学の分野でlasso （Least absolute shrinkage and selection operator）と呼ばれることが、ここで使うgraphical lassoの名前の所以である。
【００４８】
そこで、L1ノルムを付与した関数f(Λ;Ｓ,ρ)を次のように定義する。
【数８】
ここで、|Λi,j|は、行列Λの(i,j)要素の絶対値をあらわす。
【００４９】
すると、この関数f(Λ;Ｓ,ρ)を用いて、Λを求めることは、下記ような最適化問題を解くことに帰着される。
【数９】
【００５０】
この最適化問題を解くための好適な方法は、graphical lassoというアルゴリズムである。これの詳しい説明は、J. Friedman, T. Hastie, and R. Tibahirani, Sparse inverse covariance estimation with the graphical lasso. Biostatics, 9(3): 432-441, 2008 などにある。
【００５１】
正規化項
【数１０】
の存在によって、Λの多くの行列要素は、正定値性を失わない範囲で厳密にゼロになる。定数ρは、直感的に言うと、どの程度相関係数を無視するかの尺度であり、0 < ρ < 1の範囲でユーザが与える。graphical lassoの場合の好適な１つの値は、ρ = 0.3である。ρが1に近いと、得られる近接グラフは非常に小さくなり、ρ = 0だと、基本的に全てのセンサが結合した完全グラフが得られる。
【００５２】
なお、graphical lassoのアルゴリズムについては、後で、もう少し詳しく説明する。
【００５３】
このようにして、ステップ４０４で、テスト・データＤの精度行列Λと、参照用正常データＤ〜の精度行列Λ〜とが求められると、ステップ４０６での負のエントロピー（ネゲントロピー）の行列の計算と、ステップ４０８での異常度スコアの計算に移行する。ステップ４０６及びステップ４０８は、i = 1,2,...,pについて順次実行される。
【００５４】
ところで、テスト・データＤと、参照用正常データＤ〜との間に本質的な違いがなければ、それぞれから得られる精度行列ΛとΛ〜とが表す近接グラフの様子には大差はないはずである。すなわち、精度行列を求めるためのアルゴリズムに鑑みると、近接グラフ自体が、図１に示すような時系列データの見かけ上の相違とは関係なく求まるはずである。逆に言えば、ある変数xiに着目したとき、その周りの近接グラフに大きな変化がおきていれば、その変数に何か異常が生じている公算が高い。
【００５５】
さて、いまやテスト・データＤの精度行列Λと、参照用正常データＤ〜の精度行列Λ〜とが得られているので、
テスト・データＤの確率モデルp(x) = N(x|0,Λ-1)と、
テスト・データＤ〜の確率モデルp~(x) = N(x|0,Λ〜-1)とが得られる。そこで、これらを勘案して、変数xiの異常度を計算するに当たっては、xi以外の全ての変数x1,x2,...,xi-1, xi+1,...,xpを与えたときのxiの条件付き分布p(xi|rest)と、p~(xi|rest)の違いを評価するのが自然である。ここで、restというのは、xi以外の全ての変数x1,x2,...,xi-1, xi+1,...,xpのことである。そこで、記法の便宜上、xi以外の全ての変数x1,x2,...,xi-1, xi+1,...,xpをzi∈Ｒp-1とすると、ネゲントロピーの関数は、下記のような式となる。
【数１１】
ここで、< .. >Dという表示は、データＤでの経験分布をあらわす。
これらの式が、前述のガウス分布の式から見導かれることは、この分野の当業者に明らかであろう。これらの式に基づき、ステップ４０６で、個々のi=1,...,pにつき、負のエントロピーが計算される。
【００５６】
これらの結果に基づき、下記の式を用いて、処理プログラムは、i=1,2,...,pにつき、異常度eiの計算を行う。
【数１２】
この結果は、見て取れるように、相対エントロピー、すなわちカルバック・ライブラーの距離である。
そして、ei ≡ max(d1,d2)により、i番目のデータの異常度が求められる。なお、前述の２００７年１２月１８日出願の米国特許出願第１１／９５９０７３号明細書に記載の技法では、相違度行列Ｄの要素di,jが、- ln|Si,j|で定義されていたので、相関係数の符合反転を検知することができなかったが、本発明の技法には、そのような制約はない。
【００５７】
いくつかの実験で、本発明の技法は、ノイズが大きく、いくつかのセンサに断線が生じているような状況でも、高い精度で異常度が検知できることが分かった。これは、疎な精度行列としてテスト・データＤと、参照用正常データＤ〜を比較した場合、ノイズや断線の効果が、かなりの程度抑えられるからであろうと考察される。
【００５８】
＜grahical lassoの計算に関する補足説明＞
前述のように、grahical lassoにより解くべき問題は、下記の式を最大にする行列Λを求めることであった。
【数１３】
そこで、この式をΛで偏微分する。このとき、下記の公式に留意する。
【数１４】
【００５９】
すると、Λでの偏微分の結果の勾配の式は、下記のようになる。
【数１５】
ここで、sign(Λ)とは、符号関数sign()を、行列の要素に適用した結果を要素とする行列である。すなわち、行列sign(Λ)の(i,j)要素は、Λi,j≠0のときsign(Λi,j)であり、Λi,j=0であるなら、[-1,1]の何らかの値を返すように定義される。
【００６０】
上記勾配の式をゼロ行列と等値してΛを求めることになるが、Λはその性質から、対角要素は全て正であり、よって、sign(Λ)の対角要素は全て1となる。
【００６１】
これから、Σi,j = Si,j + ρと、対角要素は厳密に求まる。ここで、Σ≡Λ-1と定義した。
【００６２】
以上のように対角要素は厳密に求まるが、非対角要素は簡単には求まらない。そこで、ブロック勾配法という技法が使われる。ブロック勾配法については、"O. Banarjee, L.E.Ghaoui, and G. Natsoulis, Convex optimization techniques for fitting sparse gaussian graphical models. In in Proceedings of ICML, pages 89-96. Press, 2006."及び"J. Friedman, T. Hastie, H. Hofling and Tibshinari, Pathwise coordinate optimization. Annals of Applied Statistics. I(2):302-332"などにも記載がある。
ブロック勾配法では、次のように行列が分解される。
【数１６】
【００６３】
ここで、Ｌ及びＷは、(p-1)×(p-1)の行列であり、λとσはスカラーである。従って、wとlは、(p-1)次元のベクトル、ということになる。
【００６４】
ブロック勾配法は、ここでは、一つの列ベクトルであるwとlを未知、他のλ、σ、Ｌ、Ｗを既知として最適化問題を解くものである。
【００６５】
このように決めたwに対して、数１５がゼロになるという条件は、下記の式と等価である。βは、あるベクトル・パラメータ変数である。
【数１７】
但し、
【数１８】
である。また、sは、相関係数行列Ｓにおける、wに対応する位置のベクトルである。この式から見て取れるように、βは、L1正規化項付きの２次計画問題を解くことで求まる。よってβは、一般的な２次計画ソルバーを使っても解けるが、好適な解法があり、それについては、後述する。
【００６６】
数１５がゼロに等しいという条件の、Σ≡Λ-1におけるwに対応する位置のベクトルでの部分の式は、w - s - sign(l) = 0と書ける。
【００６７】
そこで、ΣΛを計算してみると、次のようになる。
【数１９】
この式の右上のブロックから得られる l = -λＷ-1wを使い、また、β≡Ｗ-1wとおけば、
Ｗβ - s +ρsign(β) = 0
数１７の偏微分を実行することによって、この式に等しくなることが確認できる。
【００６８】
上の説明では、第p変数についてブロック勾配法を適用するものであったが、より一般に次のようなステップであると言える。
(1) 最初に例えば、Σ = Ｓ + ρ１pと置く。
(2) i = 1から始めて、第i変数について、行と列を並べ替えて、i行目が第p行目、i列目が第p列目に来るようにする。
(3) lとλを最適に決めた後、行と列の配列を戻す。
(4) 元の行列と違いが小さければ終了。そうでなければ、iを１つ増やす。iがpを超えたなら、iを1にする。
(5) ステップ(2)に戻る。
【００６９】
このとき、lとλを最適に決めるステップであるが、数１９の右下部分から、wTl+σλ=1が得られ、これと、l = -λＷ-1wとから、下記の式が得られる。
【数２０】
【００７０】
ここで、結果的に使われる式には、Ｗ-1が陽に表れていないことに留意されたい。すなわち、相関係数行列がランク落ちしている場合、一般的にＷ-1を求めることができない。その意味で、逆行列を使わないこの計算方法は、ランク落ちしている相関係数行列にも対応できる点で、有利である。
【００７１】
さて、数１７に対応する最適化問題を書き下す。
【数２１】
この目的関数をg(β;Ｗ,ρ)と書き、βiで偏微分すると、
【数２２】
が成り立つ。これに基づく最適化条件を次の２つの場合に分けて考える。
【００７２】
(1) βi > 0のとき、数２２が0という条件は、βm (m≠i)を与えた下で、
【数２３】
と解ける。但し、
【数２４】
と置いた。数２３で、Ai - ρ < 0 であれば、もともとの条件βi > 0と矛盾する。このときは、
【数２５】
のように単調であることから、定義域の左端、すなわち、βi → 0が解となる。
【００７３】
(2) βi < 0のとき、数２２が0という条件は、
【数２６】
と、解ける。但し、Ai + ρ > 0であれば、もともとの条件βi < 0と矛盾する。このときは、
【数２７】
のように単調であることから、定義域の左端、すなわち、βi → 0が解となる。
【００７４】
以上、２つの場合を纏めると、数２１の最適化問題は、p-1個のiそれぞれについて、
【数２８】
という置換を繰り返せばよい。添え字iが最後の変数まで行ったら、また最初に戻って、収束するまで繰り返す。
【００７５】
なお、Σi,j = Si,j + ρという式から見てとれるように、Ｗi,iは必ず正であり、この解に特異性はない。また、上述したgraphical lassoの２つの記述法は、簡単な反復公式で記述でき、計算効率が高い。らに、前述のように、Σがランク落ちしていも、その逆行列を計算に使う必要がないので、安定して計算を実行することができる。
【００７６】
以上の例は、自動車のセンサから取得したデータを異常検知する場合であったが、これには限定されず、発電所や、プラントなど、測定データが多次元時系列データとあらわされる、任意の場合に適用可能である。
【図面の簡単な説明】
【００７７】
【図１】参照用正常時系列データと、テスト用時系列データの波形の例を示す図である。
【図２】本発明を実施するためのハードウェアのブロック図である。
【図３】変数間の近傍関係を示す図である。
【図４】本発明の処理のフローチャートを示す図である。
【技術分野】
【０００１】
この発明は、センサなどから取得された時系列データを解析するためのシステム、方法及びプログラムに関する。
【背景技術】
【０００２】
従来より、製造装置や監視装置などに取り付けたセンサから取得した時系列データを監視することにより、異常の有無を監視するという処理が行われている。
【０００３】
一方、近年、センサを取り付けた自動車を走行させ、そのようなセンサから得られた時系列データを解析する試みも行われている。自動車は、エンジンの冷却水の温度センサ、エンジンの吸入空気の温度センサ、オイル温度センサ、燃料噴射装置用の吸気管内圧力センサ、ターボチャージャ用の過給圧センサ、スロットルポジョンセンサ、ステアリング舵角センサ、車高センサ、液面センサ、電磁による回転速度センサ、ノックセンサ、加速度センサ、流量センサ、酸素センサ、希薄空燃比センサなど多くのセンサを有し、最新の自動車では、センサの数は、数百にも達する。
【０００４】
このような複数のセンサの時系列出力データを用いて行いたいことの１つに、異常診断がある。すなわち、図１(b)に示すのが、そのような異常診断を行いたいテスト・データとする。このようなテストデータが、
Ｄ = {x(t)|x(t)∈Ｒp, t = 1,2,...,N}で表されるとする。ここで、Ｒpは、p次元実数空間をあらわし、よって、x(t)は、実数成分をもつp次元ベクトルである。その１つの成分が、１つのセンサに対応する。tは、時間のインデックスで、ある増分毎の時間に対応する。よって、x(1),x(2),...,x(N)は、時系列に沿った、N個のベクトルである。
【０００５】
一方、参照用正常データが、テストデータと同一フォーマットで用意されているとする。これは、例えば、自動車の正常走行中に測定したセンサの値とする。参照用正常データは、
Ｄ〜 = {x~(t)|x~(t)∈Ｒp, t = 1,2,...,N}で表されるとする。
【０００６】
すると、異常診断を行うということは、Ｄ〜に対して、Ｄを比較して、どのセンサからの信号が異常を示しているかを判別する、ということである。
【０００７】
このような比較を行うためのアプローチとして、従来より知られている方法には、下記のようなものがある。
(1) 時間相関が強い極限からのアプローチとして、自己回帰モデルのような時系列モデリングを利用する技法、例えば、特開平６−１０９４９８号公報に開示されているような技法である。
(2) 独立同一分布の仮定に基づく統計解析の手法を利用する技法。
(3) システムの運動法則に基づくカオス時系列解析などを用いる技法であり、例えば、特開平８−２７８８１５号公報に開示されているような技法である。
【０００８】
しかし、特に自動車のセンサのデータの場合、(a) 非常に動的であり、(b) 時系列信号同士にしばしば強い相関がある、(c) センサ毎に事前知識は与えられない、という特徴がある。すなわち、信号が非常に動的で複雑であるため、データが異なれば、時系列の様相ががらりと変わる。従って、同一センサからの信号であっても、「重ねて比べてみる」ことが意味をなさない。つまり、時系列からの定常性が仮定できない。このような理由から、伝統的な時系列予測の手法を用いて時系列予測モデルを作り、予測からのずれの検出という技法で問題を解くことが難しい。
【０００９】
とはいえ、センサのデータは、独立同一分布に従うとはみなせず、しかも、その確率分布は一般に、複雑なものになる。さらに、センサ間には一般に依存関係があり、各センサを多変量系として扱わざるを得ない。ここで依存関係というのは、例えば、「アクセルを踏むと加速する」という関係で、これは、スロットルポジョンセンサからの出力と、回転速度センサや加速度センサからの出力とが関連している、というようなことである。ところが、多変量系では一般に、統計的漸近論の成立に必要なデータ数が、センサの個数であるpに指数的に依存する。従って、N→∞の漸近論が、多くの現実的な状況では使えない。
【００１０】
一方、自動車、発電設備などの場合、センサを設置される対象は非常に複雑な物理系となっており、システム全体の運動方程式を立てることは、ほぼ絶望である。
【００１１】
このような事情から、上記従来技術を、自動車のセンサのデータのような多変量系と見なされるデータに適用することはできない。
【００１２】
そこで本願発明者の一人を含む発明者らは、本出願に譲渡された、２００７年１月３０日出願の米国特許出願第１１／６６８７４５号明細書において、センサ信号の相関異常の解明を行う発明として、相関行列を元に各センサを、多次元尺度構成法でユークリッド空間に埋め込み、埋め込まれた空間での斥力ポテンシャルを異常度と結びつける方法を提案した。しかし、この技法は、ノイズ除去機能が十分でなく、埋め込み過程が不安定になりがちで、また、算出される異常度が規格化されていない、という点で、上述の問題に当てはめるに十分でない。
【００１３】
さらに、本願発明者の一人を含む発明者らは、本出願に譲渡された、２００７年１２月１８日出願の米国特許出願第１１／９５９０７３号明細書において、近傍性保存原理というアイデアに基づいて、上記の困難を解決することを意図した異常度検出手法を提案した。この手法の特徴は、センサ信号同士の類似度の破れから異常を見出すこと、各センサについて、近傍にあるセンサ以外は考慮しないこと、及び、近傍グラフのタイトさの変化から、異常度を計算することである。
【００１４】
この技法は、自動車データの解析や中規模の発電設備の異常検知に利用され、実用的な成功を収めている。しかし、この技法の限界は、センサの近傍をうまく見出す手段をもたないことである。例えば、他のどれとも独立であるような変数では、近傍と呼べる仲間はいない。一方、測定系が多重化されているようなセンサでは、その多重度に応じた仲間を見出せるはずである。しかし、事前にそのようなセンサ個別の相違を解析する手段は未知であり、変数毎の個性を無視して、一律に近傍数ｋを決めるしかなかった。そのため、とりわけ他と独立性が強いようなセンサで擬陽性（問題でないのに、問題とみなさる）が多発するという問題があった。
【００１５】
さらに、異常度の定義がヒューリスティックスに基づいていたため、単一のセンサの陽性反転に基づくエラーが検出できないとか、ノイズが多い状況では検知能力が十分でない、という問題があった。
【００１６】
【特許文献１】特開平６−１０９４９８号公報
【特許文献２】特開平８−２７８８１５号公報
【特許文献３】２００７年１月３０日出願の米国特許出願第１１／６６８７４５号明細書
【特許文献４】２００７年１２月１８日出願の米国特許出願第１１／９５９０７３号明細書
【非特許文献１】J. Friedman, T. Hastie, and R. Tibahirani, Sparse inverse covariance estimation with the graphical lasso. Biostatics, 9(3): 432-441, 2008
【非特許文献２】O. Banarjee, L.E.Ghaoui, and G. Natsoulis, Convex optimization techniques for fitting sparse gaussian graphical models. In in Proceedings of ICML, pages 89-96. Press, 2006.
【非特許文献３】J. Friedman, T. Hastie, H. Hofling and Tibshinari, Pathwise coordinate optimization. Annals of Applied Statistics. I(2):302-332
【発明の開示】
【発明が解決しようとする課題】
【００１７】
この発明の目的は、各センサの周りの近接グラフ構造をデータから自動的に推測する機能を提供することによって、上記従来技術の問題点を解決することにある。
【課題を解決するための手段】
【００１８】
本発明において、入力として与えられるのは、下記のものである。
(1) テスト・データＤ = {x(t)|x(t)∈Ｒp, t = 1,2,...,N}
(2) 正規化定数ρ ここで、0 < ρ < 1であり、その値は、ユーザが適宜選択する。
(3) 参照用正常データＤ〜 = {x~(t)|x~(t)∈Ｒp, t = 1,2,...,N}の相関係数行列Ｓ〜
【００１９】
これによって出力されるのは、各センサの異常度e1, e2,...,epである。
【００２０】
本発明に従うアルゴリズムの概要は、次のとおりである。以下の処理はすべて、コンピュータのプログラムによって、自動的に実行されることに留意されたい。まず、corr()という関数表示で、相関係数行列を求める手続きをあらわす。すると、最初のステップは、下記のとおり、テスト・データＤの相関係数行列Ｓを求めることである。
１．Ｓ ← corr(Ｄ)
２．次のステップは、Ｄ、及びＤ〜に対する精度行列の計算である。精度行列とは基本的に、相関係数行列の逆行列であるが、特に本発明においては、L1ノルムの正規化項を付加した式を最適化するように解くことによって、精度行列が疎行列になるような解を得る。このような最適化の式を解くために好適なアルゴリズムの１つに、graphical lassoがある。それを、glasso()という関数表示であらわすと、
Λ ← glasso(ρ,Ｓ)、Λ〜 ← glasso(ρ,Ｓ〜)で、テスト・データＤの精度行列Λと、参照用正常データＤ〜の精度行列Λ〜を得る。
３，次のステップは、i = 1,2,...,pに対する特徴量行列の計算である。特徴量行列は、次のようにあらわされる。
【数１】
但し、li(A|B)というのは、negentropy(A,B)のことである。
４．次のステップは、下記の式によって、i = 1,2,...,pに対する異常度を計算することである。
ei = max{li(Λ|Ｓ) - li(Λ〜|Ｓ),li(Λ|Ｓ〜) - li(Λ〜|Ｓ〜)}
これは、相対エントロピー、またはカルバック・ライブラー距離として知られているものである。
【発明の効果】
【００２１】
以上のように、この発明によれば、精度行列が疎行列になるような処理の工夫を施して、好適にはgraphic lassoを用いてこの最適化の式を解き、各センサについて、カルバック・ライブラー距離が異常度をあらわすようにしたことによって、近傍性保存原理に基づく時系列データの異常検出手法に、各センサの周りの近接グラフ構造をデータから自動的に推測する機能が組み込まれ、さらには、局所的な統計モデルの直接比較に基づいて、理論的に首尾一貫した異常度が与えられる。
【００２２】
これによって、擬陽性の少ない、より適切なセンサ時系列データ異常判別技法が提供される。
【発明を実施するための最良の形態】
【００２３】
以下、図面に従って、本発明の実施例を説明する。これらの実施例は、本発明の好適な態様を説明するためのものであり、発明の範囲をここで示すものに限定する意図はないことを理解されたい。また、以下の図を通して、特に断わらない限り、同一符号は、同一の対象を指すものとする。
【００２４】
図２を参照すると、本発明の一実施例に係るシステム構成及び処理を実現するためのコンピュータ・ハードウェアのブロック図が示されている。図２において、システム・バス２０２には、ＣＰＵ２０４と、主記憶（ＲＡＭ）２０６と、ハードディスク・ドライブ（ＨＤＤ）２０８と、キーボード２１０と、マウス２１２と、ディスプレイ２１４が接続されている。ＣＰＵ１０４は、好適には、３２ビットまたは６４ビットのアーキテクチャに基づくものであり、例えば、インテル社のＰｅｎｔｉｕｍ（商標）４、インテル社のＣｏｒｅ（商標） ２ ＤＵＯ、ＡＭＤ社のＡｔｈｌｏｎ（商標）などを使用することができる。主記憶２０４は、好適には、５１２ＫＢ以上の容量、より好ましくは、１ＧＢ以上の容量をもつものである。
【００２５】
ハードディスク・ドライブ２０８には、個々に図示しないが、オペレーティング・システム及び本発明に係る処理プログラムなどが、予め格納されている。オペレーティング・システムは、Ｌｉｎｕｘ（商標）、マイクロソフト社のＷｉｎｄｏｗｓ Ｖｉｓｔａ、Ｗｉｎｄｏｗｓ ＸＰ（商標）、Ｗｉｎｄｏｗｓ（商標）２０００、アップルコンピュータのＭａｃ ＯＳ（商標）などの、ＣＰＵ２０４に適合する任意のものでよい。
【００２６】
ハードディスク・ドライブ２０８にはさらに、センサの基準用の正常の時系列データと、センサのテスト用時系列データとが個別のファイルとして格納されている。異なる状況で記録された複数のデータが、個別のファイルとして保存されていてよい。すると、ユーザは、テスト時に、保存されているファイルのうちの１つを選択することによって、所望の時系列データの異常検知を行うことができる。
【００２７】
同様に、基準用の正常の時系列データのファイルも、複数用意してハードディスク・ドライブ２０８に保存していてもよい。
【００２８】
キーボード２１０及びマウス２１２は、オペレーティング・システムが提供するグラフィック・ユーザ・インターフェースに従い、ディスプレイ２１４に表示されたアイコン、タスクバー、ウインドウなどのグラフィック・オブジェクトを操作するために使用される。キーボード２１０及びマウス２１２はまた、異常検知のためのテスト用センサ時系列データが格納されたファイルを指定するためにも使用される。
【００２９】
ディスプレイ２１４は、これには限定されないが、好適には、１０２４×７６８以上の解像度をもち、３２ビットtrue colorのＬＣＤモニタである。ディスプレイ２１４は、時系列データの波形や、異常検知の結果を表示したりするために使用される。
【００３０】
ハードディスク・ドライブ２０８にはさらに、本発明に係る異常検知を行うためのプログラムが格納されている。このプログラムは、Ｃ＋＋、Ｃ＃、Ｊａｖａ（商標）、Ｐｅｒｌ、Ｒｕｂｙなどの既存の任意のプログラム言語で書くことができる。オペレーティング・システムとして、Ｗｉｎｄｏｗｓ Ｖｉｓｔａ、Ｗｉｎｄｏｗｓ ＸＰ（商標）、Ｗｉｎｄｏｗｓ（商標）２０００などを使用する場合には、Ｗｉｎ３２ ＡＰＩの機能を利用して、ＧＵＩも含むアプリケーション・プログラムとして実装することができる。しかし、本発明に係る異常検知を行うためのプログラムは、ＣＵＩとしても実装することが可能である。
【００３１】
図３は、センサからのデータの間の関連性を示す図である。図示されているように、本発明の異常検知の仕組みでは、iを1からpまでの間の任意の整数としたとき、センサからのデータxiにつき、データxi以外のデータx1, ..., xj, ...,xi-1,xi+1,...,xpとの間の関係をより際立たせるようなアルゴリズムが適用される。すなわち、弱い関係同士のデータの関連性は無視してしまい、かなり関係の大きいデータ同士の関係性のみが残されるようになされる。このような関係性を行列表示すると、ほとんどの成分が0の疎行列になる。換言すると、このような関係性をあらわす行列である、精度行列が、疎行列であるような解を求めることが、本発明の１つの目標である。図３では、データxiに対して、xjとはある程度関係があり、xi+1とは強い関係があり、それ以外の点線であらわす関係は、無関係とみなすことを示している。
【００３２】
なお、精度行列については、C.M. Bishop, "Pattern Matching and Machine Learning". Springer Verlagの2.3.1節などを参照されたい。
【００３３】
図４は、本発明に係る処理の概要を示すフローチャートである。このフローチャートを実行するプログラムは、ハードディスク・ドライブ２０８に保存され、ユーザの操作により、オペレーティング・システムの働きで、主記憶２０６にロードされて実行可能となる。
【００３４】
ステップ４０２では、手続きcorr()によって、テスト・データＤの相関係数行列が計算される。
テスト・データＤの定義を再掲すると、 {x(t)|x(t)∈Ｒp, t = 1,2,...,N}である。ここで、時系列データの次元数pは、例えば50程度で、数百になることもある。Nは、時分割の数で、一般的には100から1000程度である。時分割の刻みの適切な間隔は、適用例によって異なるが、ここでは0.1秒である。但し、後の計算量に効いてくるのは、Nよりも寧ろ、次元数pである。というのは、pが直接、相関係数行列の行数及び列数になるからである。このデータは、ハードディスク・ドライブ２０８に保存されており、処理のため、本発明の処理プログラムにより、主記憶２０６に読み出される。
【００３５】
テスト・データＤの相関係数行列は、次のようにして計算される。先ず、x(t)のi番目の成分を、xi(t)とあらわすことにする。ここで、i∈{1,..,p}である。
【００３６】
そこで、処理プログラムは、xi(t)の、t = 1,2,...,Nについての平均miと標準偏差σiとを計算する。そして、処理プログラムは、平均miと標準偏差σiを用いて、次の式により、xi(t)を標準化する。
【数２】
【００３７】
こうしておいて、処理プログラムは、次の式により、相関係数行列Ｓを求める。
【数３】
【００３８】
なお、ここでは参照用正常データＤ〜 = {x~(t)|x~(t)∈Ｒp, t = 1,2,...,N}の相関係数行列Ｓ〜は、所与であると想定しているが、もしまだ計算されていないなら、ここで同様に計算する。
【００３９】
ステップ４０４で、処理プログラムは、下記のように、関数glasso()を用いて、テスト・データＤの精度行列Λと、参照用正常データＤ〜の精度行列Λ〜とを計算する。glassoとは、graphical lassoアルゴリズムを意味する。
Λ = glasso(ρ,Ｄ)
Λ〜 = glasso(ρ,Ｄ〜)
【００４０】
精度行列Λとは、定義によれば、相関係数行列のSの逆行列である。しかし、Ｓの逆行列を計算することは一般に不可能である。なぜならＳは一般に正則ではないからである。これは例えば、ある程度強い相関がある変数の対が存在する場合にそうなる。センサー系の時系列データではそのような状況は常に起こっていると考えなければならない。仮に、幸運にもＳが正則だったとしても、単純に逆行列を計算して精度行列Λを求めればよいかというと否である。それでは、精度行列Λが疎行列にならず、近傍を自動で選択するという本発明の狙いを満たさないためである。そのための工夫は後で説明するとして、疎行列である精度行列Λが求まったとしよう。すると、グラフィカルモデリングの理論によれば、精度行列Λの(i,j)成分(i≠j)が0であれば、xiとxjは、これら以外の全ての変数を与えたときに条件付独立である。すなわち、
【数４】
【００４１】
本発明は、この条件付き独立性を、近傍探索に利用する。すなわち、統計的に独立でないと判断されれば近傍とみなし、そうでなければ近傍でないとする。
例えば、x2について、精度行列Λ2,j ≠ 0となるjが4と9だけであれば、x2の近傍がx4とx9だけ、ということが分かる。従って、多くの行列要素がゼロであるような疎行列として精度行列Λを最適に決定することは、近傍の自動選択と等価である。
【００４２】
さて、前にも述べたように、センサ・データ自体が多変量正規分布に従うという保証はない。そこで寧ろ、近傍か否かを決める基準として、多変量ガウスモデルを援用することにする。そこで今は、平均ゼロを前提としているので、データに当てはめようとしているモデルは、平均ゼロ、精度行列Λのp次元ガウス分布である。それは、下記の式のようにあらわされる。
【数５】
ここで、detは行列式をあらわす。
【００４３】
ここで、特定の(i,j)、例えば、(i,j) = (1,2)に着目する。そして、x1とx2以外に適当な数字として、例えば、0をいれてみる。すると、指数関数の肩は、次のようになる。
【数６】
【００４４】
従って、Λ1,2 = 0であれば、x1とx2を結ぶ交差項はなく、すると、それらは、他の全てのデータを止めた時に独立となることが分かる。
【００４５】
そこで、精度行列Λの疎性（sparsity）の議論に移る。疎性を一旦離れて、精度行列をデータから素朴に最尤推定することを試みる。すると、テスト・データＤの対数尤度は、次のようになる。
【数７】
【００４６】
ここで、trは行列の対角和である。従って、尤度を最大にする精度行列を求めるためには、上記の式の{ }内を最大化すればよいことになる。
【００４７】
しかし生憎、このようにして求めた精度行列は、一般に疎にならない。そこで、本発明は、L1ノルムに基づく正規化項を付与して最尤推定すると疎行列が得られるという性質を利用する。このようにL1正規化項を付して疎な解を得る手法が統計学の分野でlasso （Least absolute shrinkage and selection operator）と呼ばれることが、ここで使うgraphical lassoの名前の所以である。
【００４８】
そこで、L1ノルムを付与した関数f(Λ;Ｓ,ρ)を次のように定義する。
【数８】
ここで、|Λi,j|は、行列Λの(i,j)要素の絶対値をあらわす。
【００４９】
すると、この関数f(Λ;Ｓ,ρ)を用いて、Λを求めることは、下記ような最適化問題を解くことに帰着される。
【数９】
【００５０】
この最適化問題を解くための好適な方法は、graphical lassoというアルゴリズムである。これの詳しい説明は、J. Friedman, T. Hastie, and R. Tibahirani, Sparse inverse covariance estimation with the graphical lasso. Biostatics, 9(3): 432-441, 2008 などにある。
【００５１】
正規化項
【数１０】
の存在によって、Λの多くの行列要素は、正定値性を失わない範囲で厳密にゼロになる。定数ρは、直感的に言うと、どの程度相関係数を無視するかの尺度であり、0 < ρ < 1の範囲でユーザが与える。graphical lassoの場合の好適な１つの値は、ρ = 0.3である。ρが1に近いと、得られる近接グラフは非常に小さくなり、ρ = 0だと、基本的に全てのセンサが結合した完全グラフが得られる。
【００５２】
なお、graphical lassoのアルゴリズムについては、後で、もう少し詳しく説明する。
【００５３】
このようにして、ステップ４０４で、テスト・データＤの精度行列Λと、参照用正常データＤ〜の精度行列Λ〜とが求められると、ステップ４０６での負のエントロピー（ネゲントロピー）の行列の計算と、ステップ４０８での異常度スコアの計算に移行する。ステップ４０６及びステップ４０８は、i = 1,2,...,pについて順次実行される。
【００５４】
ところで、テスト・データＤと、参照用正常データＤ〜との間に本質的な違いがなければ、それぞれから得られる精度行列ΛとΛ〜とが表す近接グラフの様子には大差はないはずである。すなわち、精度行列を求めるためのアルゴリズムに鑑みると、近接グラフ自体が、図１に示すような時系列データの見かけ上の相違とは関係なく求まるはずである。逆に言えば、ある変数xiに着目したとき、その周りの近接グラフに大きな変化がおきていれば、その変数に何か異常が生じている公算が高い。
【００５５】
さて、いまやテスト・データＤの精度行列Λと、参照用正常データＤ〜の精度行列Λ〜とが得られているので、
テスト・データＤの確率モデルp(x) = N(x|0,Λ-1)と、
テスト・データＤ〜の確率モデルp~(x) = N(x|0,Λ〜-1)とが得られる。そこで、これらを勘案して、変数xiの異常度を計算するに当たっては、xi以外の全ての変数x1,x2,...,xi-1, xi+1,...,xpを与えたときのxiの条件付き分布p(xi|rest)と、p~(xi|rest)の違いを評価するのが自然である。ここで、restというのは、xi以外の全ての変数x1,x2,...,xi-1, xi+1,...,xpのことである。そこで、記法の便宜上、xi以外の全ての変数x1,x2,...,xi-1, xi+1,...,xpをzi∈Ｒp-1とすると、ネゲントロピーの関数は、下記のような式となる。
【数１１】
ここで、< .. >Dという表示は、データＤでの経験分布をあらわす。
これらの式が、前述のガウス分布の式から見導かれることは、この分野の当業者に明らかであろう。これらの式に基づき、ステップ４０６で、個々のi=1,...,pにつき、負のエントロピーが計算される。
【００５６】
これらの結果に基づき、下記の式を用いて、処理プログラムは、i=1,2,...,pにつき、異常度eiの計算を行う。
【数１２】
この結果は、見て取れるように、相対エントロピー、すなわちカルバック・ライブラーの距離である。
そして、ei ≡ max(d1,d2)により、i番目のデータの異常度が求められる。なお、前述の２００７年１２月１８日出願の米国特許出願第１１／９５９０７３号明細書に記載の技法では、相違度行列Ｄの要素di,jが、- ln|Si,j|で定義されていたので、相関係数の符合反転を検知することができなかったが、本発明の技法には、そのような制約はない。
【００５７】
いくつかの実験で、本発明の技法は、ノイズが大きく、いくつかのセンサに断線が生じているような状況でも、高い精度で異常度が検知できることが分かった。これは、疎な精度行列としてテスト・データＤと、参照用正常データＤ〜を比較した場合、ノイズや断線の効果が、かなりの程度抑えられるからであろうと考察される。
【００５８】
＜grahical lassoの計算に関する補足説明＞
前述のように、grahical lassoにより解くべき問題は、下記の式を最大にする行列Λを求めることであった。
【数１３】
そこで、この式をΛで偏微分する。このとき、下記の公式に留意する。
【数１４】
【００５９】
すると、Λでの偏微分の結果の勾配の式は、下記のようになる。
【数１５】
ここで、sign(Λ)とは、符号関数sign()を、行列の要素に適用した結果を要素とする行列である。すなわち、行列sign(Λ)の(i,j)要素は、Λi,j≠0のときsign(Λi,j)であり、Λi,j=0であるなら、[-1,1]の何らかの値を返すように定義される。
【００６０】
上記勾配の式をゼロ行列と等値してΛを求めることになるが、Λはその性質から、対角要素は全て正であり、よって、sign(Λ)の対角要素は全て1となる。
【００６１】
これから、Σi,j = Si,j + ρと、対角要素は厳密に求まる。ここで、Σ≡Λ-1と定義した。
【００６２】
以上のように対角要素は厳密に求まるが、非対角要素は簡単には求まらない。そこで、ブロック勾配法という技法が使われる。ブロック勾配法については、"O. Banarjee, L.E.Ghaoui, and G. Natsoulis, Convex optimization techniques for fitting sparse gaussian graphical models. In in Proceedings of ICML, pages 89-96. Press, 2006."及び"J. Friedman, T. Hastie, H. Hofling and Tibshinari, Pathwise coordinate optimization. Annals of Applied Statistics. I(2):302-332"などにも記載がある。
ブロック勾配法では、次のように行列が分解される。
【数１６】
【００６３】
ここで、Ｌ及びＷは、(p-1)×(p-1)の行列であり、λとσはスカラーである。従って、wとlは、(p-1)次元のベクトル、ということになる。
【００６４】
ブロック勾配法は、ここでは、一つの列ベクトルであるwとlを未知、他のλ、σ、Ｌ、Ｗを既知として最適化問題を解くものである。
【００６５】
このように決めたwに対して、数１５がゼロになるという条件は、下記の式と等価である。βは、あるベクトル・パラメータ変数である。
【数１７】
但し、
【数１８】
である。また、sは、相関係数行列Ｓにおける、wに対応する位置のベクトルである。この式から見て取れるように、βは、L1正規化項付きの２次計画問題を解くことで求まる。よってβは、一般的な２次計画ソルバーを使っても解けるが、好適な解法があり、それについては、後述する。
【００６６】
数１５がゼロに等しいという条件の、Σ≡Λ-1におけるwに対応する位置のベクトルでの部分の式は、w - s - sign(l) = 0と書ける。
【００６７】
そこで、ΣΛを計算してみると、次のようになる。
【数１９】
この式の右上のブロックから得られる l = -λＷ-1wを使い、また、β≡Ｗ-1wとおけば、
Ｗβ - s +ρsign(β) = 0
数１７の偏微分を実行することによって、この式に等しくなることが確認できる。
【００６８】
上の説明では、第p変数についてブロック勾配法を適用するものであったが、より一般に次のようなステップであると言える。
(1) 最初に例えば、Σ = Ｓ + ρ１pと置く。
(2) i = 1から始めて、第i変数について、行と列を並べ替えて、i行目が第p行目、i列目が第p列目に来るようにする。
(3) lとλを最適に決めた後、行と列の配列を戻す。
(4) 元の行列と違いが小さければ終了。そうでなければ、iを１つ増やす。iがpを超えたなら、iを1にする。
(5) ステップ(2)に戻る。
【００６９】
このとき、lとλを最適に決めるステップであるが、数１９の右下部分から、wTl+σλ=1が得られ、これと、l = -λＷ-1wとから、下記の式が得られる。
【数２０】
【００７０】
ここで、結果的に使われる式には、Ｗ-1が陽に表れていないことに留意されたい。すなわち、相関係数行列がランク落ちしている場合、一般的にＷ-1を求めることができない。その意味で、逆行列を使わないこの計算方法は、ランク落ちしている相関係数行列にも対応できる点で、有利である。
【００７１】
さて、数１７に対応する最適化問題を書き下す。
【数２１】
この目的関数をg(β;Ｗ,ρ)と書き、βiで偏微分すると、
【数２２】
が成り立つ。これに基づく最適化条件を次の２つの場合に分けて考える。
【００７２】
(1) βi > 0のとき、数２２が0という条件は、βm (m≠i)を与えた下で、
【数２３】
と解ける。但し、
【数２４】
と置いた。数２３で、Ai - ρ < 0 であれば、もともとの条件βi > 0と矛盾する。このときは、
【数２５】
のように単調であることから、定義域の左端、すなわち、βi → 0が解となる。
【００７３】
(2) βi < 0のとき、数２２が0という条件は、
【数２６】
と、解ける。但し、Ai + ρ > 0であれば、もともとの条件βi < 0と矛盾する。このときは、
【数２７】
のように単調であることから、定義域の左端、すなわち、βi → 0が解となる。
【００７４】
以上、２つの場合を纏めると、数２１の最適化問題は、p-1個のiそれぞれについて、
【数２８】
という置換を繰り返せばよい。添え字iが最後の変数まで行ったら、また最初に戻って、収束するまで繰り返す。
【００７５】
なお、Σi,j = Si,j + ρという式から見てとれるように、Ｗi,iは必ず正であり、この解に特異性はない。また、上述したgraphical lassoの２つの記述法は、簡単な反復公式で記述でき、計算効率が高い。らに、前述のように、Σがランク落ちしていも、その逆行列を計算に使う必要がないので、安定して計算を実行することができる。
【００７６】
以上の例は、自動車のセンサから取得したデータを異常検知する場合であったが、これには限定されず、発電所や、プラントなど、測定データが多次元時系列データとあらわされる、任意の場合に適用可能である。
【図面の簡単な説明】
【００７７】
【図１】参照用正常時系列データと、テスト用時系列データの波形の例を示す図である。
【図２】本発明を実施するためのハードウェアのブロック図である。
【図３】変数間の近傍関係を示す図である。
【図４】本発明の処理のフローチャートを示す図である。
【特許請求の範囲】
【請求項１】
コンピュータにより、多次元時系列データの異常度を計算するためのシステムであって、
第１の多次元時系列データと、該第１の多次元時系列データと同次元の第２の多次元時系列データとを入力する手段と、
前記第１の多次元時系列データの共分散行列の逆行列である第１の精度行列と、前記第２の多次元時系列データの共分散行列の逆行列である第２の精度行列とを、各々の精度行列が疎行列になるように、計算する手段と、
前記第１の精度行列に基づく多変量確率分布モデルと、第２の精度行列に基づく多変量確率分布モデルによって、前記多次元時系列データの各次元毎の負のエントロピーを計算する手段と、
前記計算された負のエントロピーに基づく、相対エントロピー距離に基づき、各次元毎の異常度を計算する手段とを有する、
多次元時系列データの異常度計算システム。
【請求項２】
前記第１の第多次元時系列データが参照用正常データであり、前記第２の第多次元時系列データが、異常を検知すべきテスト・データである、請求項１のシステム。
【請求項３】
前記多変量確率分布モデルが、多変量ガウスモデルである、請求項１のシステム。
【請求項４】
前記精度行列が疎行列になるように計算する手段は、L1ノルムの正規化項を付加した式を最適化することにより計算を行う、請求項１のシステム。
【請求項５】
前記精度行列が疎行列になるように計算する手段は、graphical lassoのアルゴリズムを用いる、請求項４のシステム。
【請求項６】
前記多次元時系列データが、自動車の複数のセンサから取得されたデータであり、前記多次元の個々の要素が、個々のセンサに対応する、請求項１のシステム。
【請求項７】
コンピュータにより、多次元時系列データの異常度を計算するための方法であって、
第１の多次元時系列データと、該第１の多次元時系列データと同次元の第２の多次元時系列データとを入力するステップと、
前記第１の多次元時系列データの共分散行列の逆行列である第１の精度行列と、前記第２の多次元時系列データの共分散行列の逆行列である第２の精度行列とを、各々の精度行列が疎行列になるように、計算するステップと、
前記第１の精度行列に基づく多変量確率分布モデルと、第２の精度行列に基づく多変量確率分布モデルによって、前記多次元時系列データの各次元毎の負のエントロピーを計算するステップと、
前記計算された負のエントロピーに基づく、相対エントロピー距離に基づき、各次元毎の異常度を計算するステップとを有する、
多次元時系列データの異常度計算方法。
【請求項８】
前記第１の第多次元時系列データが参照用正常データであり、前記第２の第多次元時系列データが、異常を検知すべきテスト・データである、請求項７の方法。
【請求項９】
前記多変量確率分布モデルが、多変量ガウスモデルである、請求項７の方法。
【請求項１０】
前記精度行列が疎行列になるように計算するステップは、L1ノルムの正規化項を付加した式を最適化することにより計算を行う、請求項７の方法。
【請求項１１】
前記精度行列が疎行列になるように計算するステップは、graphical lassoのアルゴリズムを用いる、請求項１０のシステム。
【請求項１２】
前記多次元時系列データが、自動車の複数のセンサから取得されたデータであり、前記多次元の個々の要素が、個々のセンサに対応する、請求項７の方法。
【請求項１３】
コンピュータにより、多次元時系列データの異常度を計算するためのプログラムであって、
前記コンピュータに、
第１の多次元時系列データと、該第１の多次元時系列データと同次元の第２の多次元時系列データとを入力するステップと、
前記第１の多次元時系列データの共分散行列の逆行列である第１の精度行列と、前記第２の多次元時系列データの共分散行列の逆行列である第２の精度行列とを、各々の精度行列が疎行列になるように、計算するステップと、
前記第１の精度行列に基づく多変量確率分布モデルと、第２の精度行列に基づく多変量確率分布モデルによって、前記多次元時系列データの各次元毎の負のエントロピーを計算するステップと、
前記計算された負のエントロピーに基づく、相対エントロピー距離に基づき、各次元毎の異常度を計算するステップとを実行させる、
プログラム。
【請求項１４】
前記第１の第多次元時系列データが参照用正常データであり、前記第２の第多次元時系列データが、異常を検知すべきテスト・データである、請求項１３のプログラム。
【請求項１５】
前記多変量確率分布モデルが、多変量ガウスモデルである、請求項１３のプログラム。
【請求項１６】
前記精度行列が疎行列になるように計算するステップは、L1ノルムの正規化項を付加した式を最適化することにより計算を行う、請求項１３のプログラム。
【請求項１７】
前記精度行列が疎行列になるように計算するステップは、graphical lassoのアルゴリズムを用いる、請求項１６のプログラム。
【請求項１８】
前記多次元時系列データが、自動車の複数のセンサから取得されたデータであり、前記多次元の個々の要素が、個々のセンサに対応する、請求項１３のプログラム。
【請求項１】
コンピュータにより、多次元時系列データの異常度を計算するためのシステムであって、
第１の多次元時系列データと、該第１の多次元時系列データと同次元の第２の多次元時系列データとを入力する手段と、
前記第１の多次元時系列データの共分散行列の逆行列である第１の精度行列と、前記第２の多次元時系列データの共分散行列の逆行列である第２の精度行列とを、各々の精度行列が疎行列になるように、計算する手段と、
前記第１の精度行列に基づく多変量確率分布モデルと、第２の精度行列に基づく多変量確率分布モデルによって、前記多次元時系列データの各次元毎の負のエントロピーを計算する手段と、
前記計算された負のエントロピーに基づく、相対エントロピー距離に基づき、各次元毎の異常度を計算する手段とを有する、
多次元時系列データの異常度計算システム。
【請求項２】
前記第１の第多次元時系列データが参照用正常データであり、前記第２の第多次元時系列データが、異常を検知すべきテスト・データである、請求項１のシステム。
【請求項３】
前記多変量確率分布モデルが、多変量ガウスモデルである、請求項１のシステム。
【請求項４】
前記精度行列が疎行列になるように計算する手段は、L1ノルムの正規化項を付加した式を最適化することにより計算を行う、請求項１のシステム。
【請求項５】
前記精度行列が疎行列になるように計算する手段は、graphical lassoのアルゴリズムを用いる、請求項４のシステム。
【請求項６】
前記多次元時系列データが、自動車の複数のセンサから取得されたデータであり、前記多次元の個々の要素が、個々のセンサに対応する、請求項１のシステム。
【請求項７】
コンピュータにより、多次元時系列データの異常度を計算するための方法であって、
第１の多次元時系列データと、該第１の多次元時系列データと同次元の第２の多次元時系列データとを入力するステップと、
前記第１の多次元時系列データの共分散行列の逆行列である第１の精度行列と、前記第２の多次元時系列データの共分散行列の逆行列である第２の精度行列とを、各々の精度行列が疎行列になるように、計算するステップと、
前記第１の精度行列に基づく多変量確率分布モデルと、第２の精度行列に基づく多変量確率分布モデルによって、前記多次元時系列データの各次元毎の負のエントロピーを計算するステップと、
前記計算された負のエントロピーに基づく、相対エントロピー距離に基づき、各次元毎の異常度を計算するステップとを有する、
多次元時系列データの異常度計算方法。
【請求項８】
前記第１の第多次元時系列データが参照用正常データであり、前記第２の第多次元時系列データが、異常を検知すべきテスト・データである、請求項７の方法。
【請求項９】
前記多変量確率分布モデルが、多変量ガウスモデルである、請求項７の方法。
【請求項１０】
前記精度行列が疎行列になるように計算するステップは、L1ノルムの正規化項を付加した式を最適化することにより計算を行う、請求項７の方法。
【請求項１１】
前記精度行列が疎行列になるように計算するステップは、graphical lassoのアルゴリズムを用いる、請求項１０のシステム。
【請求項１２】
前記多次元時系列データが、自動車の複数のセンサから取得されたデータであり、前記多次元の個々の要素が、個々のセンサに対応する、請求項７の方法。
【請求項１３】
コンピュータにより、多次元時系列データの異常度を計算するためのプログラムであって、
前記コンピュータに、
第１の多次元時系列データと、該第１の多次元時系列データと同次元の第２の多次元時系列データとを入力するステップと、
前記第１の多次元時系列データの共分散行列の逆行列である第１の精度行列と、前記第２の多次元時系列データの共分散行列の逆行列である第２の精度行列とを、各々の精度行列が疎行列になるように、計算するステップと、
前記第１の精度行列に基づく多変量確率分布モデルと、第２の精度行列に基づく多変量確率分布モデルによって、前記多次元時系列データの各次元毎の負のエントロピーを計算するステップと、
前記計算された負のエントロピーに基づく、相対エントロピー距離に基づき、各次元毎の異常度を計算するステップとを実行させる、
プログラム。
【請求項１４】
前記第１の第多次元時系列データが参照用正常データであり、前記第２の第多次元時系列データが、異常を検知すべきテスト・データである、請求項１３のプログラム。
【請求項１５】
前記多変量確率分布モデルが、多変量ガウスモデルである、請求項１３のプログラム。
【請求項１６】
前記精度行列が疎行列になるように計算するステップは、L1ノルムの正規化項を付加した式を最適化することにより計算を行う、請求項１３のプログラム。
【請求項１７】
前記精度行列が疎行列になるように計算するステップは、graphical lassoのアルゴリズムを用いる、請求項１６のプログラム。
【請求項１８】
前記多次元時系列データが、自動車の複数のセンサから取得されたデータであり、前記多次元の個々の要素が、個々のセンサに対応する、請求項１３のプログラム。
【図１】
【図２】
【図３】
【図４】
【図２】
【図３】
【図４】
【公開番号】特開２０１０−７８４６７（Ｐ２０１０−７８４６７Ａ）
【公開日】平成２２年４月８日（２０１０．４．８）
【国際特許分類】
【出願番号】特願２００８−２４７３８０（Ｐ２００８−２４７３８０）
【出願日】平成２０年９月２６日（２００８．９．２６）
【出願人】（３９０００９５３１）インターナショナル・ビジネス・マシーンズ・コーポレーション (4,084)
【氏名又は名称原語表記】ＩＮＴＥＲＮＡＴＩＯＮＡＬ ＢＵＳＩＮＥＳＳ ＭＡＳＣＨＩＮＥＳ ＣＯＲＰＯＲＡＴＩＯＮ
【Ｆターム（参考）】
【公開日】平成２２年４月８日（２０１０．４．８）
【国際特許分類】
【出願日】平成２０年９月２６日（２００８．９．２６）
【出願人】（３９０００９５３１）インターナショナル・ビジネス・マシーンズ・コーポレーション (4,084)
【氏名又は名称原語表記】ＩＮＴＥＲＮＡＴＩＯＮＡＬ ＢＵＳＩＮＥＳＳ ＭＡＳＣＨＩＮＥＳ ＣＯＲＰＯＲＡＴＩＯＮ
【Ｆターム（参考）】
[ Back to top ]