JP4953239B2 - 観測対象の異常を検出する技術 - Google Patents

Description

本発明は、観測対象の異常を検出する技術に関する。特に、本発明は、観測対象から得られた観測値の時系列データに基づいて、観測対象の異常を検出する技術に関する。

近年、工場設備の稼働状況、自動車の走行状態、あるいは経済活動の指標など、様々な観測対象から様々な観測値を取得することができる。観測対象は複雑化しており、観測値がどのような値のときに観測対象がどのような状態にあるのかを明確に把握することは難しい場合がある。たとえば、観測値がある閾値を超えた場合に観測対象に異常が発生したと判断する方法も考えられるが、その閾値を適切に設定するのは容易でない。１つの観測対象から多数の観測値が並行して取得される場合には、そのそれぞれの性質について充分な予備知識を得ることは難しく、閾値の設定は一層困難になる。また、観測対象の状態をリアルタイムで把握することは、観測値の解析に要する時間も限定されるので、更に難しい。

M. Brand. Fast online SVD revisions for lightweight recommender systems. In Proc. SIAM International Conference on Data Mining, 2003. G. H. Golub and C. F. V. Loan. Matrix computations (3rd ed.). Johns Hopkins University Press, Baltimore, MD, 1996. T. Id´e, "Why does subsequence time-series clustering produce sine waves?", In Proceedings Of the 10th European Conference on Principles and Practice of Knowledge Discovery in Databases (PKDD), Sep 18-22, 2006; Lecture Notes in Artificial Intelligence 4213, Springer, pp.311-322. T. Id´e and K. Inoue, "Knowledge discovery from heterogeneous dynamic systems using change-point correlations", In Proceedings of 2005 SIAM International Conference on Data Mining (SDM 05), pages 571-575, 2005. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes. Cambridge University Press, 1989. S. Roweis. EM algorithms for PCA and SPCA. In M. I. Jordan, M. J. Kearns, and S. A. Solla, editors, Advances in Neural Information Processing Systems, volume 10. The MIT Press, 1998. Gilbert Strang , "Linear Algebra and Its Applications", Academic Press, 1976. William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling, "Numerical Recipes in C: The Art of Scientific Computing", 2nd ed., Cambridge University Press, 1992, Section 11.3.

このように観測値の性質について予備知識が無い場合に有効な方法として、観測値の性質が変化したことを検出する方法が提案されている。たとえば、ウェーブレット変換を用いる方法が知られているが、この方法においてもパラメータの調整にある程度の予備知識が必要となる場合がある。また、観測値が多次元正規分布に従うと仮定して、その分布関数のパラメータの変化を検出する方法も用いられているが、観測値が多次元正規分布に従うという仮定は不適切な場合もある。

これに対し、特異スペクトル変換方法を利用した変化点解析によれば、観測値の性質について予備知識が無くとも、観測値の変化を適切に捉えることができる（非特許文献４を参照。）。ただし、この方法においては、特異スペクトル解析をするべき行列の次元数が大きくなると、現実的な計算時間で完了しない場合がある。この行列の次元数は、観測対象となる期間の長さに応じて大きくなることから、充分に長い観測期間について観測して観測対象の変化を精度良く捉えようとすると、リアルタイム用途には適用しにくい。

そこで本発明は、上記の課題を解決することのできるシステム、方法およびプログラムを提供することを目的とする。この目的は特許請求の範囲における独立項に記載の特徴の組み合わせにより達成される。また従属項は本発明の更なる有利な具体例を規定する。

上記課題を解決するために、本発明の一側面においては、観測対象を観測して得られた観測値の時系列データに基づいて、観測対象の異常を検出するシステムであって、時系列データのうち、異常検出の基準となる複数の期間のそれぞれにおいて観測された複数の部分時系列データを、それぞれ列ベクトルとして配列して第１行列を生成する第１生成部と、時系列データのうち、異常検出の対象となる複数の期間のそれぞれにおいて観測された複数の部分時系列データを、それぞれ列ベクトルとして配列して第２行列を生成する第２生成部と、第２行列の特異ベクトルを算出する特異ベクトル算出部と、第１行列に第１行列の転置行列を右から乗じて積行列を算出する行列積算出部と、特異ベクトルを基底として含む積行列のクリロフ部分空間の正規直交基底の各基底を列ベクトルとして配列した行列、および当該行列の転置行列を用いて３重対角化された、積行列の３重対角化行列の対角要素および副対角要素を算出する要素算出部と、３重対角化行列の固有ベクトルを算出する固有ベクトル算出部と、それぞれの固有ベクトルと特異ベクトルとの間の内積の線形結合に基づいて、異常検出の基準となる期間に対する異常検出の対象となる期間における観測対象の変化の程度を算出して、観測対象の異常を示す指標として出力する変化度算出部とを備えるシステムを提供する。また、当該システムにより観測対象の異常を検出する方法、および、当該システムとして情報処理装置を機能させるプログラムを提供する。
なお、上記の発明の概要は、本発明の必要な特徴の全てを列挙したものではなく、これらの特徴群のサブコンビネーションもまた、発明となりうる。

以下、発明を実施するための最良の形態（以下、実施形態と称す）を通じて本発明を説明するが、以下の実施形態は特許請求の範囲にかかる発明を限定するものではなく、また実施形態の中で説明されている特徴の組み合わせの全てが発明の解決手段に必須であるとは限らない。

図１は、検出システム３０周辺の構成概要を示す。観測対象を観測して得られた観測値の時系列データはデータベース装置２０に格納される。観測対象は、たとえば製鉄工場の製鉄用設備であり、観測値は、たとえばその製鉄設備で順次製造される帯鋼の張力である。観測対象としては、その他走行中の自動車であってもよい。この場合の観測値は、自動車の速度やエンジンの回転数のような連続値のほか、シフト位置や各種機能のＯＮ／ＯＦＦなどの離散値であってもよい。データベース装置２０は、観測値をリアルタイムで取得して格納してもよい。即ち、データベース装置２０は、観測対象から新たに観測値が観測される毎に、その新たな観測値を、既に格納している時系列データの末尾に追加してもよい。

検出システム３０は、データベース装置２０に格納されているこの時系列データに基づいて、観測対象の異常を検出する。たとえば製鉄用設備の例では、帯鋼の張力の推移に何らかの変化が検出されると、観測対象である製鉄用設備に何らかの異常が発生したと判断する。このような異常判断において、観測対象の特徴について充分な予備知識が得られない場合には、観測値と所定の閾値との比較などの方法によって異常を高精度に検出することは難しい。これに対し、本実施形態に係る検出システム３０は、観測対象についての予備知識を必要とすることなく、観測対象に生じた異常を高速な計算手法により高精度に検出することを目的とする。これにより、工場用設備や自動車などの複雑な観測対象につき、その設計情報等の機密性の高い情報を必要とすることなく、リアルタイムな異常検出を実現することができる。

図２は、データベース装置２０に記憶される時系列データの具体例を示す。この図に示す時系列データは、自動車各部を観測した観測値の時系列データであり、ある自動車シミュレータにより生成されたものである。図２（ａ）から図２（ｅ）までのそれぞれは、互いに異なる部分を観測した観測値の時系列データを示す。それぞれの時系列データは、１００ｍｓ間隔で観測された５８０個の観測値を含む。また、これらの時系列データは、回転数または速度などの物理量であってもよいし、シフト位置などの設定値であってもよい。この図においては便宜上、それぞれの観測値を分散が１で平均が３（無単位）となるように正規化して示している。

データベース装置２０は、このような各時系列データについて、当該時系列データに含まれるそれぞれの観測値を、その観測時刻に対応付けて格納している。そして、データベース装置２０は、検出システム３０からの読み出し指示に応じて、格納している観測値を検出システム３０に対し出力する。
以下、このように多数取得される時系列データのうちある１つに基づいて、その時系列データを観測した観測対象の異常を検出する機能について説明する。

図３は、検出システム３０の機能構成を示す。検出システム３０は、第１生成部３００と、第２生成部３１０と、特異ベクトル算出部３２０と、行列積算出部３３０と、要素算出部３４０と、固有ベクトル算出部３５０と、変化度算出部３６０とを有する。まず、これらの各部材と、ハードウェア資源との関係を述べる。検出システム３０は、後述の情報処理装置５００により実現される。情報処理装置５００に読み込まれたプログラムは、後述のＣＰＵ１０００およびＲＡＭ１０２０を協働させて、情報処理装置５００を、第１生成部３００と、第２生成部３１０と、特異ベクトル算出部３２０と、行列積算出部３３０と、要素算出部３４０と、固有ベクトル算出部３５０と、変化度算出部３６０として機能させる。

第１生成部３００は、ある時点を基準として、その基準の時点以降の各時刻に対応する観測値をデータベース装置２０から検索して、時系列データとしてデータベース装置２０から読み出す。そして、第１生成部３００は、その時系列データのうち、異常検出の基準となる複数の期間のそれぞれにおいて観測された複数の部分時系列データを、それぞれ列ベクトルとして配列して第１行列を生成する。生成した第１行列は、行列積算出部３３０に対し出力される。この行列をＨ_１と呼ぶ。

第１生成部３００は、生成した第１行列全体をＲＡＭ１０２０に記憶したうえで、行列積算出部３３０からの要求に応じて読み出して、行列積算出部３３０に対し出力してもよい。これに代えて、第１生成部３００は、行番号および列番号を指定した要求を行列積算出部３３０から受けて、その行番号および列番号により指定される第１行列中の要素を生成し、その要素のみを行列積算出部３３０に対し出力してもよい。この場合には、第１生成部３００は第１行列全体をＲＡＭ１０２０に記憶する必要がないので、ＲＡＭ１０２０の必要容量を小さくすることができる。

第２生成部３１０は、第１生成部３００と同様にデータベース装置２０から時系列データを読み出す。そして、第２生成部３１０は、異常検出の対象となる複数の期間のそれぞれにおいて観測された複数の部分時系列データを、それぞれ列ベクトルとして配列した第２行列を生成する。生成した第２行列は特異ベクトル算出部３２０に対し出力される。この第２行列を行列Ｈ_２と呼ぶ。以下、図４を参照して、第１行列および第２行列を対比して示す。

図４は、時系列データから生成する第１行列および第２行列の一例を示す。図中の折れ線グラフは時系列データを示す。図中の矩形の領域は、部分時系列データを生成する対象となる期間を模式的に示している。即ち、矩形の領域のそれぞれが各部分時系列データを構成する観測値が観測された期間に対応する。第１行列（Ｈ_１）は、異常検出の基準となる複数の期間のそれぞれにおける部分時系列データを配列したものであり、好ましくは観測対象から観測された観測値の履歴に基づいて生成される。そして、この複数の期間とは、ある基準の時点に所定の差分期間を順次加えた複数の時点のそれぞれを開始時点とした、互いに長さがｗとして等しい複数の期間であり、その数はｎ個である。たとえば、時系列データに含まれる観測値が、１，２，３，４，５，…，ｋ，…の場合に、ｗ＝３としたときの部分時系列データは、（１，２，３），（２，３，４），…，（ｋ，ｋ＋１，ｋ＋２）などである。

また、基準の時点とは、リアルタイムで新たな観測値が観測される毎に、その観測値の観測時点から予め定められた期間遡った時点である。また、長さｗの単位は、ここでは、定期的に観測される観測値の数である。この場合の、所定の差分期間とは、たとえば、観測値の観測間隔１つ分である。これに代えて、ｗは、観測期間の時間の長さを示すものであり、差分期間も時間の長さとして予め定められてもよい。第２行列（Ｈ_２）についても同様に、この行列は、異常検出の対象となる複数の期間のそれぞれにおける部分時系列データを配列したものである。そして、この複数の期間とは、たとえば、上記の基準の時点よりも所定の期間（γ）後の時点を基準として、当該基準の時点に所定の差分期間を順次加えた複数の時点のそれぞれを開始時点とした、互いに長さがｗとして等しい複数の期間である。この期間の数もまたｎ個である。

本実施形態に係る検出システム３０は、第２行列Ｈ_２によって表される時系列データの断片を、第１行列Ｈ_１によって表される時系列データの断片と比較することによって、新たな観測値が観測された時点の近傍、たとえば現在時刻よりも少し遡ったある時刻ｔにおいて観測対象に異常が生じたことを検出する。このように、観測期間にｗというある程度の幅をもたせることと、観測期間を少しずつシフトした複数の部分系列を滑走窓式に取り出すことで、異常検出の信頼性を高められることが知られている。また、このように取り出した部分系列を列ベクトルとして配列した第１行列Ｈ_１および第２行列の対比については、特異ベクトルの比較が有効であることが知られており、このアイデアはＳｉｎｇｕｌａｒ Ｓｐｅｃｔｒｕｍ Ｔｒａｎｓｆｏｒｍａｔｉｏｎ（ＳＳＴ）と呼ばれている（上述の非特許文献４を参照。）。

たとえば、第１行列Ｈ_１の左特異ベクトルを、対応する特異値が大きい方から順にｋ個算出し、それらによって構成される超平面と、第２行列Ｈ_２の左特異ベクトルとがなす角度を求めることができれば、観測値の変化の程度を示す指標として利用することができるとも考えられる。しかしながら、この角度算出の前提として必要な左特異ベクトルの算出は、計算に時間がかかるため現実的な時間で終了しない場合がある。たとえば、特異ベクトルの算出に用いる特異値分解（ＳＶＤ：ｓｉｎｇｕｌａｒ ｖａｌｕｅ ｄｅｃｏｍｐｏｓｉｔｉｏｎ）において、Ｈｏｕｓｅｈｏｌｄｅｒ法を用いると、行列のサイズが大きい場合には膨大な計算時間を要して現実的な時間で終了しない。一方、行列のサイズを小さくすると異常検出の精度は低下する。疎行列の行列演算についての高速な解法が知られているが（非特許文献２を参照）、上記第１行列Ｈ_１は密行列でありそのまま適用することはできない。その他、行列演算の近似解法として直交反復法が知られているが（非特許文献２を参照）、その速度は充分なものではない。

これに対し、本実施形態に係る検出システム３０は、第１行列Ｈ_１に現れる観測値の特徴を充分に反映しつつも、その行数よりも充分に小さい次元数の部分空間に基づき、観測値の傾向の変化を捉える指標を算出する。これにより、第１行列Ｈ_１の特異ベクトルを直接に計算する必要が無いので、計算に要する時間を大幅に削減して処理を効率化できる。以下、具体的に説明する。
なお、以降の説明の前提として、特異ベクトルまたは固有ベクトルは、Ｌ２ノルムが１に規格化されていると仮定する。また、窓幅ｗを入力パラメターとし、これに対してｎとγを、それぞれｎ＝ｗおよびγ＝ｎ／２と選んでおく（結果はｎおよびγの選択に敏感ではないので、他の値でも構わない）。

図３に戻る。特異ベクトル算出部３２０は、第２行列Ｈ_２の左特異ベクトルのうち、最大の特異値に対応するものを算出する。算出されるこの左特異ベクトルをμとする。行列積算出部３３０は、第１行列Ｈ_１にこの第１行列の転置行列を右から乗じて積行列を算出する。この積行列をρとおく。即ちρ＝Ｈ_１Ｈ_１ ^Ｔである。この積行列ρの固有ベクトルは、もとの行列Ｈ_１の特異ベクトルと等しくなるから（非特許文献７参照）、以降この固有ベクトルを近似的に求めていく方針で説明をすすめる。次に、要素算出部３４０は、この積行列ρを３重対角化した３重対角化行列Ｔの対角要素α_ｓおよび副対角要素β_ｓを算出する。この３重対角化行列Ｔは、左特異ベクトルμを基底として含むこの積行列ρのクリロフ部分空間の正規直交基底の各基底を列ベクトルｑ_ｉとして配列した行列Ｑ、および、その行列の転置行列Ｑ^Ｔ用いて３重対角化されたものと一致する。クリロフ部分空間の次元数ｒを用いて、対角要素の各々はα_０からα_ｒ−１と表され、副対角要素のそれぞれはβ_０からβ_ｒ−１と表される。

固有ベクトル算出部３５０は、この３重対角化行列Ｔの固有ベクトルｘ^（ｉ）を算出する。好ましくは、固有ベクトル算出部３５０は、この固有ベクトルのうち、固有値の大きい方から順に、クリロフ部分空間の次元数ｒより少ないｋ個の固有ベクトルｘ^（１）からｘ^（ｋ）を算出する。そして、変化度算出部３６０は、それぞれの固有ベクトルと特異ベクトルμとの間の内積の線形結合に基づいて、異常検出の基準となる期間に対する異常検出の対象となる期間における観測対象の程度を示す指標値ｚを算出する。算出されるこの指標値は、第１行列Ｈ_１の特異ベクトルにより構成される超平面に対する、第２行列Ｈ_２の特異ベクトルμがなす角の大きさを示すものである。そして、変化度算出部３６０は、算出したその指標値ｚを、観測対象の異常を示す指標の指標値として出力する。この指標は、たとえば、製鉄用設備の例では、帯鋼の張力の変化を示す指標であり、これにより、製鉄用設備の異常を検出することができる。

次に、図５から図６を参照して、変化度を算出する処理につき更に詳細に説明するが。その説明に先立って、算出するべき指標値ｚの定義について述べる。指標ｚは、観測値の変化度を示すものである。変化度は、第１行列Ｈ_１および第２行列Ｈ_２の相違の程度を示すものであり、直感的には、第１行列Ｈ_１の特異ベクトルがなす超平面に対する、第２行列Ｈ_２の特異ベクトルのなす角度の大きさとして把握される。ここでは、このような角度と同様の性質を有しつつも、ベクトルの内積の線形結合によって表すことのできる指標を指標ｚとして導入する。具体的には、第１行列Ｈ_１の特異ベクトルがなす空間と、第２行列Ｈ_２の特異ベクトルがなす空間との間の距離の２乗として解釈できる量であり、以下の式（１）のように定義される。なお、空間同士の距離の概念については非特許文献２を参照されたい。

ただし、Ｋは内積若しくはそれに基づく量であり、以下の式（２）のように表される。

ｕ^（１）からｕ^（ｋ）は第１行列Ｈ_１の特異ベクトルを示す。この指標値を正確に算出するためには、特異ベクトルの算出処理が必要となり、既に述べたようにその計算の所要時間は極めて大きくなって実用的でない場合がある。以下、この式（１）および式（２）をそのまま計算することなく、この指標ｚの値を高精度かつ効率的に算出する処理の詳細を述べる。

図５は、順次観測される時系列データに基づき検出システム３０が変化度を算出する処理のフローチャートを示す。検出システム３０は、新たな観測値が観測される毎に、あるいは、順次観測値が観測される状況において定期的に、以下の処理を行う。まず、第１生成部３００は、現在の時点から予め定められた期間遡った時点を基準とした複数の期間について、第１行列Ｈ_１を生成する（Ｓ５００）。この第１行列Ｈ_１は、各観測値を要素とした要素数ｗの部分時系列データを、それぞれ列ベクトルとしてｎ本配列したものである。第１生成部３００は、この処理を観測値の取得に応じて順次行うことで、前回の第１行列の生成に用いた少なくとも１つの期間と比較して開始点がより後の複数の期間のそれぞれについて、次の第１行列を順次生成してゆく。

次に、第２生成部３１０は、第２行列Ｈ_２を生成する（Ｓ５１０）。この第２行列Ｈ_２は、第１行列Ｈ_１と比較して、観測開始の基準時点をγ遅らせたものである。本実施形態では上述のようにγ＝ｎ／２とおいている。第２生成部３１０は、観測対象から新たな観測値が観測される毎に上記の処理を行うことで、前回の第２行列の生成に用いた少なくとも１つの期間と比較して開始点がより後の複数の期間のそれぞれについて、次の第２行列を順次生成してゆく。

次に、特異ベクトル算出部３２０は、最大の特異値に対応する、第２行列Ｈ_２の左特異ベクトルを算出する（Ｓ５２０）。この算出処理は、たとえば直交反復法（非特許文献２を参照。）により実現される。直交反復法は、与えられた解の初期値に対し、反復処理によりその解を正解に徐々に近づけることにより、左特異ベクトルの近似解を精度良く算出するものである。したがって、解の初期値が正解に近い方が、計算の所要時間が小さくて済み、計算精度も高くなり易い。このため、好ましくは、特異ベクトル算出部３２０は、当該次の第２行列の左特異ベクトルを、前回の第２行列について算出した特異ベクトルに基づくベクトルを解の初期値とした直交反復法により算出する。前回の第２行列と今回の第２行列とでは、含んでいる観測値の大部分が共通であり、特異ベクトルも類似していると考えられるからである。更に好ましくは、特異ベクトル算出部３２０は、次の第２行列の特異ベクトルを、前回の第２行列について算出した特異ベクトルに、当該特異ベクトルとは方向が異なり当該特異ベクトルより短い他のベクトルを加えたベクトルを解の初期値とした直交反復法により算出してもよい。これにより、前回の特異ベクトルには含まれていなかった成分を解の初期値に僅かに含ませることができ、今回算出するべき特異ベクトルが前回の特異ベクトルには含まれていなかった成分を含む場合であっても、精度よい解を迅速に求めさせることができる。

次に、行列積算出部３３０は、積行列ρを算出する（Ｓ５３０）。この積行列ρを用いると、左特異ベクトルｕは、次の式（３）として示す変分方程式の解として求められる。次のステップＳ５４０の説明に先立ち、行列ρのクリロフ部分空間の次元数ｒが、積行列ρの行数ｗより充分に少なくても、この部分空間が行列ρの特異ベクトルの特徴を強く反映したものとなることを説明する。

R(u) はレイリー商と呼ばれる量である（詳しくは非特許文献７を参照。）。ここで次の問題を考える。「wより小さい任意の正整数s に対し、s +1次元の近似空間が、s 次元の近似空間よりもより「濃く」最大固有状態を含むような近似空間を構成せよ。」この問題の解を、1 次元空間を開始点として検討する。さしあたり任意の、L2ノルムが規格化された基底μを考える。上記の問題の要請を満たすためには、ｓｐａｎ｛μ｝に他の基底を付け加えてｓｐａｎ｛μ，Δ｝のような２次元空間を構成した場合に、これがＲの最急勾配方向を含む必要がある。但しｓｐａｎ｛ｘ，ｙ，ｚ…｝は、ベクトルｘ、ｙ、ｚにより張られる空間を示す。Ｒの勾配は以下の式（４）のように表される。

これにより、Δとしてρμを選択してｓｐａｎ｛μ，ρμ｝のような２次元空間を構成すれば、それがＲの最急勾配方向を含むことがわかる。以降、３次以上の次元の空間についても同様であるから、以下の式（５）によって表されるｒ次元空間は、Ｒの最大化の観点から最善のｒ次元部分空間であることがわかる。

言い換えると、第１行列Ｈ_１の列空間（w 次元）の部分空間としてのr次元空間のとり方は多くありえるが、その多くの選択肢のうちで、μをひとつの基底とする制約の下で、ρの最大固有状態を最も「濃く」含むものがＫ_ｒ（μ，ρ）である。最大固有状態ばかりではなく、第２、第３固有状態についても同様のことが言える。この部分空間は、数学では、（μとρから誘導される）クリロフ部分空間と呼ばれる。

次に、要素算出部３４０は、積行列ρを３重対角化した行列の対角要素および副対角要素を算出する（Ｓ５４０）。この３重対角化行列は、行列Ｑ、および、その行列Ｑの転置行列Ｑ^Ｔを用いて３重対角化されたものと一致する。行列Ｑは、左特異ベクトルμを先頭の基底として含む、積行列ρのクリロフ部分空間の正規直交基底の各基底を列ベクトルｑ_ｓとして配列した行列である。クリロフ部分空間の次元数をｒ（＜ｗ）とおくと、その正規直交基底はｒ個であるから、行列Ｑはｗ行ｒ列の非正方行列となり、以下の式（６）のように表される。

正規直交基底の定め方は様々であるから、行列Ｑは一般には一意に定まらない。しかしながら、第２行列Ｈ_２の左特異ベクトルμを１つの基底として含む正規直交基底を求めることとすれば、クリロフ行列［μ，ρμ，…，ρ^ｒ−１μ］のＱＲ分解により正規直交基底を求めることができ、ＱＲ分解の一意性より行列Ｑは一意に定まる。ここで、行列Ｑは非正方行列であるが、これが仮に正方行列であると仮定する。すると、行列ρについての固有方程式ρu＝λuは、クリロフ部分空間Ｋ_ｒ（μ，ρ）へのuの射影であるｘ＝Ｑ^Ｔｕを用いることで、近似的に以下の式（７）に変換される。

この式の左辺の係数行列Ｑ^ＴρＱは、ＱＲ分解の性質から、３重対角化行列となることが知られている（式（８）。非特許文献２参照）。Ｑが縦長の長方形行列であることから、その行数は積行列ρよりも少ない。

積行列ρは、原理的には、その積行列の正規直交基底の各基底を列ベクトルとして配列した正方行列Ｑを右から乗じ、かつ、その正方行列Ｑの転置行列Ｑ^Ｔを左から乗じることにより３重対角化されるが、実際には、下記の通り、行列Ｑをあらわに計算する必要はなく、Ｑとρの積を計算する必要もない。（なお、実際には行列Ｑは非正方行列である。このため、行列Ｔは、以上の導出過程を非正方行列Ｑについて近似的に適用して、即ち、正方行列に代えて非正方行列Ｑを右から乗じ、正方行列の転置に代えて非正方行列Ｑ^Ｔの転置を左から乗じることにより、近似的に３重対角化されるものである。ただし、これは行列Ｔとは如何なるものであるかを定義するものであり、実際には要素算出部３４０はＱとρの積を算出しない。）Ｑにより3重対角化された対角要素と副対角要素の算出は、積行列ρ、正規直交基底の１つであるμ、３重対角化行列の行数ｒを、ランチョスアルゴリズムに与えることによって実現される（ランチョスアルゴリズムについては非特許文献２を参照。）。

次に、固有ベクトル算出部３５０は、対角要素および副対角化要素を算出したこの３重対角化行列Ｔの固有ベクトルｘ^（ｉ）を算出する（Ｓ５５０）。好ましくは、固有ベクトル算出部３５０は、クリロフ部分空間の次元数ｒより少ないｋ個の固有ベクトルｘ^（１）からｘ^（ｋ）を算出する。3重対角行列の固有ベクトルは、きわめて効率よく計算できることが知られている。たとえば、QR反復アルゴリズムは、その反復過程において、3重対角行列の構造を崩すことなしに（つまり余計な非ゼロ成分を生成せずに）、固有値分解された行列に変形することができるため、3重対角行列に対しては非常に効率よく動作する。具体的な算出処理については、たとえば非特許文献５にｔｑｌｉルーチンとして記述されているので参照されたい。また、積行列ρはｗ行ｗ列であるのに対し、3重対角化行列Ｔはｒ行ｒ列と小さく、その上、固有値分解に特に有利な3重対角という構造をもっているから、その固有ベクトルは極めて高速に算出することが可能である。

次に、変化度算出部３６０は、算出したそれぞれの固有ベクトルｘ^（ｉ）と特異ベクトルμとに基づいて、異常検出の基準となる期間に対する異常検出の対象となる期間における観測対象の変化の程度を算出する（Ｓ５６０）。この詳細について説明する。いま、ここまでの説明によれば、固有ベクトル｛ｘ^（１），…，ｘ^（ｋ）｝は第１行列Ｈ_１の基底が張る空間における特異ベクトルではなく、そのクリロフ部分空間における特異ベクトルを表している。したがって、元の空間の特異ベクトルに基づいて、上述の式（１）に示す指標を算出しようとすれば、クリロフ部分空間の特異ベクトルを元の空間の特異ベクトルに変換する必要があるとも考えられる。この変換は、変換後の特異ベクトルを｛ｕ^（１），…，ｕ^（ｋ）｝とおくと、以下の式（９）のように表される。

ただし、ｘ_ｉ ^（α）は、ｒ次元のベクトルｘ^（α）の第ｉ番目の要素を表す。そして、このように求めた｛ｕ^（１），…，ｕ^（ｋ）｝に基づき、上述の式（１）における内積に基づくＫの値は、以下の式（１０）のように表される。

このように、固有ベクトルｘ^（ｉ）に基づき観測対象の変化の程度を示す指標を算出するためには、クリロフ部分空間の各基底ｑ_ｉを求めたうえで、式（１０）のような計算をすればよいとも考えられる。しかしながら、このような計算は不要である。これは、クリロフ部分空間の正規直交基底のうち１つを、第２行列Ｈ_２の左特異ベクトルμと設定したからである。これにより、基底ｑ_ｉ ^Ｔｑ_ｊ＝δ_ｉ，ｊとなり（δ_ｉ，ｊは、クロネッカーのデルタであり、ｉ＝ｊのときは１、ｉ≠ｊのときは０を採る。）、これに基づき式（１０）は式（１１）のように変形される。ただし、ｑ_１＝μである。

このように、式（１１）の右辺にはｑ_ｉが含まれていない。また、この式（１１）に基づき上述の式（１）を変形すると、以下の式（１２）のようになり、この式（１２）にもｑ_ｉが含まれていない。

従って、変化度算出部３６０は、この式（１２）に対し、固有ベクトル｛ｘ^（１），…，ｘ^（ｋ）｝をそれぞれ代入することにより、何らその他の変換処理を必要とすることなく、変化度を示す指標値ｚを算出することができる。この際、変換が必要であるとすれば必要な各基底ｑ_ｉや積行列ρは不要である。

図６は、Ｓ５４０における処理の具体例を示す。この処理は、積行列ρ、クリロフ部分空間の次元数ｒ、および、クリロフ部分空間の正規直交基底の１つであるμを入力として、積行列ρを３重対角化した行列Ｔの対角要素α_ｓおよび副対角要素β_ｓ（０≦ｓ≦ｒ−１）を求めることを目的とする。なお、計算の途中で一時変数ｓ、γ_０からγ_ｒ−１とｑ_０からｑ_ｒ−１を用いる。

まず、要素算出部３４０は、γ_０、β_０、ｑ_０およびｓを、それぞれ、μ、１、０および０により初期化する（Ｓ６００）。具体的には、たとえば、要素算出部３４０は、メモリ中に変数γ_０の値を格納する領域を確保して、その領域にベクトルμの値を格納する。以降の処理において、変数γは、添え字として変数ｓを伴って参照され、その他の変数を添え字として伴わないので、要素算出部３４０は、確保したこの領域を、変数γ_０からγ_ｒ−１のそれぞれの値を格納するための共通の領域として用いる。また、要素算出部３４０は、メモリ中に変数β_０の値を格納する領域を確保して、その領域に定数１を格納する。変数βについては添え字として１からｒ−１のそれぞれの値を伴って最終結果として出力されるものであるから、要素算出部３４０は、変数β_１からβ_ｒ−１のそれぞれの値を格納するための領域を確保する。

また、要素算出部３４０は、変数ｑ_０の値を格納する領域をメモリ中に確保して、その領域に定数０を格納する。変数ｑについては、添え字としてｓ＋１、ｓ、ｓ−１を伴って参照される。ただし、Ｓ６１０においてｑ_ｓ＋１に代入された値は、Ｓ６２０においてｓがインクリメントされた後の処理においてのみｑ_ｓとして参照される。このため、実質的に変数ｑは添え字としてｓまたはｓ−１のみを伴って参照される。したがって、要素算出部３４０は、変数ｑ０から変数ｑ_ｒ−１のそれぞれについてその値を格納する領域を確保する必要はなく、変数ｑ_ｓおよび変数ｑ_ｓ−１の値を格納する２つの領域をメモリに確保すればよい。Ｓ６００においては、たとえば、要素算出部３４０は、これら２つの領域を確保したうえで、変数ｑ_ｓの値を格納するための領域に定数０を格納する。また、要素算出部３４０は、メモリ中に、変数ｓの値を格納する領域を確保して、その領域の定数０を格納する。

次に、要素算出部３４０は、ベクトル変数であるγ_ｓのそれぞれの要素を、β_ｓで割ったベクトルを生成して、ベクトルｑ_ｓ＋１とする（Ｓ６１０）。具体的には、たとえば、要素算出部３４０は、変数γ_ｓおよびβ_ｓの値をメモリから読み出して、ＣＰＵ１０００の動作により演算処理を行ってベクトルｑ_ｓ＋１を生成し、変数ｑ_ｓを格納するために確保した領域に格納する。そして、要素算出部３４０は、変数ｓをインクリメントする（Ｓ６２０）。また、これに伴って、要素算出部３４０は、変数ｑ_ｓのベクトル値によって変数ｑ_ｓ−１のベクトル値を更新する。次に、要素算出部３４０は、変数ｓが次元数ｒと等しいか判断する（Ｓ６３０）。等しければ（Ｓ６３０：ＹＥＳ）、本図の処理を終了する。

等しくなければ（Ｓ６３０：ＮＯ）、要素算出部３４０は、ベクトルｑ_ｓの転置に積行列ρを右から乗じ、その結果に対しベクトルｑ_ｓをさらに右から乗じて、その結果をα_ｓとする（Ｓ６４０）。このとき、要素算出部３４０は、α_１からα_ｒ−１のそれぞれの値を格納するための領域が未確保であればそれらをメモリに確保したうえで、現在のｓの値に対応する領域のその計算結果を格納する。

次に、要素算出部３４０は、積行列ρにベクトルｑ_ｓを右から乗じたベクトルから、ベクトルｑ_ｓのそれぞれの要素にα_ｓを乗じたベクトル、および、ベクトルｑ_ｓ−１の各要素にβ_ｓ−１を乗じたベクトルを減じる演算を行い、その結果をベクトルγ_ｓとする（Ｓ６５０）。ベクトルｑ_ｓはベクトルｑ_ｓを記録するために確保した領域からベクトル値を読み出すことによって取得される。また、ベクトルｑ_ｓ−１はベクトルｑ_ｓ−１を記録するために確保した領域からベクトル値を読み出すことによって取得される。また、α_ｓおよびβ_ｓ−１についても同様である。算出されたγ_ｓは、γ_ｓを記録するための領域に対し、前回に算出したγ_ｓの値を上書きすることによって格納される。そして、要素算出部３４０は、算出したベクトルγ_ｓの転置とγ_ｓとの積を求めてその平方根を算出し、その結果をβ_ｓとする（Ｓ６６０）。要素算出部３４０は、Ｓ６１０に戻って処理を繰り返す。

以上、図６を参照して説明したように、要素算出部３４０によれば、対角要素および副対角要素の計算を、次元数ｒに比例した回数の繰返し処理により実現できる。これにより、反復解法により近似解の精度を徐々に改善して収束判定する方法と比較して、精度の良い解を安定的な時間で計算できる。また、繰返しの過程で変数γおよびｑについては、添え字としてｓかｓ−１のみを伴って参照されるので、変数の値を格納するために必要なメモリの容量も極めて小さくすることができる。

図７は、変化度算出の所要時間を示す。図７を参照して、本実施形態に係る検出システム３０の性能について説明する。異なる５つの手法のそれぞれによって、前述の自動車シミュレータで生成した時系列データについて、その変化度を算出し、その所要時間を算出した。Ｈｏｕｓｅｈｏｌｄｅｒ法（非特許文献５を参照。）とＱＲ反復法の組合せにより、第１行列Ｈ_１の特異値を直接に計算した手法の所要時間をシンボル「Ｈｏｕｓｅ」により示す。なお、所要時間はそれぞれ対数を取って所定の値を１として正規化して示している。

また、「Ｈｏｕｓｅ」の手法のうちＨｏｕｓｅｈｏｌｄｅｒ法を適用した部分を直交反復法に置換えた手法の所要時間を、シンボル「ＦＯ」により示す。観測期間の長さｗの値が大きくなると、「Ｈｏｕｓｅ」の手法よりも所要時間が短くなる。なお、他の手法も含めて、密行列に対してある程度の精度の解を求める手法の所要時間は、これらの手法と同等であることが知られている（たとえば非特許文献１を参照。）。また、本実施形態に係る検出システム３０による計算の所要時間をシンボル「ＰＫ」により示す。本実施形態との比較のため、疎行列の演算について優れた近似手法として知られるランチョスアルゴリズムに、ＱＬ反復法（非特許文献８を参照。）を組み合わせた手法によって、第１行列Ｈ_１から特異ベクトルを直接求めた手法を、シンボル「ＬＫ」により示す。また、やはり疎行列の演算手法として知られるＥＭ−ＰＣＡ法（非特許文献６を参照）によって、第１行列Ｈ_１から特異ベクトルを直接求めた手法を、シンボル「ＥＭＫ」により示す。但し、これらの手法は、本実施形態で取り扱う行列Ｈ_１のような密行列に適用しても充分な精度が得られない。これ似たいし、本実施形態に係る検出システム３０によれば、第１行列Ｈ_１の特異ベクトルを直接求めるという発想の範囲内の手法では達成できなかった高速化を図ることができる。

図８は、算出される変化度の具体例を時系列データと対比して示す。図８（ａ）に、観測値の時系列データを示す。この時系列データには、観測値の周期的な変化に加えて人工的なノイズを与えている。図８（ｂ）に、本実施形態に係る検出システム３０による変化度の算出結果を示す。また、図８（ｃ）に、上記「Ｈｏｕｓｅ」として説明した手法による変化度の算出結果を示す。ただし、上記で説明した変数ｋおよびｒの値はそれぞれ３および５に設定している。（なお、非特許文献３に記載されているように、固有ベクトルの数を表す変数ｋの値は、観測期間中の観測値の数ｗに関わらず、３程度の小さい整数でよいことが知られている。）図８（ｂ）および（ｃ）を対比しても違いは殆ど見当たらない。このように、本実施形態に係る検出システム３０によれば、変化度算出の精度についても充分に高いことが確かめられた。

図９は、検出システム３０として機能する情報処理装置５００のハードウェア構成を示す。情報処理装置５００は、ホストコントローラ１０８２により相互に接続されるＣＰＵ１０００、ＲＡＭ１０２０、及びグラフィックコントローラ１０７５を有するＣＰＵ周辺部と、入出力コントローラ１０８４によりホストコントローラ１０８２に接続される通信インターフェイス１０３０、ハードディスクドライブ１０４０、及びＣＤ−ＲＯＭドライブ１０６０を有する入出力部と、入出力コントローラ１０８４に接続されるＲＯＭ１０１０、フレキシブルディスクドライブ１０５０、及び入出力チップ１０７０を有するレガシー入出力部とを備える。

ホストコントローラ１０８２は、ＲＡＭ１０２０と、高い転送レートでＲＡＭ１０２０をアクセスするＣＰＵ１０００及びグラフィックコントローラ１０７５とを接続する。ＣＰＵ１０００は、ＲＯＭ１０１０及びＲＡＭ１０２０に格納されたプログラムに基づいて動作し、各部の制御を行う。グラフィックコントローラ１０７５は、ＣＰＵ１０００等がＲＡＭ１０２０内に設けたフレームバッファ上に生成する画像データを取得し、表示装置１０８０上に表示させる。これに代えて、グラフィックコントローラ１０７５は、ＣＰＵ１０００等が生成する画像データを格納するフレームバッファを、内部に含んでもよい。

入出力コントローラ１０８４は、ホストコントローラ１０８２と、比較的高速な入出力装置である通信インターフェイス１０３０、ハードディスクドライブ１０４０、及びＣＤ−ＲＯＭドライブ１０６０を接続する。通信インターフェイス１０３０は、ネットワークを介して外部の装置と通信する。ハードディスクドライブ１０４０は、情報処理装置５００が使用するプログラム及びデータを格納する。ＣＤ−ＲＯＭドライブ１０６０は、ＣＤ−ＲＯＭ１０９５からプログラム又はデータを読み取り、ＲＡＭ１０２０又はハードディスクドライブ１０４０に提供する。

また、入出力コントローラ１０８４には、ＲＯＭ１０１０と、フレキシブルディスクドライブ１０５０や入出力チップ１０７０等の比較的低速な入出力装置とが接続される。ＲＯＭ１０１０は、情報処理装置５００の起動時にＣＰＵ１０００が実行するブートプログラムや、情報処理装置５００のハードウェアに依存するプログラム等を格納する。フレキシブルディスクドライブ１０５０は、フレキシブルディスク１０９０からプログラム又はデータを読み取り、入出力チップ１０７０を介してＲＡＭ１０２０またはハードディスクドライブ１０４０に提供する。入出力チップ１０７０は、フレキシブルディスク１０９０や、例えばパラレルポート、シリアルポート、キーボードポート、マウスポート等を介して各種の入出力装置を接続する。

情報処理装置５００に提供されるプログラムは、フレキシブルディスク１０９０、ＣＤ−ＲＯＭ１０９５、又はＩＣカード等の記録媒体に格納されて利用者によって提供される。プログラムは、入出力チップ１０７０及び/又は入出力コントローラ１０８４を介して、記録媒体から読み出され情報処理装置５００にインストールされて実行される。プログラムが情報処理装置５００等に働きかけて行わせる動作は、図１から図８において説明した検出システム３０における動作と同一であるから、説明を省略する。

以上に示したプログラムは、外部の記憶媒体に格納されてもよい。記憶媒体としては、フレキシブルディスク１０９０、ＣＤ−ＲＯＭ１０９５の他に、ＤＶＤやＰＤ等の光学記録媒体、ＭＤ等の光磁気記録媒体、テープ媒体、ＩＣカード等の半導体メモリ等を用いることができる。また、専用通信ネットワークやインターネットに接続されたサーバシステムに設けたハードディスク又はＲＡＭ等の記憶装置を記録媒体として使用し、ネットワークを介してプログラムを情報処理装置５００に提供してもよい。

以上、本発明を実施の形態を用いて説明したが、本発明の技術的範囲は上記実施の形態に記載の範囲には限定されない。上記実施の形態に、多様な変更または改良を加えることが可能であることが当業者に明らかである。その様な変更または改良を加えた形態も本発明の技術的範囲に含まれ得ることが、特許請求の範囲の記載から明らかである。

図１は、検出システム３０周辺の構成概要を示す。 図２は、データベース装置２０に記憶される時系列データの具体例を示す。 図３は、検出システム３０の機能構成を示す。 図４は、時系列データから生成する第１行列および第２行列の一例を示す。 図５は、順次観測される時系列データに基づき検出システム３０が変化度を算出する処理のフローチャートを示す。 図６は、Ｓ５４０における処理の具体例を示す。 図７は、変化度算出の所要時間を示す。 図８は、算出される変化度の具体例を時系列データと対比して示す。 図９は、検出システム３０として機能する情報処理装置５００のハードウェア構成を示す。

符号の説明

２０ データベース装置
３０ 検出システム
３００ 第１生成部
３１０ 第２生成部
３２０ 特異ベクトル算出部
３３０ 行列積算出部
３４０ 要素算出部
３５０ 固有ベクトル算出部
３６０ 変化度算出部
５００ 情報処理装置

Claims (9)

1. 観測対象を観測して得られた観測値の時系列データに基づいて、前記観測対象の異常を検出するシステムであって、
   前記時系列データのうち、異常検出の基準となる複数の期間のそれぞれにおいて観測された複数の部分時系列データを、それぞれ列ベクトルとして配列して第１行列を生成する第１生成部と、
   前記時系列データのうち、異常検出の対象となる複数の期間のそれぞれにおいて観測された複数の部分時系列データを、それぞれ列ベクトルとして配列して第２行列を生成する第２生成部と、
   前記第２行列の特異ベクトルを算出する特異ベクトル算出部と、
   前記第１行列に前記第１行列の転置行列を右から乗じて積行列を算出する行列積算出部と、
   前記特異ベクトルを基底として含む前記積行列のクリロフ部分空間の正規直交基底の各基底を列ベクトルとして配列した行列、および当該行列の転置行列を用いて３重対角化された、前記積行列の３重対角化行列の対角要素および副対角要素を算出する要素算出部と、
   前記３重対角化行列の固有ベクトルを算出する固有ベクトル算出部と、
   それぞれの前記固有ベクトルと前記特異ベクトルとの間の内積の線形結合に基づいて、異常検出の基準となる期間に対する異常検出の対象となる期間における前記観測対象の変化の程度を算出して、前記観測対象の異常を示す指標として出力する変化度算出部と
   を備えるシステム。
2. 前記要素算出部は、前記積行列の行数よりも少ない次元数のクリロフ部分空間の正規直交基底を用いて３重対角化された、前記積行列よりも少ない行数の前記３重対角化行列の対角要素および副対角化要素を算出する
   請求項１に記載のシステム。
3. 前記積行列は、前記積行列の正規直交基底の各基底を列ベクトルとして配列した正方行列を右から乗じ、かつ、当該正方行列の転置行列を左から乗じることにより３重対角化されるものであり、
   前記要素算出部は、前記積行列に対し、前記正方行列に代えて、前記クリロフ部分空間の正規直交基底の各基底を列ベクトルとして配列した非正方行列を右から乗じ、前記正方行列の転置行列に代えて、当該非正方行列の転置行列を左から乗じることにより、近似的に３重対角化された３重対角化行列について、対角要素および副対角要素を算出する
   請求項２に記載のシステム。
4. 前記第２生成部は、前記観測対象から新たな観測値が観測される毎に、前回に前記第２行列の生成に用いた少なくとも１つの前記期間と比較して開始点がより後の複数の期間のそれぞれについて、次の第２行列を生成し、
   前記特異ベクトル算出部は、前記次の第２行列の特異ベクトルを、前回の前記第２行列について算出した特異ベクトルに基づくベクトルを解の初期値とした直交反復法により算出する
   請求項１に記載のシステム。
5. 前記特異ベクトル算出部は、前記次の第２行列の特異ベクトルを、前回の前記第２行列について算出した特異ベクトルに、当該特異ベクトルとは方向が異なり当該特異ベクトルより短い他のベクトルを加えたベクトルを解の初期値とした直交反復法により算出する
   請求項４に記載のシステム。
6. 前記第１生成部は、前記観測対象から新たな観測値が観測される毎に、当該新たな観測値の観測時点から予め定められた期間遡った時点を基準として、当該基準の時点に所定の差分期間を順次加えた複数の時点のそれぞれを開始時点とした、互いに長さの等しい複数の期間を、異常検出の基準となる前記複数の期間として前記第１行列を生成し、
   前記第２生成部は、前記観測対象から新たな観測値が観測される毎に、当該基準の時点よりも後の時点を基準として、当該基準の時点に所定の差分期間を順次加えた複数の時点のそれぞれを開始時点とした、互いに長さの等しい複数の期間を、異常検出の対象となる前記複数の期間として前記第２行列を生成し、
   前記変化度算出部は、前記観測対象から新たな観測値が観測される毎に、前記指標を算出して、前記新たな観測値が観測された時点の近傍における前記観測対象の異常を示す指標として出力する
   請求項１に記載のシステム。
7. 前記第１生成部は、製鉄用装置により順次製造される帯鋼の張力の観測値を時系列データとして前記第１行列を生成し、
   前記第２生成部は、前記製鉄用装置により順次製造される帯鋼の張力の観測値を時系列データとして前記第２行列を生成し、
   前記変化度算出部は、帯鋼の張力の変化を示す指標を算出して、前記製鉄用装置の異常を示す指標として出力する
   請求項１に記載のシステム。
8. 観測対象を観測して得られた観測値の時系列データに基づいて、前記観測対象の異常を検出する方法であって、
   前記時系列データのうち、異常検出の基準となる複数の期間のそれぞれにおいて観測された複数の部分時系列データを、それぞれ列ベクトルとして配列して第１行列を生成することと、
   前記時系列データのうち、異常検出の対象となる複数の期間のそれぞれにおいて観測された複数の部分時系列データを、それぞれ列ベクトルとして配列して第２行列を生成することと、
   前記第２行列の特異ベクトルを算出することと、
   前記第１行列に前記第１行列の転置行列を右から乗じて積行列を算出することと、
   前記特異ベクトルを基底として含む前記積行列のクリロフ部分空間の正規直交基底の各基底を列ベクトルとして配列した行列、および当該行列の転置行列を用いて３重対角化された、前記積行列の３重対角化行列の対角要素および副対角要素を算出することと、
   前記３重対角化行列の固有ベクトルを算出することと、
   それぞれの前記固有ベクトルと前記特異ベクトルとの間の内積の線形結合に基づいて、異常検出の基準となる期間に対する異常検出の対象となる期間における前記観測対象の変化の程度を算出して、前記観測対象の異常を示す指標として出力することと
   を備えるシステム。
9. 観測対象を観測して得られた観測値の時系列データに基づいて、前記観測対象の異常を検出するシステムとして、情報処理装置を機能させるプログラムであって、
   前記システムを、
   前記時系列データのうち、異常検出の基準となる複数の期間のそれぞれにおいて観測された複数の部分時系列データを、それぞれ列ベクトルとして配列して第１行列を生成する第１生成部と、
   前記時系列データのうち、異常検出の対象となる複数の期間のそれぞれにおいて観測された複数の部分時系列データを、それぞれ列ベクトルとして配列して第２行列を生成する第２生成部と、
   前記第２行列の特異ベクトルを算出する特異ベクトル算出部と、
   前記第１行列に前記第１行列の転置行列を右から乗じて積行列を算出する行列積算出部と、
   前記特異ベクトルを基底として含む前記積行列のクリロフ部分空間の正規直交基底の各基底を列ベクトルとして配列した行列、および当該行列の転置行列を用いて３重対角化された、前記積行列の３重対角化行列の対角要素および副対角要素を算出する要素算出部と、
   前記３重対角化行列の固有ベクトルを算出する固有ベクトル算出部と、
   それぞれの前記固有ベクトルと前記特異ベクトルとの間の内積の線形結合に基づいて、異常検出の基準となる期間に対する異常検出の対象となる期間における前記観測対象の変化の程度を算出して、前記観測対象の異常を示す指標として出力する変化度算出部と
   して機能させるプログラム。

Images