JPWO2020145215A1

JPWO2020145215A1 - 情報処理装置、情報処理方法及びプログラム

Info

Publication number: JPWO2020145215A1
Application number: JP2020565726A
Authority: JP
Inventors: 中川　淳一; 淳一中川; 秀樹南; 啓之上島; 章典竹下; 亮輔栗木
Original assignee: Nippon Steel Corp
Current assignee: Nippon Steel Corp
Priority date: 2019-01-09
Filing date: 2019-12-27
Publication date: 2021-09-30
Anticipated expiration: 2039-12-27
Also published as: WO2020145215A1; JP7036233B2

Abstract

物体の周期的な運動に係る計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定し、決定された修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、予め定められた状態に対応する修正係数のパターンを示すデータを格納する非負値の行列である教師基底行列と、データ行列における教師基底の重みを示す非負値の行列である教師重み行列と、データ行列から抽出される基底を格納する非負値の行列である抽出基底行列と、データ行列における抽出基底の重みを示す非負値の行列である抽出重み行列と、に分解し、教師重み行列と、抽出重み行列と、に基づいて、物体が予め定められた状態であるか否かを診断する。

Description

本発明は、情報処理装置、情報処理方法及びプログラムに関する。本願は、２０１９年１月９日に日本に出願された特願２０１９−２０５５号に基づき優先権を主張し、それらの内容を全てここに援用する。

周期的な運動を行う物体の診断を行う技術として、運動を行っている物体から得られる信号を計測し、その信号の特徴量を算出し、その特徴量から物体の診断を行う方法がある。特許文献１には、信号の周波数特性を自己回帰モデルで同定することで、より高精度なスペクトル特性を取得する技術が開示されている。特許文献２には、以下の技術が開示されている。まず、物体の周期的な運動から計測された計測データから、自己相関行列を生成する。自己相関行列を特異値分解することにより得られた固有値のうち、最大のものから設定された数を用いて、計測データを近似する修正された自己回帰モデルの係数を決定する。決定した係数から、修正された自己回帰モデルの周波数特性を求める。求めた周波数特性に基づいて、軸受の異常を診断する。

特開２００３−１１６８１１号公報特開２０１８−１６３１３５号公報

周期的な運動を行う物体が予め定められた状態であるか否かをより精度よく診断したいという要望がある。特許文献１、２等の従来技術では、この要望を満たすことができなかった。
本発明は、以上のような問題点に鑑みてなされたものであり、周期的な運動を行う物体が予め定められた状態であるか否かをより精度よく診断することを目的とする。

本発明の情報処理装置は、周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定手段と、前記決定手段により決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、抽出基底行列と、抽出重み行列と、に分解する分解手段と、前記分解手段により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断手段と、を有し、前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、前記決定手段は、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σ_ｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕ_ｓと、から導出される行列Ｕ_ｓΣ_ｓＵ_ｓ ^Ｔであり、前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、前記教師基底行列は、教師基底を格納する非負値の行列であり、前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、前記抽出基底行列は、抽出基底を格納する非負値の行列であり、前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、前記分解手段は、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解する。

本発明の情報処理方法は、情報処理装置が実行する情報処理方法であって、周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定ステップと、前記決定ステップで決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、抽出基底行列と、抽出重み行列と、に分解する分解ステップと、前記分解ステップでの分解により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断ステップと、を含み、前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、前記決定ステップでは、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σ_ｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕ_ｓと、から導出される行列Ｕ_ｓΣ_ｓＵ_ｓ ^Ｔであり、前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、前記教師基底行列は、教師基底を格納する非負値の行列であり、前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、前記抽出基底行列は、抽出基底を格納する非負値の行列であり、前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、前記分解ステップでは、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解する。

本発明のプログラムは、周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定ステップと、前記決定ステップで決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、前抽出基底行列と、抽出重み行列と、に分解する分解ステップと、前記分解ステップでの分解により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断ステップと、を含むステップをコンピュータに実行させ、前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、前記決定ステップでは、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σ_ｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕ_ｓと、から導出される行列Ｕ_ｓΣ_ｓＵ_ｓ ^Ｔであり、前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、前記教師基底行列は、教師基底を格納する非負値の行列であり、前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、前記抽出基底行列は、抽出基底を格納する非負値の行列であり、前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、前記分解ステップでは、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解する。

図１Ａは、実験の状況を示す図である。図１Ｂは、歯車箱に取り付けられている軸受を示す図である。図２は、軸受の詳細の一例を説明する図である。図３Ａは、固有値の分布の一例を示す図である。図３Ｂは、固有値の分布の一例を示す図である。図３Ｃは、固有値の分布の一例を示す図である。図３Ｄは、固有値の分布の一例を示す図である。図３Ｅは、固有値の分布の一例を示す図である。図４Ａは、計測データの波形の一例を示す図である。図４Ｂは、計測データの波形の一例を示す図である。図４Ｃは、計測データの波形の一例を示す図である。図４Ｄは、計測データの波形の一例を示す図である。図４Ｅは、計測データの波形の一例を示す図である。図５Ａは、修正された自己回帰モデルの波形の一例を示す図である。図５Ｂは、修正された自己回帰モデルの波形の一例を示す図である。図５Ｃは、修正された自己回帰モデルの波形の一例を示す図である。図５Ｄは、修正された自己回帰モデルの波形の一例を示す図である。図５Ｅは、修正された自己回帰モデルの波形の一例を示す図である。図６Ａは、修正された自己回帰モデルの周波数特性の一例を示す図である。図６Ｂは、修正された自己回帰モデルの周波数特性の一例を示す図である。図６Ｃは、修正された自己回帰モデルの周波数特性の一例を示す図である。図６Ｄは、修正された自己回帰モデルの周波数特性の一例を示す図である。図６Ｅは、修正された自己回帰モデルの周波数特性の一例を示す図である。図７Ａは、修正された自己回帰モデルの波形の一例を示す図である。図７Ｂは、修正された自己回帰モデルの波形の一例を示す図である。図７Ｃは、修正された自己回帰モデルの波形の一例を示す図である。図７Ｄは、修正された自己回帰モデルの波形の一例を示す図である。図７Ｅは、修正された自己回帰モデルの波形の一例を示す図である。図８Ａは、修正された自己回帰モデルの周波数特性の一例を示す図である。図８Ｂは、修正された自己回帰モデルの周波数特性の一例を示す図である。図８Ｃは、修正された自己回帰モデルの周波数特性の一例を示す図である。図８Ｄは、修正された自己回帰モデルの周波数特性の一例を示す図である。図８Ｅは、修正された自己回帰モデルの周波数特性の一例を示す図である。図９は、係数のパターンの一例を示す図である。図１０は、学習により特定された係数のパターンの一例を説明する図である。図１１は、各パターンの重みの一例を説明する図である。図１２は、教師基底の一例を示す図である。図１３Ａは、実験対象の軸受の一例を説明する図である。図１３Ｂは、実験対象の軸受の一例を説明する図である。図１３Ｃは、実験対象の軸受の一例を説明する図である。図１４Ａは、データ行列の分解結果の一例を説明する図である。図１４Ｂは、データ行列の分解結果の一例を説明する図である。図１５Ａは、データ行列の分解結果の一例を説明する図である。図１５Ｂは、データ行列の分解結果の一例を説明する図である。図１６Ａは、データ行列の分解結果の一例を説明する図である。図１６Ｂは、データ行列の分解結果の一例を説明する図である。図１７Ａは、データ行列の分解結果の一例を説明する図である。図１７Ｂは、データ行列の分解結果の一例を説明する図である。図１８は、診断システムのシステム構成等の一例を示す図である。図１９は、情報処理装置のハードウェア構成の一例を示す図である。図２０は、学習処理の一例を示すフローチャートである。図２１は、診断処理の一例を示すフローチャートである。図２２Ａは、教師基底と抽出基底との重みの一例を示す図である。図２２Ｂは、教師基底と抽出基底との重みの一例を示す図である。図２２Ｃは、教師基底と抽出基底との重みの一例を示す図である。図２２Ｄは、教師基底と抽出基底との重みの一例を示す図である。図２２Ｅは、教師基底と抽出基底との重みの一例を示す図である。図２３Ａは、教師基底と抽出基底との重みの一例を示す図である。図２３Ｂは、教師基底と抽出基底との重みの一例を示す図である。図２３Ｃは、教師基底と抽出基底との重みの一例を示す図である。図２３Ｄは、教師基底と抽出基底との重みの一例を示す図である。図２３Ｅは、教師基底と抽出基底との重みの一例を示す図である。図２４Ａは、教師基底と抽出基底との重みの一例を示す図である。図２４Ｂは、教師基底と抽出基底との重みの一例を示す図である。図２４Ｃは、教師基底と抽出基底との重みの一例を示す図である。図２４Ｄは、教師基底と抽出基底との重みの一例を示す図である。図２５Ａは、教師基底と抽出基底との重みの一例を示す図である。図２５Ｂは、教師基底と抽出基底との重みの一例を示す図である。図２５Ｃは、教師基底と抽出基底との重みの一例を示す図である。図２５Ｄは、教師基底と抽出基底との重みの一例を示す図である。図２６Ａは、教師基底と抽出基底との重みの一例を示す図である。図２６Ｂは、教師基底と抽出基底との重みの一例を示す図である。図２６Ｃは、教師基底と抽出基底との重みの一例を示す図である。図２６Ｄは、教師基底と抽出基底との重みの一例を示す図である。図２７Ａは、教師基底と抽出基底との重みの一例を示す図である。図２７Ｂは、教師基底と抽出基底との重みの一例を示す図である。図２７Ｃは、教師基底と抽出基底との重みの一例を示す図である。図２７Ｄは、教師基底と抽出基底との重みの一例を示す図である。

＜実施形態１＞
以下、本発明の一実施形態について図面に基づいて説明する。

（本実施形態の処理の概要）
本実施形態では、診断システムが鉄道台車に利用される軸受から計測された振動の信号に基づいて、その軸受の状態を診断する。
診断システムは、軸受の振動に応じた信号を計測する。診断システムは、計測された信号を自己回帰モデルで表現し、自己回帰モデルの係数に関する条件式における行列を特異値分解し、その行列の固有値を求める。診断システムは、求めた固有値のうち最大のものから設定された数のみを用いて、自己回帰モデルの係数を決定する。そして、診断システムは、決定した係数のパターンに基づいて、軸受の状態を診断する。

（実験１）
鉄道台車に利用されている軸受の状態診断精度の評価のために行った実験について説明する。
図１Ａは、本実験の状況を説明する図である。図１Ａの状況は、実験室に鉄道台車における車輪の駆動機構が設置されている状況である。駆動モータ１０１は、駆動モータ軸１０２、歯車継手１０３、小歯車軸１００を介して、小歯車１０５を回転させる。小歯車１０５は、回転することで、大歯車１０７を回転させる。それに合わせて、車軸１０９が回転し、車軸１０９に接続される車輪が回転することになる。本実験では、車軸１０９には、車輪の代わりにダイナモ１１０が接続されている。ダイナモ１１０は、疑似的な走行負荷を与えるために接続されている発電機である。本実験では、ダイナモ１１０を回転させるための負荷を、走行の際の負荷として想定する。

小歯車１０５が嵌められた小歯車軸１００と大歯車１０７が嵌められた車軸１０９とは、歯車箱１０４に固定されている。そして、歯車箱１０４には、小歯車軸との間に軸受１０６が装着されており、車軸１０９との間に軸受１０８が装着されている。図１Ｂは、歯車箱１０４に取り付けられている軸受１０６、１０８を側面から見た様子を示す図である。

歯車箱１０４上の振動計測位置１１１は、振動計測装置が設置される位置である。振動計測装置は、加速度センサ、レーザ変位計等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を出力する。振動計測位置１１１に設置された振動計測装置が振動計測位置１１１に伝達される振動を計測する。振動計測装置は、計測した振動を示す信号を外部の情報処理装置等に有線又は無線による通信を介して出力する。振動計測装置の計測周波数は、２０４８００Ｈｚである。ただし、他の例として、振動計測装置の計測周波数は、４０９６００Ｈｚ、１０２４００Ｈｚ等の他の周波数であることとしてもよい。本実験では、振動計測装置は、振動計測位置１１１に設置されるとするが、軸受１０６に起因する振動が伝達される位置であるならば、任意の位置に設置することとしてもよい。

ここで、軸受１０６について、説明する。図２は、軸受１０６の一例の構成を説明する図である。軸受１０６は、外輪、内輪、転動体、保持器の４つの部分から構成される。図２には、軸受１０６の外輪、内輪、転動体、保持器の概要が示されている。軸受１０６は、外輪と内輪との間に保持器により保持された転動体が挟まれる構成となっている。
本実験では、小歯車１０５に装着されている軸受１０６として、疵のない正常な軸受、内輪に疵のある軸受、保持器に疵のある軸受のそれぞれを利用した。それぞれの軸受１０６毎に、駆動モータ１０１を駆動させ、振動計測位置１１１に設置された振動計測装置が振動を計測することで、軸受１０６の状態診断を行うための計測データを取得した。本実験では、駆動モータ１０１は、軸受１０６の内輪の回転数が２９５４．７ｒｐｍ（ｒｏｕｎｄｐｅｒｍｉｎｕｔｅ）となるように駆動した。

以上のように、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合について、振動計測位置１１１における振動の計測データを取得した。

次に、振動計測装置が計測した計測データを用いた軸受１０６の状態診断の方法を説明する。以下では、振動計測装置が計測した計測データを、計測データｙとする。

本実験では、計測データｙを、修正された自己回帰モデルで近似し、この修正された自己回帰モデルの係数から特徴量を算出することとした。

ある時刻ｋ（ｋ：１≦ｋ≦Ｍの整数）における計測データｙの値をｙ_ｋとする。Ｍは、計測データｙがどの時刻までのデータを含むかを示す数であり、予め設定されている。Ｍは、例えば、後述する自己相関（式５で定義されるＲ_ｊｌ）が統計量として十分な信頼性を得るように、事前に試行錯誤して決定すればよい。ｙ_ｋを近似する自己回帰モデルは、例えば、以下の式１のようになる。式１に示すように、自己回帰モデルとは、時系列データにおけるある時刻ｋ（ｍ＋１≦ｋ≦Ｍ）のデータの予測値ｙ＾_ｋ（式において＾はｙの上に付けて表記）を、時系列データにおけるその時刻よりも前の時刻ｋ−ｌ（１≦ｌ≦ｍ）のデータの実績値ｙ_ｋ−ｌを用いて表す式である。

式１におけるαは、自己回帰モデルの係数である。また、ｍは、自己回帰モデルにおいてある時刻ｋにおける計測データｙの値であるｙ_ｋを、その時刻よりも前の過去幾つのデータを用いて近似するかを示すＭ未満の整数である。ここでは、ｍを１５００とする。

続いて、最小二乗法を用いて、自己回帰モデルによる予測値ｙ＾_ｋが実測値であるｙ_ｋに近似するための条件式を求める。自己回帰モデルによる予測値ｙ＾_ｋが実測値であるｙ_ｋに近似するための条件として、式１で与えられるｙ＾_ｋとｙ_ｋとの二乗誤差を最小化することを考える。即ち、本実験では、自己回帰モデルによる予測値ｙ＾_ｋをｙ_ｋに近似するために最小二乗法を用いる。以下の式２は、計測データと自己回帰モデルによる予測値との二乗誤差を最小にするための条件式である。

式２より、以下の式３の関係を満たす。

また、式３を変形（行列表記）することで、以下の式４のようになる。

式４におけるＲ_ｊｌは計測データｙの自己相関と呼ばれるもので、以下の式５で定義される値である。このときの｜ｊ−ｌ｜を時差という。

式４を基に、以下の自己回帰モデルの係数に関する関係式である式６を考える。式６は、自己回帰モデルによる計測データの予測値と、その予測値に対応する時刻における計測データと、の誤差を最小化する条件から導出される方程式である。式６は、ユール・ウォーカー（Ｙｕｌｅ−Ｗａｌｋｅｒ）方程式と呼ばれるものである。また、式６は自己回帰モデルの係数から成るベクトルを変数ベクトルとする線形方程式である。式６における左辺の定数ベクトルは、時差が１からｍまでの計測データの自己相関を成分とするベクトルである。以下では、式６における左辺の定数ベクトルを自己相関ベクトルとする。また、式６における右辺の係数行列は、時差が０からｍ−１までの計測データの自己相関を成分とする行列である。以下では、式６における右辺の係数行列を自己相関行列とする。

また、式６における右辺の自己相関行列（Ｒ_ｊｌで構成されるｍ×ｍの行列）を、以下の式７のように、自己相関行列Ｒと表記する。

一般に、自己回帰モデルの係数を求める際には、式６を係数αについて解くという方法が用いられる。式６では、自己回帰モデルで導出されるある時刻ｋにおける計測データの予測値ｙ＾_ｋが、その時刻ｋにおける計測データの実績値ｙ_ｋにできるだけ近づくように係数αが導出される。よって、自己回帰モデルの周波数特性に、各時刻における計測データの実績値ｙ_ｋに含まれる多数の周波数成分が含まれる。したがって、例えば、計測データｙに含まれるノイズが多い場合には、軸受１０６の振動に係る信号を抽出できなかったり、軸受１０６の故障の態様に応じた特徴を抽出することができなかったりするという問題が生じる。

そこで、本発明者らは、自己回帰モデルの係数αに乗算される自己相関行列Ｒに着目し、鋭意検討した結果、自己相関行列Ｒの固有値の一部を用いて、計測データｙに含まれるノイズの影響が低減され、軸受の振動に係る信号成分が強調される（ＳＮ比を高める）ように自己相関行列Ｒを書き換えればよいという着想に至った。

以下で、このことの具体例を説明する。
自己相関行列Ｒを特異値分解する。自己相関行列Ｒの成分は、対称であるので、自己相関行列Ｒを特異値分解すると以下の式８のように、直交行列Ｕと、対角行列Σと、直交行列Ｕの転置行列との積となる。

式８の行列Σは、式９に示すように、対角成分が自己相関行列Ｒの固有値となる対角行列である。対角行列Σの対角成分を、σ_１１、σ_２２、・・・、σ_ｍｍとする。また、行列Ｕは、各列成分ベクトルが自己相関行列Ｒの固有ベクトルとなる直交行列である。直交行列Ｕの列成分ベクトルを、ｕ_１、ｕ_２、・・・、ｕ_ｍとする。自己相関行列Ｒの固有ベクトルｕ_ｊに対する固有値がσ_ｊｊという対応関係が有る。自己相関行列Ｒの固有値は、自己回帰モデルによる計測データの予測値の時間波形に含まれる各周波数の成分の強度に反映する変数である。

自己相関行列Ｒの特異値分解の結果得られる対角行列Σの対角成分であるσ_１１、σ_２２、・・・、σ_ｍｍの値は、数式の表記を簡略にするために降順とする。これらの自己相関行列Ｒの固有値のうち、最大のものから１以上且つｍ未満の設定された数である使用固有値数ｓ個の固有値を用いて、以下の式１０のように、行列Ｒ’を定義する。行列Ｒ’は、自己相関行列Ｒの固有値のうち使用固有値数ｓ個の固有値用いて自己相関行列Ｒを近似した行列である。

式１０における行列Ｕ_ｓは、式８の直交行列Ｕの左からｓ個の列成分ベクトル（使用される固有値に対応する固有ベクトル）により構成されるｍ×ｓ行列である。つまり、行列Ｕ_ｓは、直交行列Ｕから左のｍ×ｓの成分を切り出して構成される部分行列である。また、Ｕ_ｓ ^ＴはＵ_ｓの転置行列であり、式８の行列Ｕ^Ｔの上からｓ個の行成分ベクトルにより構成されるｓ×ｍ行列である。式１０における行列Σ_ｓは、式８の対角行列Σの左からｓ個の列と上からｓ個の行により構成されるｓ×ｓ行列である。つまり、行列Σ_ｓは、対角行列Σから左上のｓ×ｓの成分を切り出して構成される部分行列である。
行列Σ_ｓ及び行列Ｕ_ｓを行列成分表現すれば、以下の式１１のようになる。

自己相関行列Ｒの代わりに行列Ｒ’を用いることで、式６の関係式を、以下の式１２のように書き換える。

式１２を変形することで、係数αを求める式１３が得られる。式１３によって求められた係数αを用いて、式１により予測値ｙ＾_ｋを算出するモデルを「修正された自己回帰モデル」とする。

これまで対角行列Σの対角成分であるσ_１１、σ_２２、・・・、σ_ｍｍの値を降順として説明したが、係数αの算出過程において対角行列Σの対角成分は降順である必要はない。対角行列Σの対角成分は降順でない場合には、行列Ｕ_ｓは直交行列Ｕから左のｍ×ｓの成分を切り出すのではなく、使用される固有値に対応する列成分ベクトル（固有ベクトル）を切り出して構成される部分行列である。また、対角行列Σの対角成分は降順でない場合には、行列Σ_ｓは対角行列Σから左上のｓ×ｓの成分を切り出すのではなく、使用される固有値を対角成分とするように切り出される部分行列である。

式１３は、修正された自己回帰モデルの係数の決定に利用される方程式である。式１２、１３の行列Ｕ_ｓは、自己相関行列Ｒの特異値分解により得られる直交行列の部分行列であって、利用される固有値に対応する固有ベクトルを列成分ベクトルとする行列である第３の行列である。また、式１２、１３の行列Σ_ｓは、自己相関行列Ｒの特異値分解により得られる対角行列の部分行列であって、利用される固有値を対角成分とする行列である第２の行列である。そして、式１２、１３の行列Ｕ_ｓΣ_ｓＵ_ｓ ^Ｔは、行列Σ_ｓと行列Ｕ_ｓとから導出される行列である第１の行列である。

式１３の右辺を計算することにより、修正された自己回帰モデルの係数αが求まる。以上に、修正された自己回帰モデルの係数の導出方法の一例について説明してきたが、その基となる自己回帰モデルの係数の導出については直感的に分かり易いように予測値に対して最小二乗法を用いる方法で説明した。しかしながら、一般的には確率過程という概念を用いて自己回帰モデルを定義し、その係数を導出する方法が知られている。その場合に、自己相関は確率過程（母集団）の自己相関で表現される。この確率過程の自己相関は時差の関数として表されるものである。従って、本実施形態における計測データの自己相関は、確率過程の自己相関を近似するものであれば他の計算式で算出した値に代えても良い。例えば、Ｒ_２２〜Ｒ_ｍｍは時差が０の自己相関であるが、これらをＲ_１１に置き換えてもよい。

本実験では、以下の４つの場合のそれぞれについて、得られた計測データ毎に、式５と式７を用いて自己相関行列Ｒを求め、式８で表される特異値分解を行うことによって自己相関行列Ｒの固有値を求め、自己相関行列Ｒの固有値の分布を求めた。４つの場合は、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、および軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合である。

次に、求めた係数αで特定される修正された自己回帰モデルの周波数特性を求める方法を説明する。
式１で与えられる計測データの予測値ｙ＾_ｋの予測誤差ｘ_ｋが白色雑音になるという性質を使うと、修正された自己回帰モデルは白色雑音である予測誤差ｘ_ｋを入力とし、計測データの実測値ｙ_ｋを出力とする線形時不変システムとみなせる。従って、以下の手順で周波数特性の算出式を導出することができる。以下では、修正された自己回帰モデルを線形時不変モデルとみなしたものを、単にシステムと称することにする。予測誤差ｘ_ｋは、以下の式１４で表される。

式１４の両辺にｚ変換を施すと、以下の式１５が得られる。

式１５から、システムのインパルス応答のｚ変換である伝達関数Ｈ（ｚ）は、以下の式１６のように求まる。

システムの周波数特性は、正弦波入力に対する出力の振幅と位相との変化として現れ、インパルス応答のフーリエ変換で求められる。言い換えると、ｚが複素平面の単位円上を回転した時の伝達関数Ｈ（ｚ）が周波数特性となる。ここで、式１６におけるｚを、以下の式１７のように置くことを考える。

ここで、ｊは虚数単位、ωは角周波数、Ｔは標本間隔である。

その場合、Ｈ（ｚ）の振幅特性（システムの周波数特性）は、以下の式１８のように表すことができる。式１８におけるωＴは、０〜２πの範囲で変化するものとする。

本実験では、得られた計測データ毎に、式５と式７を用いて自己相関行列Ｒを求め、式８で表される特異値分解を行うことによって自己相関行列Ｒの固有値を求め、自己相関行列Ｒの固有値の分布を求めた。また、得られた計測データ毎に、計測データの波形を求めた。また、得られた計測データ毎に、上記手順で求められた自己相関行列Ｒの特異値分解から式１１を用いて行列Σ_ｓと行列Ｕ_ｓとを求め、式５と式１３を用いて修正された自己回帰モデルの係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形を求めた。そして、得られた計測データ毎に、式１８を用いて、修正された自己回帰モデルの周波数特性等を求めた。
求められた実験結果を、図３Ａ〜図８Ｅを用いて説明する。

図３Ａ〜図３Ｅは、計測データそれぞれについて、上記手順で求められた自己相関行列Ｒの固有値の分布を示す図である。
図３Ａのグラフは、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３Ｂのグラフは、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３Ｃのグラフは、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３Ｄのグラフは、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３Ｅのグラフは、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。

図３Ａ〜図３Ｅの各グラフは、自己相関行列Ｒを特異値分解して得られた固有値σ_１１〜σ_ｍｍを昇順に並べ替え、プロットしたグラフであり、横軸が固有値のインデックス、縦軸が固有値の値を示す。

図３Ａを見ると、他よりも顕著に高い値をもつ固有値が５つあるのが分かる。その５つの固有値の中でも特に２つの固有値が他の３つよりも高い値となっていることが分かる。
何れの場合にも、自己相関行列Ｒの固有値のうち値が他よりも顕著に高い固有値の数は、固有値全体の数よりも顕著に少ないことが分かる。

また、本実験では、得られた計測データ毎に、計測データの波形を求めた。
図４Ａ〜図４Ｅは、計測データそれぞれの波形を示す図である。
図４Ａのグラフは、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データの波形を示すグラフである。図４Ａのグラフでは、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４Ｂのグラフは、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフである。図４Ｂのグラフでは、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４Ｃのグラフは、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフである。図４Ｃのグラフでは、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４Ｄのグラフは、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフである。図４Ｄのグラフでは、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４Ｅのグラフは、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフである。図４Ｅのグラフでは、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。

図４の各グラフを見ると、例えば、図４Ａ及び図４Ｃが似通った波形をしているのが分かる。また、例えば、図４Ｄ及び図４Ｅが似通った波形をしているのも分かる。そのため、図４Ａ〜図４Ｅのグラフから軸受１０６に生じた異常の状態を診断することは困難である。

図５Ａ〜図５Ｅは、使用固有値数ｓを１として求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示す図である。

図５Ａのグラフは、使用固有値数ｓを１として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図５Ａのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５Ｂのグラフは、使用固有値数ｓを１として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図５Ｂのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５Ｃのグラフは、使用固有値数ｓを１として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図５Ｃのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５Ｄのグラフは、使用固有値数ｓを１として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図５Ｄのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５Ｅのグラフは、使用固有値数ｓを１として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図５Ｅのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図５Ａ〜図５Ｅの各グラフは、図４Ａ〜図４Ｅの各グラフと比べると、ノイズ成分が低減されていることが分かる。

図６Ａ〜図６Ｅは、使用固有値数ｓを１として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図６Ａ〜図６Ｅの各グラフは、図５Ａ〜図５Ｅの各波形についての周波数特性を示すグラフである。

図６Ａのグラフは、使用固有値数ｓを１として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図６Ａのグラフでは、横軸が周波数、縦軸が信号強度を示す。図６Ｂのグラフは、使用固有値数ｓを１として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図６Ｂのグラフでは、横軸が周波数、縦軸が信号強度を示す。図６Ｃのグラフは、使用固有値数ｓを１として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図６Ｃのグラフでは、横軸が周波数、縦軸が信号強度を示す。図６Ｄのグラフは、使用固有値数ｓを１として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図６Ｄのグラフでは、横軸が周波数、縦軸が信号強度を示す。図６Ｅのグラフは、使用固有値数ｓを１として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。る。図６Ｅのグラフでは、横軸が周波数、縦軸が信号強度を示す。

図６Ａのグラフを見ると、周波数８２４Ｈｚの部分にピークが立っているのが分かる。また、図６Ｂのグラフを見ると、周波数１６４９Ｈｚの部分にピークが立っているのが分かる。また、図６Ｃのグラフを見ると、周波数１６４６Ｈｚの部分にピークが立っているのが分かる。

図６Ｄのグラフを見ると、周波数１９８５３Ｈｚの部分にピークが立っているのが分かる。また、図６Ｅのグラフを見ると、周波数２３００７Ｈｚの部分にピークが立っているのが分かる。

図６Ａ〜図６Ｅの各グラフにおけるピークが立っている周波数を見ると、軸受１０６に疵がない場合は、軸受１０６に疵がある場合に比べて低い周波数でピークが立っていることが分かる。即ち、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６における疵の有無を判断することが可能であることが分かる。また、図６Ｄよりも、図６Ｅの方が高い周波数であることが分かる。また、図６Ｄ〜図６Ｅの各グラフにおけるピークが立っている周波数は、図６Ａ〜図６Ｃの各グラフにおけるピークが立っている周波数よりも顕著に高いことが分かる。そのため、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６の転動体又は外輪の疵の有無を判断することが可能であることが分かる。

以上より、自己相関行列Ｒの１個の固有値に基づいて求められた係数αで特定される修正された自己回帰モデルの周波数特性に基づいて、軸受１０６の状態診断を行うことができるという知見を得ることができた。

しかし、図６Ｂと図６Ｃとのグラフに示されるように、軸受１０６の内輪に疵がある場合と、軸受１０６の保持器に疵がある場合と、では、似通った値の周波数において周波数特性のピークが立っている。そのため、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６に内輪の又は保持器の疵の有無を判断することは可能であるが、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を判断することは困難である。

そこで、本実験では、更に、使用固有値数ｓの値を５に増やして、自己相関行列Ｒの固有値のうち、最大のものから５つを用いて係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形、周波数特性を求めた。

図７Ａ〜図７Ｅは、使用固有値数ｓを５として求められた係数αで特定される修正された自己回帰モデルの波形を示す図である。

図７Ａのグラフは、使用固有値数ｓを５として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図７Ａのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７Ｂのグラフは、使用固有値数ｓを５として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図７Ｂのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７Ｃのグラフは、使用固有値数ｓを５として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図７Ｃのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７Ｄのグラフは、使用固有値数ｓを５として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図７Ｄのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７Ｅのグラフは、使用固有値数ｓを５として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図７Ｅのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。

図８Ａ〜図８Ｅは、使用固有値数ｓを５として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図８Ａ〜図８Ｅの各グラフは、図７Ａ〜図７Ｅの各波形についての周波数特性を示すグラフである。

図８Ａのグラフは、使用固有値数ｓを５として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図８Ａのグラフでは、横軸が周波数、縦軸が信号強度を示す。図８Ｂのグラフは、使用固有値数ｓを５として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図８Ｂのグラフでは、横軸が周波数、縦軸が信号強度を示す。図８Ｃのグラフは、使用固有値数ｓを５として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図８Ｃのグラフでは、横軸が周波数、縦軸が信号強度を示す。図８Ｄのグラフは、使用固有値数ｓを５として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図８Ｄのグラフでは、横軸が周波数、縦軸が信号強度を示す。図８Ｅのグラフは、使用固有値数ｓを５として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図８Ｅのグラフでは、横軸が周波数、縦軸が信号強度を示す。

図８Ａのグラフを見ると、周波数８２１Ｈｚと１０４７Ｈｚの部分にピークが立っているのが分かる。また、図８Ｂのグラフを見ると、周波数１６４８Ｈｚ、２０２１Ｈｚ、６４７４Ｈｚの部分にピークが立っているのが分かる。また、図８Ｃのグラフを見ると、周波数１６４６Ｈｚ、３２９１Ｈｚ、１７４６Ｈｚの部分にピークが立っているのが分かる。図８Ｄのグラフを見ると、周波数１９８５３Ｈｚの部分にピークが立っているのが分かる。また、図８Ｅのグラフを見ると、周波数２２７８５Ｈｚ、２３０２０Ｈｚの部分にピークが立っているのが分かる。

図８Ａ〜図８Ｅの各グラフにおけるピークと、図６Ａ〜図６Ｅの各グラフにおけるピークと、を見比べてみると、図６Ａ〜図６Ｅのグラフでピークが立っている周波数近辺においては、図８Ａ〜図８Ｅのグラフでもピークが立っているのが分かる。図８Ａ〜図８Ｅのグラフでは、更に、図６Ａ〜図６Ｅのグラフで見られなかったピークが見られるようになっている。

このように、係数αを求めるのに利用される自己相関行列Ｒの固有値を増やすことで、修正された自己回帰モデルによる予測値において、図６Ａ〜図６Ｅでピークが立っていた周波数以外の周波数の成分の信号強度が増加したことが分かる。即ち、係数αを求めるのに利用される自己相関行列Ｒの固有値は、求められた係数αで特定される修正された自己回帰モデルによる計測データの予測値の時間波形に含まれる各周波数の成分の強度と相関があることが見出された。

また、図８Ｂ、図８Ｃのグラフを見てみると、図６Ｂ、図６Ｃのグラフと同様に、１６４０Ｈｚ付近にピークが立っているのが分かる。しかし、図８Ｂでは、更に、２０２１Ｈｚ、６４７４Ｈｚにおいてもピークが立っており、図８Ｃでは、更に、３２９１Ｈｚ、１７４６Ｈｚにおいてもピークが立っていることが分かる。このように、図８Ｂ〜図８Ｃでは、内輪に疵のある軸受１０６と保持器に疵のある軸受１０６とで、ピークが立っている周波数に顕著な違いが見て取れる。即ち、使用固有値数ｓを５にすることで、軸受の疵が、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を見分けることができる。

このように本発明者らは、使用固有値数ｓが１の場合には見分けることが困難であった疵の存在する部位を、使用固有値数ｓを増やすことで、見分けることができるようになるという知見を得た。一方で、本発明者らは、使用固有値数ｓを増やし過ぎると修正された自己回帰モデルが、実際の計測データに近づいていくことによって、軸受１０６の状態診断に有用でないノイズ成分等についても信号強度が増加されてしまい、式１８で得られる周波数特性にもノイズ成分に起因するピークが生じるという知見を得た。即ち、自己相関行列Ｒの固有値のうち、相対的に値の高い一部の固有値を用いて係数αを求めることで、計測データに含まれる成分のうち、物体の状態診断に有用な成分を近似した修正された自己回帰モデルを求めることができると考えられる。

以上より、係数αを求めるのに利用される自己相関行列Ｒの固有値の数である使用固有値ｓを調整することで、より精度よく、軸受１０６の状態診断を行うことができるという知見を得ることができた。

（実験２）
また、発明者らは、実験１と同様の実験状況において、以下のような実験を行った。軸受１０６が疵のない正常な状態である場合、内輪に疵がある状態である場合、保持器に疵がある状態である場合、転動体に疵がある状態である場合、外輪に疵がある状態である場合それぞれについて、以下の作業をｑ回行った。ｑは５０としたが、３０や１００等の他の値としてもよい。即ち、軸受１０６について、連続したＭ個の計測データを取得し、取得した連続したＭ個の計測データ毎に、式１１を用いて行列Σ_ｓと行列Ｕ_ｓとを求め、式５と式１３を用いて修正された自己回帰モデルの係数αを求めた。使用固有値数ｓは、３とした。

これにより、軸受１０６が正常な状態である場合、内輪に疵がある状態である場合、保持器に疵がある状態である場合、転動体に疵がある状態である場合、外輪に疵がある状態である場合の５つの場合それぞれについて、ｑ個の修正された自己回帰モデルの係数αを取得した。

軸受１０６が取り得る状態は、疵のない正常な状態、内輪に疵がある状態、保持器に疵がある状態、転動体に疵がある状態、および外輪に疵がある状態の５つの状態であるとする。以下では、軸受１０６が取り得る状態の数をｐとする。ｐは、これらの５つの状態に対応して５とした。

図９は、軸受１０６が正常な状態である場合、内輪に疵がある状態である場合、保持器に疵がある状態である場合、転動体に疵がある状態である場合、軸受１０６の外輪に疵がある場合それぞれについての、上記手順で求められた係数αのパターンを示す図である。係数αのパターンとは、係数αの成分を並べた際の形状を示す情報である。ｍの値は、５００とした。そのため、係数αは、α_１〜α_５００の５００個の成分から成る。図９に示されるパターンそれぞれは、α_ｉ（１＜＝ｉ＜＝５００）をインデックスｉ順に並べたパターンである。このインデックスｉは、式１において係数の成分がかかる計測データｙの計測時点を示すインデックスとみなすことができる。即ち、図９に示されるパターンそれぞれは、係数αの各成分を、式１において各成分がかかる計測データｙの計測時点に基づいて並べたパターンの一例である。

係数αのパターンには、例えば、係数αの成分をインデックス順に並べたベクトル、係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αのパターンには、例えば、各成分に同じ値が加算された係数αの成分をインデックス順に並べたベクトル、各成分に同じ値が加算された係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、各成分に同じ値が加算された係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αのパターンには、例えば、各成分に同じ値が掛けられた係数αの成分をインデックス順に並べたベクトル、各成分に同じ値が掛けられた係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、各成分に同じ値が掛けられた係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αの成分を並べた際の形状を示す情報である係数αのパターンは、ベクトルでなく、係数αの成分を並べたグラフの画像であるとしてもよい。

図９に示されるパターンを見てみると、以下のようなことが見てとれる。軸受１０６の外輪に疵がある場合、ｑ個の係数αのパターンそれぞれは、上下対称な２つの波線で囲まれた領域を示すようなパターンとなっている。また、軸受１０６の転動体に疵がある場合、ｑ個の係数αのパターンそれぞれは、複数の波線が重なり合っているようなパターンとなっている。また、軸受１０６の内輪に疵がある場合、ｑ個の係数αのパターンそれぞれは、４つのピークのある山形のパターンが連続しているパターンとなっている。また、軸受１０６に疵がない場合、ｑ個の係数αのパターンそれぞれは、徐々に減衰していく波型のパターンとなっている。また、軸受１０６の保持器に疵がある場合、ｑ個の係数αのパターンそれぞれは、２つのピークのある山形のパターンが連続しているパターンとなっている。

このように、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合それぞれについて、修正された自己回帰モデルの係数αが、特徴的なパターンを示すことが見て取れる。

そこで、発明者らは、修正された自己回帰モデルの係数αのパターンが、どのようなパターンを示すのかを特定することで、軸受１０６の診断を行うことができるという知見を得た。

また、本実験では、以下のようにして、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合それぞれについて、修正された自己回帰モデルの係数αの代表的なパターンを学習により特定した。

まず、各列が図９に示される係数αそれぞれを示すｍ×ｄのサイズの行列Ｙを、以下の式１９のように生成する。ｄは、ｐ×ｑであり、２５０となる。

式１９のα_{ｋ，（ｉ，ｊ）}は、ｉに対応する状態の軸受１０６から計測されたｑ個の計測データのうちｊ番目の計測データから決定された係数αの第ｋ成分α_ｋを示す。本実験では、外輪に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるα_{ｋ，（０，ｊ）}（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、転動体に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるα_{ｋ，（１，ｊ）}（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、内輪に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるα_{ｋ，（２，ｊ）}（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、疵がない正常な状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるα_{ｋ，（３，ｊ）}（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、保持器に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるα_{ｋ，（４，ｊ）}（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。ここでは、ｐは５である。式１９では、表記の都合で、α_{ｋ，（４，ｊ）}をα_{ｋ，（ｐ−１，ｊ）}とする。また、式１９では、α_{ｋ，（２，ｊ）}およびα_{ｋ，（３，ｊ）}の表記を省略する。

本実験では、この行列Ｙを非負値行列因子分解することで、軸受１０６の状態それぞれについて、修正された自己回帰モデルの係数αの代表的なパターンを特定する。行列Ｙを非負値行列因子分解するためには、行列Ｙの各列の係数αのパターンを変えずに、行列Ｙ自体の各成分を、非負値にするよう調整する必要がある。そこで、例えば、行列Ｙの０未満の成分のうち、最小のものの絶対値を、行列Ｙの各成分に加算することで、行列Ｙの各成分が非負値となるようにする。以下では、各成分が非負値である行列を、非負値行列とする。

以下の式２０に示すように、この行列Ｙを、非負値行列である行列Ａと行列Ｐとに非負値行列因子分解する。行列Ａは、ｍ×ｐのサイズの行列である。行列Ｐは、ｐ×ｄのサイズの行列である。行列Ａの列数は、非負値行列因子分解における基底ベクトルの数であり、本実験では、軸受１０６が取り得る状態の数としてｐ（５）とする。これにより、行列Ａの各列には、行列Ｙから導出されたｐ（５）個の係数αのパターン（基底ベクトルともいう）が格納されることとなる。

以下の式２１のようにコスト関数Ｊを定義する。

式２１の右辺第１項は、行列Ａと行列Ｐとの積と、行列Ｙと、が、これらの残差の二乗和を最小にするという意味で近似的に一致するための条件を示す。式２１の右辺第１項の右下のＦｒ記号は、フロベニウスノルムであることを示す記号である。即ち、式２１の右辺第１項内の｜｜ＡＰ−Ｙ｜｜は、フロベニウスノルムであり、行列（ＡＰ−Ｙ）の各成分の絶対値の２乗の合計の平方根となる。式２１の右辺第２項は、行列Ｐのスパース性を要求する条件を示す。この条件を付加することによって、非負値行列因子分解が過学習となることを防ぐことができる。

式２１の右辺第２項におけるインデックスｉは、行列の行を示すインデックスである。また、式２１の右辺第２項におけるインデックスｋは、行列の列を示すインデックスである。式２１の右辺第３項は、行列Ａの異なる列に格納されるパターン同士の直交性を近似的に要求する条件を示す。式２１の右辺第３項によって、行列Ａのｐ（５）個の列に格納されるパターンが可能な限り線形従属となることを防ぐ効果が得られる。行列Ａの列に格納されるパターンが線形従属となれば、あるパターンを行列Ａの列に格納されるパターンの線形結合で近似するときの係数（重みともいう）が一通りではなくなるため、パターンの同定が難しくなる。また、式２１の右辺第３項における行列Ａ^Ｔは、行列Ａの転置行列であり、行列Ｉは、単位行列である。また、式２１の右辺第３項におけるインデックスｋ、ｋ⁻（式中では、上にバーが引かれたｋとして表現）は、それぞれ行列の行と、列と、を示すインデックスである。

ここで、コスト関数Ｊとして、式２１の右辺第１項のみを用いても、非負値行列因子分解は機能するので、コスト関数Ｊとして、式２１の右辺第１項のみを用いてもよい。また、上述の過学習や線形従属といった問題が生じる場合に、コスト関数Ｊとして、式２１の右辺第１項に式２１の右辺第２項、第３項を追加した式を用いることとしてもよい。そして、右辺第２項、第３項を追加する場合、軸受１０６の各状態と、図１４Ａ〜図１４Ｂ等で後述する方法で決定される学習パターンと、を良好に対応付けることが出来るように、例えば、ケーススタディにより係数βと係数γとを予め決定しておく。本実施形態では、βを０．０３５、γを０とした。

このコスト関数Ｊを最小化するように、行列Ｙを、行列Ａと行列Ｐと分解する。即ち、以下の式２２で定義されるコスト関数Ｊの最小化問題を解くこととなる。

式２２における第１の制約条件は、行列Ａと行列Ｐとが非負値行列であることを示す条件である。式２２における第２の制約条件は、行列Ａを基準化（ｎｏｒｍａｌｉｚｅ）するための条件である。また、式２２の第２の制約条件におけるインデックスｋは、行列の行と、列と、を示すインデックスである。
本実験では、コスト関数Ｊを最小化するために、以下のような処理を行う。即ち、以下の式２３と式２４とを用いて行列Ａと行列Ｐとの更新を繰り返すことで、コスト関数Ｊを最小化する。また、コスト関数Ｊを最小化するために、以下に記載する方法以外の方法（例えば、最急降下法、ニュートン法、確率的勾配降下法等の最適化手法）を用いてもよい。

式２３のη^Ｐは、行列Ｐの更新に用いられる緩和係数である。式２４のη^Ａは、行列Ａの更新に用いられる緩和係数である。
式２３におけるコスト関数Ｊの行列Ｐによる偏微分は、以下の式２５の通りとなる。

式２５の行列Ｅ_１は、サイズがｐ×ｄの全ての成分が１の行列である。
式２４におけるコスト関数Ｊの行列Ａによる偏微分は、以下の式２６の通りとなる。

式２６の行列Ｅ_２は、サイズがｐ×ｐの全ての成分が１の行列である。
式２５を用いて、式２３を書き換えると以下の式２７が求まる。

式２７のインデックスｉは、行列の行を示すインデックスである。式２７のインデックスｋは、行列の列を示すインデックスである。

行列Ｐは、非負値であるため、式２７の値は、非負値である必要がある。式２７の右辺の第３項は、必ず非負値となる項である。即ち、式２７の右辺から第３項を除いた式の値が０であるならば、必ず式２７は、非負値となる。そのため、以下の式２８に示すように、式２７の右辺から第３項を除いた式が０であるという条件を満たすように、行列Ｐの更新を行うこととする。

式２８のインデックスｉは、行列の行を示すインデックスである。式２８のインデックスｋは、行列の列を示すインデックスである。
式２８から、式２３における緩和係数η^Ｐの値が、以下の式２９のように求まる。

式２９のインデックスｉは、行列の行を示すインデックスである。式２９のインデックスｋは、行列の列を示すインデックスである。
式２９を用いて、式２３を書き換えると、以下の式３０のようになる。

式３０のインデックスｉは、行列の行を示すインデックスである。式３０のインデックスｋは、行列の列を示すインデックスである。この式３０を用いて、行列Ｐの各成分を更新することで、行列Ｐを更新していくこととなる。
また、式２６を用いて、式２４を書き換えると以下の式３１が求まる。

式３１のインデックスｋは、行列の行を示すインデックスである。式２７のインデックスｊは、行列の列を示すインデックスである。
行列Ａは、非負値であるため、式３１の値は、非負値である必要がある。式３１の右辺の第３項は、必ず非負値となる項である。即ち、式３１の右辺から第３項を除いた式の値が０であるならば、必ず式２７は、非負値となる。そのため、以下の式３２に示すように、式３１の右辺から第３項を除いた式が０であるという条件を満たすように、行列Ａの更新を行うこととする。

式３２のインデックスｋは、行列の行を示すインデックスである。式３２のインデックスｊは、行列の列を示すインデックスである。
式３２から、式２４における緩和係数η^Ａの値が、以下の式３３のように求まる。

式３３のインデックスｋは、行列の行を示すインデックスである。式３３のインデックスｊは、行列の列を示すインデックスである。
式３３を用いて、式２４を書き換えると、以下の式３４のようになる。

式３４のインデックスｋは、行列の行を示すインデックスである。式３４のインデックスｊは、行列の列を示すインデックスである。この式３４を用いて、行列Ａの各成分を更新することで、行列Ａを更新していくこととなる。
また、式２２における第２の制約条件を満たすようにするため、以下の式３５に示すように、行列Ａを更新する。

式３５のｋは、行列の行を示すインデックスである。式３５のｊは、行列の行又は列を示すインデックスである。
コスト関数Ｊの値が、収束するまで、式３０を用いて行列Ｐを更新し、式３４、式３５を用いて行列Ａを更新することを繰り返す。例えば、各更新段階において、更新後の行列Ｐ及び行列Ａから、それぞれ更新前の行列Ｐ及び行列Ａを引いた差である行列△Ｐ及び行列△Ａを求める。そして、行列△Ｐ及び行列△Ａそれぞれのフロベニウスノルムを計算する。各更新段階における行列△Ｐ及び行列△Ａのフロベニウスノルムの、それぞれ最初の更新段階における行列△Ｐ及び行列△Ａのフロベニウスノルムに対する比が、予め定められた閾値以下となる更新段階をもって、コスト関数Ｊの値が収束したとする。そして、その際の行列Ａと行列Ｐとをコスト関数Ｊを最小化させる行列とする。

本実験では、以上のように、コスト関数Ｊを最小化させる行列Ａと行列Ｐとを求めることとした。ただし、他の例として、行列Ｐと行列Ａとの更新前後で、コスト関数Ｊの値の変動が予め定められた閾値以下となることが、予め定められた回数連続したら、コスト関数Ｊの値が収束したとすることとしてもよい。そして、その際の行列Ａと行列Ｐとをコスト関数Ｊを最小化させる行列とすることとしてもよい。

コスト関数Ｊが収束した際の行列Ａの各列が示すパターンを確認した。図１０に、コスト関数Ｊが収束した際の行列Ａの各列が示すパターンの一例を示す。行列Ａ列には、列毎に、ｍ（５００）個の値が含まれ、係数αのパターンを示す。
図１０の一番上のグラフは、行列Ａの１列目のｍ個の値（Ａ_ｉ，１（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の上から２つ目のグラフは、行列Ａの２列目のｍ個の値（Ａ_ｉ，２（１＜＝ｉ＜＝ｍ））を、行を示すインデックス（ｉ）順に並べたパターンを示す。図１０の上から３つ目のグラフは、行列Ａの３列目のｍ個の値（Ａ_ｉ，３（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の下から２つ目のグラフは、行列Ａの４列目のｍ個の値（Ａ_ｉ，４（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の一番下のグラフは、行列Ａの５列目（ｐ列目）のｍ個の値（Ａ_ｉ，５（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。

また、コスト関数Ｊが収束した際の行列Ｐの各成分の値を確認した。式２０から得られる以下の式３６に示されるように、行列Ｐの各列（第ｊ列Ｐ_ｋ，ｊ（１＜＝ｋ＜＝ｐ））は、行列Ｙの各列（第ｊ列Ｙ_ｉ，ｊ（１＜＝ｉ＜＝ｍ））に示されるパターンを行列Ａの各列に示されるパターン（基底ベクトル）の線形結合で近似するときの係数（重み）を示している。

図１１に、コスト関数Ｊが収束した際の行列Ｐの各行の値を示すグラフを示す。図１１のグラフの横軸は、行列Ｐの列のインデックスを示す。図１１のグラフの縦軸は、行列Ｐの各成分の値を示す。

図１１のライン１１０１は、行列Ｐの１行目のｄ（２５０）個の値を結んだラインである。ライン１１０１は、横軸が０以上、５０未満の範囲（１列目からｑ（５０）列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。

図１１のライン１１０２は、行列Ｐの２行目のｄ個の値を結んだラインである。ライン１１０２は、横軸が５０以上、１００未満の範囲（ｑ＋１列目から２ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。

図１１のライン１１０３は、行列Ｐの４行目のｄ個の値を結んだラインである。ライン１１０３は、横軸が１００以上、１５０未満の範囲（２ｑ＋１列目から３ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。

図１１のライン１１０４は、行列Ｐの５行目のｄ個の値を結んだラインである。ライン１１０４は、横軸が１５０以上、２００未満の範囲（３ｑ＋１列目から４ｑ列目までの範囲）では、値が約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。

図１１のライン１１０５は、行列Ｐの３行目のｄ個の値を結んだラインである。ライン１１０５は、横軸が２００以上、２５０未満の範囲（４ｑ＋１列目から５ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。

即ち、行列Ｐの各列の値は、１列目からｑ（５０）列目までの範囲では、１行目の値が約１２となり、２〜５行目の値が約０となる。また、行列Ｐの各列の値は、ｑ＋１列目から２ｑ列目までの範囲では、２行目の値が約１２となり、１行目、３〜５行目の値が約０となる。また、行列Ｐの各列の値は、２ｑ＋１列目から３ｑ列目までの範囲では、４行目の値が約１２となり、１〜３行目、５行目の値が約０となる。また、行列Ｐの各列の値は、３ｑ＋１列目から４ｑ列目までの範囲では、５行目の値が約１２となり、１〜４行目の値が約０となる。また、行列Ｐの各列の値は、４ｑ＋１列目から５ｑ列目までの範囲では、３行目の値が約１２となり、１、２、４、５行目の値が約０となる。
そのため、行列ＡＰの１列目からｑ列目までの各列のデータが示すパターンは、図１０の一番上のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の一番上のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰのｑ＋１列目から２ｑ列目までの各列のデータが示すパターンは、図１０の上から２つ目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の上から２つ目のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの２ｑ＋１列目から３ｑ列目までの各列のデータが示すパターンは、図１０の下から２つ目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の下から２つ目のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの３ｑ＋１列目から４ｑ列目までの各列のデータが示すパターンは、図１０の一番下のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の一番下のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの４ｑ＋１列目から５ｑ列目までの各列のデータが示すパターンは、図１０の上から３番目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の上から３つ目のグラフに示されるパターンに類似するパターンとなる。

行列Ｙの１列目からｑ列目までの各列のデータが示すパターンは、図９に示す外輪に疵がある軸受１０６に対応するパターンである。また、行列Ｙのｑ＋１列目から２ｑ列目までの各列のデータが示すパターンは、図９に示す転動体に疵がある軸受１０６に対応するパターンである。また、行列Ｙの２ｑ＋１列目から３ｑ列目までの各列のデータが示すパターンは、図９に示す内輪に疵がある軸受１０６に対応するパターンである。また、行列Ｙの３ｑ＋１列目から４ｑ列目までの各列のデータが示すパターンは、図９に示す疵のない正常な軸受１０６に対応するパターンである。また、行列Ｙの４ｑ＋１列目から５ｑ列目までの各列のデータが示すパターンは、図９に示す保持器に疵がある軸受１０６に対応するパターンである。

行列Ｙと行列Ａとを見比べてみると、双方の各列のデータが示すパターンは、類似するパターンとなっていることが分かる。即ち、行列Ａの各列のデータが示すパターンは、それぞれ、各状態の軸受１０６に対応する係数αのパターンとなっていることが分かる。行列Ａの１列目のデータは、外輪に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの２列目のデータは、転動体に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの３列目のデータは、保持器に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの４列目のデータは、内輪に疵のある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの５列目のデータは、疵のない正常な状態の軸受１０６に対応する係数αの代表的なパターンを示す。

以上より、行列Ｙを生成し、生成した行列Ｙを、式２２を満たすように、非負値行列因子分解により学習することで、行列Ａの各列（基底ベクトル）として、各状態の軸受１０６に対応する係数αのパターンが特定されたことが確認された。

（実験３）
また、発明者は、軸受１０６の状態を診断する方法の着想を得た。この着想の概要について説明する。状態が既知な軸受１０６についての計測データから、予めその状態に対応する修正された自己回帰モデルの係数αのパターンを示す非負値の行列を、教師基底行列として特定する。そして、状態が不知な軸受１０６についての計測データから修正された自己回帰モデルの係数αを予め定められた個数だけ求めて、求めたその個数の係数αそれぞれのパターンを示す非負値のデータを各列に格納する行列を、教師基底行列と、教師基底行列の重み行列と、この非負値のデータを各列に格納する行列から抽出される基底を示す抽出基底行列と、抽出基底行列の重み行列と、に分解する。そして、教師基底行列の重み行列と、抽出基底行列の重み行列と、を比較することで、状態が不知な軸受１０６の状態が、教師基底行列に対応する状態であるか否かを診断する。以上が、発明者が得た着想の概要である。

そこで、発明者らは、この着想で得られた方法の実効性を検証するため、実験１と同様の実験状況において、以下のような実験を行った。本実験では、使用固有値数ｓは、３とした。また、ｍは、５００とした。

実験２で説明した方法で、疵のない正常な状態の軸受１０６に対応する係数αのパターンを示すデータ（式２０の行列Ａの正常な軸受１０６に対応する列のデータ）を、教師基底に属する基底ベクトルとして取得した。取得した正常な状態の軸受１０６に対応する教師基底に属するパターンは、図１２に示すようなパターンとなった。図１２のグラフは、非負値となるように修正された正常な状態の軸受１０６に対応する係数αを示す。図１２のグラフの横軸は、係数αのインデックスを示し、縦軸は、対応するデータの値を示す。

そして、発明者らは、軸受１０６が疵のない正常な状態である場合、図１３Ａ〜図１３Ｃに示す各状態にある場合それぞれと、について、以下の作業をｑ’回行った。図１３Ａ〜図１３Ｃに示す各状態とは、図１３Ａに示す外輪に幅３ｍｍ・深さ約０．１ｍｍの疵がある状態と、図１３Ｂに示す外輪に幅０．５ｍｍ・深さ０．１ｍｍ以下の疵がある状態と、図１３Ｃに示す外輪に予め定められた大きさ（径が１ｍｍ等）以下の微小な疵がある状態とのそれぞれである。ｑ’は、２００としたが、１００等の他の値としてもよい。各場合について、軸受１０６についての連続したＭ個の計測データを取得し、取得した連続したＭ個の計測データに基づいて、式１１を用いて行列Σ_ｓと行列Ｕ_ｓとを求め、式５と式１３を用いて修正された自己回帰モデルの係数αを求める作業をｑ’回行った。

ｑ’が大きいほど診断結果の信頼性が向上するため、ｑ’は、大きい方が望ましい。ｑ’の適正な値は、診断の対象、診断結果に要求される信頼性の程度、計測データの種類等に応じ、実験等の手段を駆使して定めればよい。本実験３では、ｑ’は、１０以上とするのが望ましく、５０以上とするのがより望ましく、１００以上とするのがさらに望ましい。

これにより、疵のない正常な状態の軸受１０６と、図１３Ａ〜図１３Ｃに示す各状態にある軸受１０６それぞれと、について、ｑ’個の修正された自己回帰モデルの係数αを取得した。

そして、軸受１０６が正常な場合、図１３Ａ〜図１３Ｃに示す各状態にある場合それぞれについて、以下の式３７で示すようなデータ行列Ｙ_ｄを生成した。

式３７の行列の各列のデータは、取得されたｑ’個の修正された自己回帰モデルの係数αそれぞれを示すデータである。即ち、α_ｉ，ｊは、ｑ’個の修正された自己回帰モデルの係数αのうちのｊ個目の係数αにおけるｉ番目の要素を示す。
本実験では、このデータ行列Ｙ_ｄを、非負値の複数の行列に分解する。データ行列Ｙ_ｄを非負値の複数の行列に分解するためには、行列Ｙの各列の係数αのパターンを変えずに、行列Ｙ自体の各成分を、非負値にするよう調整する必要がある。そこで、例えば、行列Ｙの０未満の成分のうち、最小のものの絶対値を、データ行列Ｙ_ｄの各成分に加算することで、データ行列Ｙ_ｄの各成分が非負値となるようにする。データ行列Ｙ_ｄの各列には、修正された自己回帰モデルにおける係数αのパターンを示す非負値のデータが格納されていることとなる。

以下の式３８に示すように、非負値の行列に修正したデータ行列Ｙ_ｄを、複数の非負値行列に分解する。より具体的には、データ行列Ｙ_ｄを、教師基底行列Ｆと、教師基底に対する重みを示す教師重み行列Ｇと、データ行列Ｙ_ｄから抽出される基底である抽出基底を示す抽出基底行列Ｑと、抽出基底に対する重みを示す抽出基底重み行列Ｗと、に分解する。教師基底行列Ｆと、教師重み行列Ｇと、抽出基底行列Ｑと、抽出基底重み行列Ｗは、何れも非負値行列である。本実験では、教師基底は、疵のない正常な状態の軸受１０６に対応する係数αのパターンを示す基底ベクトルの集合である。抽出基底は、疵のない正常な状態の軸受１０６に対応する係数αのパターン以外の係数αのパターンを示す基底ベクトルの集合である。疵のない正常な状態の軸受１０６に対応する係数αのパターン以外の係数αのパターンは、例えば、状態が不知の軸受１０６に対応する係数αのパターンを含む。

Ｙ_ｄをＦとＧとＱとＷとに分解する処理は、以下のような処理とみなすことができる。即ち、Ｙ_ｄの各列に格納された係数αのパターンのデータを、教師基底に属するパターンと、それ以外のパターンと、の重み付けの足し合わせとして分解する処理とみなすことができる。

Ｆは、ｍ×１のサイズの行列である。Ｇは、１×ｑ’のサイズの行列である。Ｑは、ｍ×（データ行列から抽出される基底ベクトルの数を示す数ｎ（本実験では、１とした））のサイズの行列である。Ｗは、ｎ×ｑ’のサイズの行列である。
Ｇにおける各列の要素の値は、データ行列Ｙ_ｄのその列に格納された係数αのパターンを示すデータに含まれる教師基底の成分の重みを示す。Ｗにおけるｎ行目の各列の要素の値は、データ行列Ｙ_ｄのその列に格納された係数αのパターンを示すデータに含まれるＱのその列に格納された抽出基底の成分の重みを示す。

以下の式３９のようにコスト関数Ｊ_ｄを定義する。

式３９の右辺第１項のＦｒは、フロベニウスノルムであることを示す添え字である。即ち、式３９の右辺第１項は、フロベニウスノルムである。式３９の右辺第１項は、ＦとＧとの積とＱとＷとの積との和と、データ行列Ｙと、の差分の大きさを示す。式３９の右辺第２項は、教師重み行列のスパース制約を示す項である。式３９の右辺第２項のδ_１は、予め定められた実数である。式３９の右辺第２項のｋ、ｔは、それぞれ行列の行、列を示すインデックスである。

式３９の右辺第３項は、抽出重み行列のスパース制約を示す項である。式３９の右辺第３項のδ_２は、予め定められた実数である。式３９の右辺第３項のｌ、ｔは、それぞれ行列の行、列を示すインデックスである。式３９の右辺第４項は、抽出基底の正規直交制約（即ち、抽出基底行列Ｑの各列に格納されているパターンの全体が近似的に正規直交系をなすとの制約）を示す項である。即ち、式３９の右辺第４項は、正規直交性の向上に関する項である。式３９の右辺第４項のε_１は、予め定められた実数である。式３９の右辺第４項のｉ、ｉ’は、それぞれ行列の行、列を示すインデックスである。式３９の右辺第４項のＩは、単位行列である。式３９の右辺第５項は、教師基底と抽出基底との間の直交制約（即ち、教師基底行列Ｆの各列に格納されている全てのパターンと抽出基底行列Ｑの各列に格納されている全てのパターンとが近似的に直交するとの制約）を示す項である。式３９の右辺第５項のε_２は、予め定められた実数である。式３９の右辺第４項のｉ、ｉ’は、それぞれ行列の行、列を示すインデックスである。

本実験では、教師重み行列のスパース制約、抽出重み行列のスパース制約、抽出基底の正規直交制約、教師基底と抽出基底との間の直交制約の何れも付加しない。即ち、δ_１、δ_２、ε_１、ε_２それぞれは、０とする。そのため、コスト関数Ｊ_ｄは、式３９の右辺第１項のみの式となる。

Ｊ_ｄの値を適正化するように、Ｇ、Ｑ、Ｗの値を更新することで、Ｙ_ｄを、ＦとＧとＱとＷとに分解する。以下の式４０で定義されるコスト関数Ｊ_ｄの最小化問題を解くこととなる。

本実験では、コスト関数Ｊ_ｄを最小化するために、以下のような処理を行う。即ち、ＧとＱとＷとの更新を繰り返すことで、コスト関数Ｊ_ｄを最小化する。また、コスト関数Ｊ_ｄを最小化するために、以下に記載する方法以外の方法（例えば、ニュートン法、確率的勾配降下法等の最適化手法）を用いてもよい。
Ｇの更新について説明する。本実験では、以下の式４１に示すようにＧの更新を行った。

式４１のｋ、ｔは、それぞれ行列の行、列を示すインデックスである。η^Ｇ _ｋｔは、Ｇのｋ、ｔ要素についての更新のステップサイズを示す非負値の値である。
コスト関数Ｊ_ｄをＧで偏微分することで、以下の式４２が得られる。

式４２のＥ_３は、サイズが１×ｑ’の全ての成分が１の行列である。
式４１に式４２を代入することで、更新後のＧを示す以下の式４３が得られる。

行列Ｇは、非負値であるため、式４３の値は、非負値である必要がある。式４３の最右辺の末尾の項は、必ず非負値となる項である。即ち、式４３の最右辺から末尾の項を除いた式の値が０であるならば、必ず更新後のＧは、非負値となる。そのため、以下の式４４に示すように、式４３の最右辺から末尾の項を除いた式が０であるという条件を満たすように、行列Ｇの更新を行うこととする。

式４４からη^Ｇ _ｋｔは、以下の式４５のようになる。

式４１に式４２と式４５とを代入することで、更新後のＧを示す以下の式４６が得られる。

即ち、Ｙ_ｄと、Ｆと、処理時点における更新前のＧ、Ｑ、Ｗと、予め定められた係数δ_１（本実験では、０）と、に基づいて、式４６を用いて、Ｇを更新した。
次に、Ｑの更新について説明する。本実験では、以下の式４７に示すようにＱを更新した。

式４７のｗ、ｔは、それぞれ行列の行、列を示すインデックスである。η^Ｑ _ｗｔは、Ｑのｗ、ｔ要素についての更新のステップサイズを示す非負値の値である。
コスト関数Ｊ_ｄをＱで偏微分することで、以下の式４８が得られる。

式４８のＥ_４は、サイズがｎ×ｎの全ての成分が１の行列である。式４８のＥ_５は、サイズが１×１の全ての成分が１の行列である。
式４７に式４８を代入することで、更新後のＱを示す以下の式４９が得られる。

行列Ｑは、非負値であるため、式４９の値は、非負値である必要がある。式４９の最右辺の末尾の項は、必ず非負値となる項である。即ち、式４９の最右辺から末尾の項を除いた式の値が０であるならば、必ず更新後のＱは、非負値となる。そのため、以下の式５０に示すように、式４９の最右辺から末尾の項を除いた式が０であるという条件を満たすように、行列Ｑを更新することとする。

式５０からη^Ｑ _ｗｔは、以下の式５１のようになる。

４７式に４８式と５１式とを代入することで、更新後のＧを示す以下の式５２が得られる。

即ち、Ｙ_ｄと、Ｆと、処理時点における更新前のＧ、Ｑ、Ｗと、予め定められた係数ε_１、ε_２（本実験では、共に０）と、に基づいて、式５２を用いて、Ｑを更新した。
次に、Ｗの更新について説明する。本実験では、以下の式５３に示すようにＷを更新した。

式５３のｌ、ｔは、それぞれ行列の行、列を示すインデックスである。η^Ｗ _ｌｔは、Ｗのｌ、ｔ要素についての更新のステップサイズを示す非負値の値である。
コスト関数Ｊ_ｄをＷで偏微分することで、以下の式５４が得られる。

式５４のＥ_６は、サイズが１×ｑ’の全ての成分が１の行列である。
式５３に式５４を代入することで、更新後のＷを示す以下の式５５が得られる。

行列Ｗは、非負値であるため、式５５の値は、非負値である必要がある。式５５の最右辺の末尾の項は、必ず非負値となる項である。即ち、式５５の最右辺から末尾の項を除いた式の値が０であるならば、必ず更新後のＷは、非負値となる。そのため、以下の式５６に示すように、式５５の最右辺から末尾の項を除いた式が０であるという条件を満たすように、行列Ｗの更新を行うこととする。

式５６からη^Ｗ _ｌｔは、以下の式５７のようになる。

５３式に５４式と５７式とを代入することで、更新後のＷを示す以下の式５８が得られる。

即ち、Ｙ_ｄと、Ｆと、処理時点における更新前のＧ、Ｑ、Ｗと、予め定められた係数δ_２（本実験では、０）と、に基づいて、式５８を用いて、Ｇを更新した。
以上のように、本実験では、式４６を用いたＧの更新と、式５２を用いたＱの更新と、式５８を用いたＷの更新と、を繰り返す。各更新段階において、更新後のＧ、Ｑ、Ｗから、それぞれ更新前のＧ、Ｑ、Ｗを引いた差であるΔＧ、ΔＱ、ΔＷについて、これらの行列のフロベニウスノルムを計算する。そして、各更新段階におけるΔＧ、ΔＱ、ΔＷのフロベニウスノルムの、それぞれ最初の更新段階におけるΔＧ、ΔＱ、ΔＷのフロベニウスノルムに対する比が、予め定められた閾値以下となる更新段階をもって、コスト関数Ｊ_ｄの値が収束したとした。そして、Ｊ_ｄが収束した際のＧ、Ｑ、Ｗを、コスト関数Ｊ_ｄを最小化させる行列とした。これにより、データ行列Ｙ_ｄを、式３８に示すように、それぞれ非負値の行列であるＦとＧとＱとＷとに分解した。

本実験では、以上のようにして、コスト関数Ｊ_ｄを最小化させるＧとＱとＷとを求めることとした。ただし、他の例として、式４６を用いたＧの更新と、式５２を用いたＱの更新と、式５８を用いたＷの更新と、を繰り返し行い、Ｇ、Ｑ、Ｗの更新前と更新後とでコスト関数Ｊ_ｄの値の変動が予め定められた閾値以下となることが、予め定められた回数連続したら、コスト関数Ｊ_ｄの値が収束したこととしてもよい。そして、Ｊ_ｄが収束した際のＧ、Ｑ、Ｗを、コスト関数Ｊ_ｄを最小化させる行列としてもよい。

図１４Ａ〜図１４Ｂを用いて、疵のない正常な状態にある軸受１０６についての本実験の実験結果について説明する。
図１４Ａのグラフは、Ｆが示す教師基底（に属する係数αのパターン）と、正常な状態にある軸受１０６についての計測データから本実験で求めたＱが示す抽出基底（に属する係数αのパターン）と、を示すグラフである。図１４Ａのグラフの横軸は、係数αのインデックスを示し、縦軸は、対応するデータの値を示す。

図１４Ｂのグラフは、正常な状態にある軸受１０６についての計測データから本実験で求めたＧとＷとを示すグラフである。図１４Ｂのグラフの横軸は、ＧとＷとの列のインデックスを示し、縦軸は、対応する列のＧとＷとの要素の値を示す。
ＧとＷとの各要素は、データ行列Ｙ_ｄに格納されている係数αのパターンを示すデータそれぞれの中に、教師基底の成分と抽出基底の成分とのそれぞれが、どの程度の重みで存在しているかを示す重みを示す。そのため、図１４Ｂのグラフは、データ行列Ｙ_ｄ内に、教師基底の成分と抽出基底の成分とのそれぞれが、どの程度の重みで存在しているかを示すグラフである。

図１４Ａのグラフから、教師基底と抽出基底とは、似通っていることが見て取れる。また、図１４Ｂのグラフから、データ行列Ｙ_ｄ内における教師基底の成分と抽出基底の成分との重みには、顕著な違いがみられないことが分かる。図１４Ｂのグラフにおいて、ＧとＷとを、２００個の列の要素それぞれについて比較すると、１１３点においてＧの方が大きくなっており、８７点においてＷの方が大きくなっている。

図１５Ａ〜図１５Ｂを用いて、図１３Ａに示すように幅３ｍｍ・深さ約０．１ｍｍの線上の疵が外輪に存在する状態にある軸受１０６についての本実験の実験結果について説明する。

図１５Ａのグラフは、Ｆが示す教師基底（に属する係数αのパターン）と、幅３ｍｍ・深さ約０．１ｍｍの線上の疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＱが示す抽出基底（に属する係数αのパターン）と、を示すグラフである。図１５Ａのグラフの横軸は、係数αのインデックスを示し、縦軸は、対応するデータの値を示す。

図１５Ｂのグラフは、幅３ｍｍ・深さ約０．１ｍｍの線上の疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＧとＷとを示すグラフである。図１５Ｂのグラフの横軸は、ＧとＷとの列のインデックスを示し、縦軸は、対応する列のＧとＷとの要素の値を示す。

図１５Ａのグラフから、教師基底と抽出基底とは、顕著に異なることが見て取れる。また、図１５Ｂのグラフから、データ行列Ｙ_ｄ内における教師基底の成分と抽出基底の成分との重みには、顕著な違いがみられることが分かる。抽出基底の成分の重みが、教師基底の成分の重みよりも顕著に高くなっている。図１５Ｂのグラフにおいて、ＧとＷとを、２００個の列の要素それぞれについて比較すると、２００個全てにおいてＷの方が大きくなる。

図１６Ａ〜図１６Ｂを用いて、図１３Ｂに示すように幅０．５ｍｍ・深さ約０．１ｍｍ以下の線上の疵が外輪に存在する状態にある軸受１０６についての本実験の実験結果について説明する。

図１６Ａのグラフは、Ｆが示す教師基底と、幅０．５ｍｍ・深さ０．１ｍｍ以下の線上の疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＱが示す抽出基底と、を示すグラフである。図１６Ａのグラフの横軸は、係数αのインデックスを示し、縦軸は、対応するデータの値を示す。

図１６Ｂのグラフは、幅０．５ｍｍ・深さ０．１ｍｍ以下の線上の疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＧとＷとを示すグラフである。図１６Ｂのグラフの横軸は、ＧとＷとの列のインデックスを示し、縦軸は、対応する列のＧとＷとの要素の値を示す。

図１６Ａのグラフから、教師基底と抽出基底とは、大まかな波形の動きとしては、似通っているともとれる。しかし、抽出基底に属する基本ベクトルの波形は、教師基底に属する基本ベクトルの波形に比べてより細かな振動が多数存在することが分かる。また、図１６Ｂのグラフから、データ行列Ｙ_ｄ内における教師基底の成分と抽出基底の成分との重みには、顕著な違いがみられることが分かる。図１６Ｂのグラフにおいて、ＧとＷとを、２００個の列の要素それぞれについて比較すると、１９９個においてＷの方が大きくなっており、１個においてＧの方が大きくなっている。

図１７Ａ〜図１７Ｂを用いて、図１３Ｃに示すように微小な疵が外輪に存在する状態にある軸受１０６についての本実験の実験結果について説明する。

図１７Ａのグラフは、Ｆが示す教師基底と、微小な疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＱが示す抽出基底と、を示すグラフである。図１７Ａのグラフの横軸は、係数αのインデックスを示し、縦軸は、対応するデータの値を示す。

図１７Ｂのグラフは、微小な疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＧとＷとを示すグラフである。図１７Ｂのグラフの横軸は、ＧとＷとの列のインデックスを示し、縦軸は、対応する列のＧとＷとの要素の値を示す。

図１７Ａのグラフから、教師基底と抽出基底とは、似通っているともとれる。また、図１７Ｂのグラフから、データ行列Ｙ_ｄ内における教師基底の成分と抽出基底の成分との重みには、一見、顕著な違いがないようにもとれる。しかし、図１７Ｂのグラフを見ると、Ｗの値は、大半が０．１２付近に存在しており、Ｇの値は、大半が０．０３付近に存在していることが分かる。図１７Ｂのグラフにおいて、ＧとＷとを、２００個の列の要素それぞれについて比較すると、１５９個においてＷの方が大きくなっており、４１個においてＧの方が大きくなっている。図１４Ｂのグラフと比べて、ＧがＷよりも大きくなっている箇所が顕著に多いのが分かる。

図１４Ａ〜図１７Ｂに示した実験結果から、データ行列Ｙ_ｄに対応する軸受１０６の状態と、教師基底に対応する軸受１０６の状態と、が異なっている場合、データ行列Ｙ_ｄの分解により得られたＷとＧとも顕著に異なる傾向が確認された。

発明者らは、本実験の実験結果から、データ行列Ｙ_ｄの分解により得られたＧとＷとを比較することで、データ行列Ｙ_ｄに対応する軸受１０６の状態が、教師基底が示す状態であるか否かを診断することができるとの知見を得た。

そこで、本実施形態の処理は、以下のような処理である。即ち、予め定められた状態の軸受１０６に対応する修正された自己回帰モデルの係数αのパターンを示す非負値のデータを、教師基底に属する基底ベクトルとして求める。そして、診断対象の軸受１０６から計測された計測データからデータ行列Ｙ_ｄを生成し、生成したデータ行列Ｙ_ｄを、教師基底を示す行列と、教師基底の重みを示す教師重み行列と、抽出基底を示す抽出基底行列と、抽出基底の重みを示す抽出重み行列と、に分解する。そして、データ行列Ｙ_ｄの分解により得られた教師重み行列と、抽出重み行列と、に基づいて、軸受１０６の状態が、教師基底に対応する状態であるか否かを診断する処理である。

（システム構成）
図１８は、本実施形態の診断システムのシステム構成の一例を示す図である。診断システムは、周期的な運動を行う物体についての状態診断を行うシステムである。診断システムは、情報処理装置１８００、振動計測装置１８０１を含む。診断システムは、鉄道台車に利用されている軸受の状態診断を行う。本実施形態では、診断システムは、図１Ａおよび図１Ｂと同様の状況で振動計測位置１１１に設置された振動計測装置１８０１により計測された振動の計測データに基づいて、軸受１０６の状態を診断する。

情報処理装置１８００は、振動計測装置１８０１により計測された信号に基づいて、軸受の状態診断を行うパーソナルコンピュータ（ＰＣ）、サーバ装置、タブレット装置等の情報処理装置である。また、情報処理装置１８００は、電車に組み込まれたコンピュータ等であってもよい。

振動計測装置１８０１は、加速度センサ等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を有線又は無線で情報処理装置１８００等の外部の装置に出力する計測装置である。

（情報処理装置の機能構成）
図１８を参照しながら、情報処理装置１８００が有する機能の一例を説明する。
情報処理装置１８００は、取得部１８１０、決定部１８２０、学習部１８３０、分解部１８４０、診断部１８５０、出力部１８６０を含む。

≪取得部１８１０≫
取得部１８１０は、周期的な運動を行う物体について、この周期的な運動に係る計測データｙを取得する。
本実施形態では、取得部１８１０は、軸受１０６の内輪を回転させ、振動計測位置１１１に設置された振動計測装置１８０１により計測された振動の計測データｙを取得する。

≪決定部１８２０≫
決定部１８２０は、取得部１８１０により取得された計測データｙに基づいて、修正された自己回帰モデルにおける係数αを決定する。修正された自己回帰モデルは、計測データｙの実績値ｙ_ｋ−１〜ｙ_ｋ−ｍと、この実績値に対する係数α（＝α_１〜α_ｍ）と、を用いて、計測データｙの予測値ｙ＾_ｋを表す式１である。決定部１８２０は、修正された自己回帰モデルにおける係数αを、式１２を用いて決定する。式１２は、一般的に知られている自己回帰モデルにおける係数の決定に用いる式６（ユール・ウォーカー方程式）において、式７で表される自己相関行列Ｒの代わりに、自己相関行列Ｒから物体の状態診断に有用な成分を抽出した第１の行列Ｒ’（式１０）を用いる式である。式６は、時差が０からｍ−１までの計測データｙの自己相関を成分とする自己相関行列Ｒを係数行列とし、時差が１からｍまでの計測データｙの自己相関を成分とする自己相関ベクトルを定数ベクトルとする方程式である。

式６は、式１で算出される計測データｙの予測値ｙ＾_ｋと、計測データｙの予測値ｙ＾_ｋに対応する時刻ｋにおける計測データｙの実測値ｙ_ｋとの二乗誤差を最小化する条件を示す条件式として導出することができる。第１の行列Ｒ’は、自己相関行列Ｒの固有値を対角成分とする対角行列Σと、自己相関行列Ｒの固有ベクトルを列成分とする直交行列Ｕと、から導出される。自己相関行列Ｒの固有値は、自己相関行列Ｒを特異値分解（式８）することで導出される。第１の行列Ｒ’は、自己相関行列Ｒの固有値のうち１以上且つｍ未満の設定された使用固有値数ｓ個の固有値を用いて、対角行列Σの部分行列であって、使用固有値数ｓ個の固有値を対角成分とする行列Σ_ｓと、直交行列Ｕの部分行列であって、使用固有値数ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列Ｕ_ｓと、から導出される行列Ｕ_ｓΣ_ｓＵ_ｓ ^Ｔである。使用固有値数ｓ個の固有値は、自己相関行列Ｒの固有値のうち、値が最大の固有値を含むようにすればよいが、値が大きいものから順に選ばれるのが好ましい。

≪学習部１８３０≫
学習部１８３０は、複数の状態の物体について、予め取得部１８１０により取得された計測データから決定部１８２０により決定された修正された自己回帰モデルにおける係数αを学習データとして、予め定められた状態の物体に対応する修正された自己回帰モデルにおける係数αのパターンを学習により、教師基底に属する基底ベクトルとして特定する。

≪分解部１８４０≫
分解部１８４０は、診断対象の物体から計測された計測データに基づいて決定部１８２０により決定された係数αを予め定められた個数だけ取得する。分解部１８４０は、取得した係数αそれぞれのパターンを示す非負値のデータを格納するデータ行列を生成する。そして、分解部１８４０は、生成したデータ行列を、学習部１８３０により学習された教師基底を格納する教師基底行列と、データ行列内の教師基底の成分の重みを示す教師重み行列と、データ行列から抽出される基底が格納される抽出基底行列と、データ行列内の抽出基底の成分の重みを示す抽出重み行列と、に分解する。

≪診断部１８５０≫
診断部１８５０は、分解部１８４０により求められた教師重み行列と抽出重み行列とに基づいて、物体の状態が教師基底に対応する状態であるか否かを診断する。
本実施形態では、軸受１０６の状態が教師基底に対応する状態であるか否かを診断の結果とする。
≪出力部１８６０≫
出力部１８６０は、診断部１８５０による診断の結果に係る情報を出力する。

（情報処理装置のハードウェア構成）
図１９は、情報処理装置１８００のハードウェア構成の一例を示す図である。
情報処理装置１８００は、ＣＰＵ１９００、主記憶装置１９０１、補助記憶装置１９０２、入出力Ｉ／Ｆ１９０３を含む。各構成要素は、システムバス１９０４を介して、相互に通信可能に接続されている。

ＣＰＵ１９００は、情報処理装置１８００を制御する中央演算装置である。主記憶装置１９０１は、ＣＰＵ１９００のワークエリアやデータの一時的な保管場所として機能するＲＡＭ（ＲａｎｄｏｍＡｃｃｅｓｓＭｅｍｏｒｙ）等の記憶装置である。補助記憶装置１９０２は、各種のプログラム、設定データ、振動計測装置１８０１から出力される計測データ、診断情報等を記憶するＲＯＭ（ＲｅａｄＯｎｌｙＭｅｍｏｒｙ）、ＨＤＤ（ＨａｒｄＤｉｓｋＤｒｉｖｅ）、ＳＳＤ（ＳｏｌｉｄＳｔａｔｅＤｒｉｖｅ）等の記憶装置である。入出力Ｉ／Ｆ１９０３は、振動計測装置１８０１等の外部の装置との間での情報のやり取りに利用されるインターフェースである。
ＣＰＵ１９００が、補助記憶装置１９０２等に記憶されたプログラムに基づき処理を実行することによって、図１８で説明した情報処理装置１８００の機能及び図２０、２１で後述するフローチャートの処理等が実現される。

（学習処理）
図２０は、診断システムが実行する学習処理の一例を示すフローチャートである。図２０を用いて、教師基底を学習により特定する処理を説明する。
本実施形態では、情報処理装置１８００は、軸受１０６が外輪に疵のある状態である場合、軸受１０６が転動体に疵のある状態である場合、軸受１０６が内輪に疵のある状態である場合、軸受１０６に疵がなく正常な状態である場合、軸受１０６が外輪に疵のある状態である場合のそれぞれについて、予め、ｑ（５０）回、以下の処理を繰り返す。

即ち、取得部１８１０は、振動計測装置１８０１により計測された振動の計測データｙを取得する。取得部１８１０は、計測データｙとして、ランダムに設定した時刻を時刻１とした場合の時刻１〜Ｍまでの計測データであるｙ_１〜ｙ_Ｍを取得する。そして、決定部１８２０は、取得した計測データｙと、予め設定された定数Ｍと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数ｍと、に基づいて、式５と式７とを用いて自己相関行列Ｒを生成する。決定部１８２０は、予め記憶されているＭ、ｍの情報を読み出すことで、Ｍとｍとを取得する。本実施形態では、ｍの値は、５００である。また、Ｍは、ｍよりも大きい整数である。

そして、決定部１８２０は、生成した自己相関行列Ｒを特異値分解することで、式８の直交行列Ｕと対角行列Σを取得し、対角行列Σから自己相関行列Ｒの固有値σ_１１〜σ_ｍｍを取得する。決定部１８２０は、自己相関行列Ｒの複数の固有値であるσ_１１〜σ_ｍｍのうち、最大のものから使用固有値数ｓだけの固有値σ_１１〜σ_ｓｓを、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値として選択する。本実施形態では、使用固有値数ｓは、３とするが、２以下としてもよいし、４以上としてもよい。そして、決定部１８２０は、計測データｙと、固有値σ_１１〜σ_ｓｓと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、式１３を用いて、修正された自己回帰モデルの係数αを決定する。そして、決定部１８２０は、決定した係数αを、非負値化して、学習データとして補助記憶装置１９０２に記憶する。非負値化は、例えば、係数αのそれぞれの成分に、係数αに含まれる０未満の最低の成分の値に−１をかけた値を足すことにより行われる。
これにより、情報処理装置１８００は、軸受１０６が外輪に疵のある状態である場合、軸受１０６が転動体に疵のある状態である場合、軸受１０６が内輪に疵のある状態である場合、軸受１０６に疵がなく正常な状態である場合、軸受１０６が外輪に疵のある状態である場合のそれぞれについて、ｑ（５０）個ずつ修正された自己回帰モデルにおける係数αのデータを、学習データとして取得し、記憶することができる。

本実施形態では、軸受１０６が取り得る状態は、疵のない正常な状態、内輪に疵がある状態、保持器に疵がある状態、転動体に疵がある状態、および外輪に疵がある状態の５つであるとする。即ち、軸受１０６が取り得る状態の数ｐは、５となる。

Ｓ２００１において、学習部１８３０は、補助記憶装置１９０２に記憶されている学習データを取得する。

Ｓ２００２において、学習部１８３０は、Ｓ２００１で取得した学習データに基づいて、式１９の行列Ｙを生成する。学習部１８３０は、ｍ（５００）×ｄ（２５０）のサイズの行列を生成し、各列のデータを、学習データに含まれる各係数αの値とすることで、行列Ｙを生成する。また、ｍ×ｐのサイズの行列Ａを生成し、行列Ａの各成分に、予め定められた初期値を設定する。また、ｐ×ｄのサイズ行列Ｐを生成し、行列Ｐの各成分に、予め定められた初期値を設定する。

Ｓ２００３において、学習部１８３０は、Ｓ２００２で生成した行列Ｙと、行列Ａと、行列Ｐと、予め定められた係数βと、予め定められた係数γと、に基づいて、式２１を用いて、コスト関数Ｊの値を、コストとして求める。学習部１８３０は、求めたコストを、更新前コストとして、主記憶装置１９０１等に記憶する。

Ｓ２００４において、学習部１８３０は、行列Ｙと、行列Ａと、行列Ｐと、係数γと、サイズがｐ×ｐの全ての成分が１の行列である行列Ｅ_２と、に基づいて、式３４を用いて、行列Ａを更新する。また、学習部１８３０は、更新した行列Ａを、式３５を用いて、更に更新する。また、学習部１８３０は、行列Ｙと、行列Ａと、行列Ｐと、係数βと、サイズがｍ×ｐの全ての成分が１の行列である行列Ｅ_１と、に基づいて、式３０を用いて、行列Ｐを更新する。

Ｓ２００５において、学習部１８３０は、行列Ｙと、Ｓ２００３で更新した行列Ａと、Ｓ２００３で更新した行列Ｐと、係数βと、係数γと、に基づいて、式２１を用いて、コスト関数Ｊの値を、コストとして求める。学習部１８３０は、求めたコストを、更新後コストとして、主記憶装置１９０１等に記憶する。

Ｓ２００６において、学習部１８３０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であるか否かを判定する。学習部１８３０は、更新後コストと、更新前コストと、の差分が、予め定められた閾値以下であると判定した場合、Ｓ２００７の処理に進む。また、学習部１８３０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値よりも大きいと判定した場合、Ｓ２００８の処理に進む。学習部１８３０は、Ｓ２００６の処理の後、更新後コストを、新たな更新前コストとして、主記憶装置１９０１に記憶する。

Ｓ２００７において、学習部１８３０は、主記憶装置１９０１等に記憶されたカウンタの値に１を加える。学習部１８３０は、このカウンタの値を、図２０の処理の開始前に０に初期化している。
Ｓ２００８において、学習部１８３０は、主記憶装置１９０１等に記憶されたカウンタの値を、０に初期化して、Ｓ２００４の処理に進む。

Ｓ２００９において、学習部１８３０は、主記憶装置１９０１等に記憶されたカウンタの値が、予め定められた閾値以上であるか否かを判定する。学習部１８３０は、主記憶装置１９０１等に記憶されたカウンタの値が、予め定められた閾値以上であると判定した場合、コスト関数Ｊの値が極小値に収束したとして、Ｓ２０１０の処理に進む。学習部１８３０は、主記憶装置１９０１等に記憶されたカウンタの値が、予め定められた閾値未満であると判定した場合、Ｓ２００４の処理に進む。

Ｓ２０１０において、学習部１８３０は、行列Ｐに基づいて、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。

行列Ｙのある列のデータは、行列Ａと、行列Ｐにおけるその列と、をかけたデータとして求まる。即ち、行列Ｙのある列のデータは、行列Ａの各列のデータを、行列Ｐのその列における各行の値を重みとして、重ね合せたデータとなる。そこで、学習部１８３０は、Ｓ２０１０で以下のようにして、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。

学習部１８３０は、行列Ｙの各列のうち、ある１つの状態である軸受１０６に対応する係数αを示すデータを格納する列を特定する。例えば、学習部１８３０は、行列Ｙの各列のうち、外輪に疵がある状態の軸受１０６に対応する係数αを示すデータを格納する列である１列目〜ｑ列目を特定する。そして、学習部１８３０は、行列Ｐにおける特定した列（１列目〜ｑ列目）について、行毎に値を集計する。そして、学習部１８３０は、集計した値のうち、最も大きいものに対応する行を特定する。学習部１８３０は、特定した行列Ｐの行とかけ合わされる行列Ａの列を特定する。例えば、学習部１８３０は、行列Ｐの１行目を特定した場合、行列Ｐの１行目とかけ合わされる行列Ａの列は、１列目なので、１列目を特定する。そして、学習部１８３０は、行列Ａの特定した列が示すデータを、その状態である軸受１０６に対応する係数αの代表的なパターンを示すデータとして特定する。

学習部１８３０は、軸受１０６が取り得る各状態について、以上の処理を行うことで、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。

学習部１８３０は、特定した各状態の軸受１０６に対応する修正された自己回帰モデルにおける係数αの代表的なパターンを示す非負値のデータのうち、予め定められた状態の軸受１０６に対応するデータを、教師基底に属する基底ベクトルとして、補助記憶装置１９０２等に記憶する。本実施形態では、この予め定められた状態は、疵のない正常な状態であるとする。そのため、学習部１８３０は、正常な状態に対応する係数αのパターンを示すデータを、教師基底に属する基底ベクトルとして、補助記憶装置１９０２等に記憶する。

（診断処理）
本実施形態では、診断システムは、軸受１０６の状態が教師基底に対応する状態か否かを特定することで、軸受１０６の状態を診断することとする。ただし、他の例として、診断システムは、予め振動計測装置１８０１により計測された計測データに基づいて、軸受１０６の状態を診断することとしてもよいし、リアルタイムで振動計測装置１８０１により計測されて出力され続けている計測データを用いて、稼働中の軸受１０６の状態を診断することとしてもよい。

図２１は、診断処理の一例を示すフローチャートである。
Ｓ２１０１において、取得部１８１０は、振動計測装置１８０１により計測された振動の計測データｙを取得する。取得部１８１０は、計測データｙとして、振動の計測期間内からランダムに時刻を選択する。そして、取得部１８１０は、選択した時刻を時刻１として、時刻１〜Ｍまでの計測データを、ｙ_１〜ｙ_Ｍとして取得する。

Ｓ２１０２において、決定部１８２０は、Ｓ２１０１で取得した計測データｙと、予め設定された定数Ｍと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数ｍと、に基づいて、式５と式７とを用いて自己相関行列Ｒを生成する。決定部１８２０は、予め記憶されているＭ、ｍの情報を読み出すことで、Ｍとｍとを取得する。本実施形態では、ｍの値は、５００である。また、Ｍは、ｍよりも大きい整数である。

Ｓ２１０３において、決定部１８２０は、Ｓ２１０２で生成した自己相関行列Ｒを特異値分解することで、式８の直交行列Ｕと対角行列Σを取得し、対角行列Σから自己相関行列Ｒの固有値σ_１１〜σ_ｍｍを取得する。

Ｓ２１０４において、決定部１８２０は、自己相関行列Ｒの複数の固有値であるσ_１１〜σ_ｍｍのうち、最大のものから使用固有値数ｓだけの固有値σ_１１〜σ_ｓｓを、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値として選択する。本実施形態では、使用固有値数ｓは、３とするが、２以下としてもよいし、４以上としてもよい。そして、決定部１８２０は、計測データｙと、固有値σ_１１〜σ_ｓｓと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、式１３を用いて、修正された自己回帰モデルの係数αを決定する。Ｓ２１０４で決定された係数αは、計測データから決定された修正された自己回帰モデルにおける係数である修正係数の一例である。Ｓ２１０４で決定された係数αのパターンは、計測パターンの一例である。

Ｓ２１０５において、決定部１８２０は、予め定められた数ｑ’個の係数αを決定したか否かを判定する。決定部１８２０は、予め定められた数ｑ’個の係数αを決定したと判定した場合、処理をＳ２１０５に進め、予め定められた数ｑ’個以下の係数αしか決定していないと判定した場合、処理をＳ２１０１に進め、係数αの決定処理を繰り返す。

Ｓ２１０６において、分解部１８４０は、補助記憶装置１９０２から、Ｓ２０１０で記憶された教師基底を取得する。

Ｓ２１０７において、分解部１８４０は、決定部１８２０により決定されたｑ’個の係数αを用いて、式３７のように、ｍ×ｑ’のサイズのデータ行列Ｙ_ｄを生成する。即ち、データ行列Ｙ_ｄの各列が、Ｓ２１０４で決定された係数αのパターンのデータを格納することとなる。また、分解部１８４０は、データ行列Ｙ_ｄの各成分を、非負値化する。非負値化は、例えば、０未満の行列Ｙの成分のうち最低のものに−１をかけた値を、行列Ｙの各成分に足すことにより行われる。

Ｓ２１０８において、分解部１８４０は、教師基底を格納するｍ×１のサイズの教師基底行列Ｆを生成する。そして、分解部１８４０は、Ｆの各要素に、教師基底に属するデータを設定する。また、分解部１８４０は、データ行列Ｙ_ｄに含まれる教師基底の成分の重みを格納する１×ｑ’のサイズの教師重み行列Ｇを生成する。また、分解部１８４０は、データ行列Ｙ_ｄから抽出されるパターンの集合である抽出基底を格納するｍ×予め定められた数ｎのサイズの行列Ｑを生成する。本実施形態では、ｎは１とする。また、分解部１８４０は、データ行列Ｙ_ｄに含まれる抽出基底の成分の重みを格納するｎ×ｑ’のサイズの教師重み行列Ｗを生成する。
そして、分解部１８４０は、Ｇ、Ｑ、Ｗそれぞれの各成分に予め定められた初期値を設定することで、Ｇ、Ｑ、Ｗそれぞれを初期化する。

Ｓ２１０９において、分解部１８４０は、Ｙ_ｄとＦとＧとＱとＷと予め定められた係数δ_１、δ_２、ε_１、ε_２と、に基づいて、式３９を用いて、コスト関数Ｊ_ｄの値を、コストとして求める。本実施形態では、δ_１、δ_２、ε_１、ε_２それぞれの値は、０とする。分解部１８４０は、求めたコストを、更新前コストとして、主記憶装置１９０１等に記憶する。

Ｓ２１１０において、分解部１８４０は、Ｇ、Ｑ、Ｗそれぞれを更新する。より具体的には、分解部１８４０は、Ｙ_ｄと、Ｆと、更新前のＧ、Ｑ、Ｗと、予め定められた係数δ_１（本実施形態では、０）と、サイズが１×ｑ’の全ての成分が１の行列であるＥ_３と、に基づいて、式４６を用いて、Ｇを更新する。ただし、分解部１８４０は、本実施形態ではδ_１が０なので、Ｅ_３を用いずに、Ｇを更新してもよい。

また、分解部１８４０は、Ｙ_ｄと、Ｆと、更新前のＧ、Ｑ、Ｗと、予め定められた係数ε_１、ε_２（本実施形態では、共に０）と、サイズがｎ×ｎの全ての成分が１の行列であるＥ_４と、サイズが１×１の全ての成分が１の行列であるＥ_５と、に基づいて、式５２を用いて、Ｑを更新する。ただし、分解部１８４０は、本実施形態ではε_１、ε_２それぞれが０なので、Ｅ_４とＥ_５とを用いずに、Ｑを更新してもよい。

また、分解部１８４０は、Ｙ_ｄと、Ｆと、更新前のＧ、Ｑ、Ｗと、予め定められた係数δ_２（本実施形態では、０）と、サイズが１×ｑ’の全ての成分が１の行列であるＥ_６と、に基づいて、式５８を用いて、Ｗを更新する。ただし、分解部１８４０は、本実施形態ではδ_２が０なので、Ｅ_６を用いずに、Ｑを更新してもよい。

Ｓ２１１１において、分解部１８４０は、Ｙ_ｄと、Ｆと、Ｓ２１１０で更新したＧ、Ｑ、Ｗと、係数δ_１、δ_２、ε_１、ε_２と、に基づいて、式３９を用いて、コスト関数Ｊ_ｄの値を、コストとして求める。分解部１８４０は、求めたコストを、更新後コストとして、主記憶装置１９０１等に記憶する。

Ｓ２１１２において、分解部１８４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であるか否かを判定する。分解部１８４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であると判定した場合、処理をＳ２１１３に進める。また、分解部１８４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値よりも大きいと判定した場合、処理をＳ２１１４に進める。分解部１８４０は、Ｓ２１１２の処理の後、更新後コストを、新たな更新前コストとして、主記憶装置１９０１に記憶する。

Ｓ２１１３において、分解部１８４０は、主記憶装置１９０１等に記憶されたカウンタの値に１を加える。なお、分解部１８４０は、このカウンタの値を、図２１の処理の開始前に０に初期化している。

Ｓ２１１４において、分解部１８４０は、主記憶装置１９０１等に記憶されたカウンタの値を、０に初期化して、処理をＳ２１１０に進める。

Ｓ２１１５において、分解部１８４０は、主記憶装置１９０１等に記憶されたカウンタの値が、予め定められた閾値以上であるか否かを判定する。分解部１８４０は、主記憶装置１９０１等に記憶されたカウンタの値が、予め定められた閾値以上であると判定した場合、コスト関数Ｊ_ｄの値が極小値に収束したとして、処理をＳ２１１６に進める。分解部１８４０は、カウンタの値が、予め定められた閾値未満であると判定した場合、処理をＳ２１１０に進める。

Ｓ２１１６において、診断部１８５０は、処理時点までに更新されたＧとＷとに基づいて、診断対象の軸受１０６の状態が教師基底に対応する状態であるか否かを診断する。診断部１８５０は、ＧとＷとの各列について、要素同士を比較して、何れが大きいかを特定する。そして、診断部１８５０は、Ｇの要素の方がＷの要素よりも大きい列の数が、予め定められた閾値以上である場合、診断対象の軸受１０６の状態が、教師基底に対応する状態（本実施形態では、正常な状態）であると診断する。予め定められた閾値は、例えば、ＧとＷとの列数の７５％の個数（本実施形態では１５０）である。また、診断部１８５０は、Ｇの要素の方がＷの要素よりも大きい列の数が、予め定められた閾値未満である場合、診断対象の軸受１０６の状態が、教師基底に対応する状態ではないと診断する。

Ｓ２１１７において、出力部１８６０は、Ｓ２１１６での軸受１０６の診断の結果を示す情報を出力する。出力部１８６０は、例えば、表示装置に、軸受１０６の診断の結果を示す情報を表示することで出力する。また、他の例としては、出力部１８６０は、例えば、軸受１０６の診断の結果を示す情報を、補助記憶装置１９０２に記憶することで出力することとしてもよい。また、出力部１８６０は、例えば、外部のサーバ装置等の設定された送信先に、軸受１０６の診断の結果を示す情報を送信することで出力することとしてもよい。

（まとめ）
診断対象の軸受１０６の計測データから得られた係数αのパターンに含まれる、予め定められた状態に対応するパターンの成分の重みと、それ以外のパターンの成分の重みと、に基づいて、軸受１０６の状態が予め定められた状態か否かを適切に診断できるとの知見が実験３から見出された。本実施形態では、診断システムは、この知見に基づいた軸受１０６の診断処理を行う。

即ち、診断システムは、予め、疵のない正常な状態の軸受１０６に対応する修正された自己回帰モデルの係数αのパターンを、教師基底に属する基底ベクトルとして求める。そして、診断システムは、診断対象の軸受１０６から計測された計測データからｑ’個の修正された自己回帰モデルの係数αを求め、求めたｑ’個の係数αそれぞれのパターンを示すデータを格納するデータ行列Ｙ_ｄを生成する。そして、診断システムは、Ｙ_ｄを、Ｆと、教師基底の重みを示す教師重み行列Ｇと、抽出基底を示す抽出基底行列Ｑと、抽出基底の重みを示す抽出重み行列Ｗと、に分解する。診断システムは、Ｙ_ｄを分解して得られたＧとＷとに基づいて、診断対象の軸受１０６が教師基底に対応する状態であるか否かを診断する。これにより、診断システムは、軸受１０６が予め定められた状態であるか否かを適切に診断できる。

また、診断システムは、自己相関行列Ｒの固有値のうち一部の固有値を用いることで、軸受１０６の診断に有用な成分がより多く残り、有用でない成分がより残らないように、計測データを近似する修正された自己回帰モデルの係数αを決定する。診断システムは、このような係数αを用いて、軸受１０６の診断を行うことで、様々な周波数域の信号成分を含んだ信号を用いて診断を行う場合に比べて、より精度よく軸受１０６の診断を行うことができる。

（変形例１）
本実施形態では、診断システムは、鉄道台車の軸受１０６の状態診断を行うこととした。しかし、診断システムは、歯車等の回転体、振動子等の振動体、バネ等の伸縮体、生物の心臓等の周期的な運動を行う他の物体についての状態を診断することとしてもよい。

（変形例２）
本実施形態では、ｍは、予め設定された５００という数であるとした。しかし、診断対象に応じて実験を行う等して適切に計測データを近似できる過去のデータの数を特定し、特定した数をｍの値としてもよい。例えば、情報処理装置１８００は、診断対象から計測された計測データについて、ある時刻の計測データを、その時刻よりも過去の計測データを複数用いて自己回帰モデル等を利用して近似して、近似した値とその時刻の計測データとの差分が設定された閾値未満となる場合の近似に用いられた過去の計測データの数を、ｍとして決定してもよい。

（変形例３）
本実施形態では、情報処理装置１８００は、使用固有値数ｓを３とした。しかし、情報処理装置１８００は、以下のようにして、使用固有値数ｓの値を決定してもよい。即ち、情報処理装置１８００は、まず、自己相関行列Ｒの固有値全ての合計値を求める。そして、情報処理装置１８００は、自己相関行列Ｒの固有値のうち、最大のものからａ個の合計が、求めた合計値の設定された割合（例えば９０％）以上となるａのうち、最小のものを使用固有値数ｓとして決定してもよい。

（変形例４）
本実施形態では、情報処理装置１８００は、自己相関行列Ｒの固有値のうち、最大のものから、使用固有値数ｓ個の固有値を用いて、修正された自己回帰モデルにおける係数αを決定することとした。しかし、情報処理装置１８００は、自己相関行列Ｒの固有値のうち、任意のｓ個の固有値を用いて、修正された自己回帰モデルにおける係数αを決定することとしてもよい。例えば、情報処理装置１８００は自己相関行列Ｒの固有値のうち、最大のものと、３つ目に大きいものと、の２つを利用して、修正された自己回帰モデルにおける係数αを決定することとしてもよい。

（変形例５）
本実施形態では、式１９に示す行列Ｙの各列には、複数の状態の軸受１０６に対応する計測データから決定された修正された自己回帰モデルにおける係数αが、状態毎に、ｑ個ずつ格納されることとした。しかし、行列Ｙの各列には、複数の状態の軸受１０６に対応する計測データから決定された修正された自己回帰モデルにおける係数αが、状態毎に、異なる数ずつ格納されることとしてもよい。

本実施形態では、診断システムは、行列Ｙを、同じ状態の軸受１０６に対応する係数αを格納する列が連続するように、生成した。しかし、診断システムは、同じ状態の軸受１０６に対応する係数αを格納する列が連続しないように、行列Ｙを生成してもよい。

本実施形態では、診断システムは、各列に各状態の軸受１０６に対応する係数αが格納された行列Ｙを生成し、行列Ｙを非負値行列因子分解により学習することで、各状態の軸受１０６に対応する係数αの代表的なパターンを特定した。しかし、診断システムは、例えば、学習データを取得する際の処理と同様の処理で、ある状態の軸受１０６に対応する係数αを複数決定し、決定した複数の係数αの平均を求めることにより学習することで、求めた平均の係数が示すパターンを、その状態の軸受１０６に対応する係数αの代表的なパターンとして特定してもよい。

（変形例６）
本実施形態では、診断部１８５０は、Ｙ_ｄの分解により求まったＧとＷとに基づいて、図２１のＳ２１１６で説明した方法で軸受１０６の診断を行うこととした。しかし、診断部１８５０は、他の方法で、診断対象の軸受１０６の状態が教師基底に対応する状態であるか否かを診断してもよい。

例えば、診断部１８５０は、以下のようにしてもよい。即ち、診断部１８５０は、Ｇの各要素について、偏差を求めて、偏差が小さいものから予め定められた数を特定し、特定したＧの要素の平均値を求める。予め定められた数は、例えば、Ｇの全要素の７０％の数（本実施形態では、１４０）である。また、診断部１８５０は、Ｗの各要素について、偏差を求めて、偏差が小さいものから予め定められた数を特定し、特定したＷの要素の平均値を求める。予め定められた数は、例えば、Ｇの全要素の７０％の数（本実施形態では、１４０）である。そして、診断部１８５０は、求めたＧの要素の平均値がＷの要素の平均値よりも予め定められた閾値以上大きい場合、診断対象の軸受１０６が教師基底に対応する状態であると診断する。一方、診断部１８５０は、求めたＧの要素の平均値がＷの要素の平均値よりも予め定められた閾値以上大きくない場合、診断対象の軸受１０６が教師基底に対応する状態でないと診断する。

また、診断部１８５０は、以下のようにしてもよい。即ち、診断部１８５０は、Ｇの各要素のうち値が中間の予め定められた数の要素の平均値を求める。予め定められた数は、例えば、Ｇの全要素の７０％の数（本実施形態では、１４０）である。また、診断部１８５０は、Ｗの各要素のうち値が中間の予め定められた数の要素の平均値を求める。そして、診断部１８５０は、求めたＧの要素の平均値がＷの要素の平均値よりも予め定められた閾値以上大きい場合、診断対象の軸受１０６が教師基底に対応する状態であると診断する。一方、診断部１８５０は、求めたＧの要素の平均値がＷの要素の平均値よりも予め定められた閾値以上大きくない場合、診断対象の軸受１０６が教師基底に対応する状態でないと診断する。

（変形例７）
本実施形態では、データ行列Ｙ_ｄから抽出される、抽出基底に属する基底ベクトルの数（抽出基底行列Ｑの列数、抽出重み行列Ｗの行数）を１とした。ただし、データ行列Ｙ_ｄから抽出される、抽出基底に属する基底ベクトルの数を２以上としてもよい。その場合、Ｓ２１１６では、診断部１８５０は、以下のようにしてもよい。

即ち、診断部１８５０は、診断部１８５０は、Ｇと、Ｗの各行それぞれと、の各列について、要素同士を比較して、何れが大きいかを特定する。そして、診断部１８５０は、Ｇの要素が最も大きい列の数が、予め定められた閾値以上である場合、診断対象の軸受１０６の状態が、教師基底に対応する状態であると診断する。一方、診断部１８５０は、Ｇの要素が最も大きい列の数が、予め定められた閾値未満である場合、診断対象の軸受１０６の状態が、教師基底に対応する状態ではないと診断する。

（変形例８）
本実施形態では、診断システムは、各列に係数αのパターンを示す非負値のデータを格納するように行列Ｙ_ｄを生成した。ただし、他の例として、診断システムは、各列ではなく各行に係数αのパターンを示す非負値のデータを格納するように行列Ｙ_ｄを生成してもよい。その場合、診断システムは、行列Ｙ_ｄを、Ｙ_ｄ＝Ｇ^ＴＦ^Ｔ＋Ｗ^ＴＱ^Ｔのように、分解すればよい。その場合、診断システムは、コスト関数Ｊ_ｄとして、式３９における「ＦＧ＋ＱＷ−Ｙ_ｄ」の部分を「Ｇ^ＴＦ^Ｔ＋Ｗ^ＴＱ^Ｔ−Ｙ_ｄ」に置き換えた関数を用いればよい。

（変形例９）
本実施形態では、診断システムは、診断対象の軸受１０６が疵のない正常な状態であるか否かを診断することとした。ただし、診断システムは、診断対象の軸受１０６が他の状態（例えば、外輪に特定の疵の有る状態、内輪に疵の有る状態、転動体に疵の有る状態等）であるか否かを診断することとしてもよい。

（変形例１０）
本実施形態では、診断システムは、式３９で表される関数Ｊ_ｄを、コスト関数として用いることとした。ただし、診断システムは、Ｇ、Ｑ、Ｗの更新の指標となる関数であれば、他の関数をコスト関数として用いてもよい。例えば、診断システムは、式３９で表される関数Ｊ_ｄの逆数で表される関数をコスト関数として用いてもよい。その場合、診断システムは、コスト関数の値を適正化（最大化）するため、Ｇ、Ｑ、Ｗの値をコスト関数の値が増大するように更新していくこととなる。そして、診断システムは、コスト関数が収束したら、コスト関数の値が適正化（最大化）できたとして、その際のＧ、Ｑ、Ｗを、Ｙ_ｄの分解結果として取得する。

ここで、最小化とは、コスト関数Ｊ_ｄの値を最小化する解を導出する最適化アルゴリズムにおいてコスト関数Ｊ_ｄを最小化することを指す。即ち、最小化とは、コスト関数の値の厳密な最小値を求めるものに限られない。このことは、最大化についても同じである。

（変形例１１）
本実施形態では、診断システムは、図２０の学習処理において、行列Ｐと行列Ａとの更新前後で、コスト関数Ｊの値の変動が予め定められた閾値以下となることが、予め定められた回数連続するまで、行列Ｐと行列Ａとの更新を繰り返した。これにより、診断システムは、コスト関数Ｊが最小化した際の行列Ｐと行列Ａとを求めることとした。ただし、他の例として、診断システムは、以下のようにしてもよい。

即ち、診断システムは、式３０を用いて行列Ｐを更新し、式３４、式３５を用いて行列Ａを更新することを繰り返す。そして、診断システムは、各更新段階において、更新後の行列Ｐ及び行列Ａから、それぞれ更新前の行列Ｐ及び行列Ａを引いた差である行列△Ｐ及び行列△Ａを求める。診断システムは、行列△Ｐ及び行列△Ａそれぞれのフロベニウスノルムを計算する。診断システムは、各更新段階における行列△Ｐ及び行列△Ａのフロベニウスノルムの、それぞれ最初の更新段階における行列△Ｐ及び行列△Ａのフロベニウスノルムに対する比が、予め定められた閾値以下となる更新段階をもって、コスト関数Ｊの値が収束したとする。そして、診断システムは、その際の行列Ａと行列Ｐとをコスト関数Ｊを最小化させる行列として求める。

（変形例１２）
本実施形態では、診断システムは、図２１の診断処理において、Ｇ、Ｑ、Ｗの更新前と更新後とでコスト関数Ｊ_ｄの値の変動が予め定められた閾値以下となることが、予め定められた回数連続するまで、Ｇ、Ｑ、Ｗそれぞれの更新を繰り返すことで、コスト関数Ｊ_ｄの値が収束した際のＧ、Ｑ、Ｗを求めることとした。ただし、他の例として、診断システムは、以下のようにしてもよい。

即ち、診断システムは、式４６を用いたＧの更新と、式５２を用いたＱの更新と、式５８を用いたＷの更新と、を繰り返し行う。診断システムは、各更新段階において、更新後のＧ、Ｑ、Ｗから、それぞれ更新前のＧ、Ｑ、Ｗを引いた差であるΔＧ、ΔＱ、ΔＷについて、これらの行列のフロベニウスノルムを計算する。診断システムは、各更新段階におけるΔＧ、ΔＱ、ΔＷのフロベニウスノルムの、それぞれ最初の更新段階におけるΔＧ、ΔＱ、ΔＷのフロベニウスノルムに対する比が、予め定められた閾値以下となる更新段階をもって、コスト関数Ｊ_ｄの値が収束したとしてもよい。そして、診断システムは、Ｊ_ｄが収束した際のＧ、Ｑ、Ｗを、コスト関数Ｊ_ｄを最小化させる行列として求めてもよい。

＜実施形態２＞
本実施形態では、コスト関数Ｊ_ｄにおける係数δ_１、δ_２、ε_１、ε_２のうちの少なくとも１つが０以外の値である場合の処理について説明する。

（実験４）
発明者らは、実験３と同様の実験状況において、以下のような実験を行った。即ち、疵のない正常な状態にある軸受１０６と、図１３Ｃに示す外輪に微小な疵のある状態の軸受１０６と、について実験３で説明した実験と同様の実験を、式３９に示すコスト関数Ｊ_ｄのδ_１、δ_２、ε_２それぞれについて、値を変更しつつ行った。本実験では、使用固有値数ｓは、３とした。また、ｍは、５００とした。

本実験では、δ_２、ε_１、ε_２それぞれを０として、δ_１の値を、０、０．０００２５、０．０００５、０．０００７５、０．００１それぞれに変更しつつ、実験３で説明した実験と同様の実験を行った。

また、本実験では、δ_１、ε_１、ε_２それぞれを０として、δ_２の値を、０、０．００００５、０．０００２５、０．０００５それぞれに変更しつつ、実験３で説明した実験と同様の実験を行った。

また、本実験では、δ_１、δ_２、ε_１それぞれを０として、ε_２の値を、０、０．０２５、０．０５、０．１それぞれに変更しつつ、実験３で説明した実験と同様の実験を行った。

δ_２、ε_１、ε_２それぞれを０として、δ_１の値を変更しつつ、実験３で説明した実験と同様の実験を行った結果について説明する。δ_１は、コスト関数Ｊ_ｄにおける教師重み行列Ｇのスパース性の向上に関する制約の項に係る係数である。以下では、スパース性の向上に関する制約を、スパース制約とする。即ち、δ_１の値は、Ｙ_ｄの分解におけるＧのスパース制約の度合を示す係数である。

図２２Ａ〜図２２Ｅに、外輪に微小な疵のある状態の軸受１０６についての実験結果を示す。また、図２３Ａ〜図２３Ｅに、正常な状態の軸受１０６についての実験結果を示す。図２２Ａ〜図２２Ｅ、および図２３Ａ〜図２３Ｅの各グラフは、データ行列Ｙ_ｄの分解の結果求まったＧとＷとを示すグラフである。横軸は、ＧとＷとの各列のインデックスを示す。縦軸は、ＧとＷとの対応する列の要素の値を示す。

図２２Ａ、図２３Ａのグラフは、δ_１が０の場合の実験結果を示すグラフである。図２２Ｂ、図２３Ｂのグラフは、δ_１が０．０００２５の場合の実験結果を示すグラフである。図２２Ｃ、図２３Ｃのグラフは、δ_１が０．０００５の場合の実験結果を示すグラフである。図２２Ｄ、図２３Ｄのグラフは、δ_１が０．０００７５の場合の実験結果を示すグラフである。図２２Ｅ、図２３Ｅのグラフは、δ_１が０．００１の場合の実験結果を示すグラフである。

図２２Ａ〜図２２Ｅの各グラフを見ると、δ_１が０．０００２５、０．０００５、０．０００７５それぞれの場合、δ_１が０の場合に比べて、Ｇの値とＷの値とは、より分離しており、正常な状態と微小な疵のある状態との違いがより表れており、診断により有利になっていることが見て取れる。ただし、δ_１が０．００１の場合、δ_１が０．０００７５の場合に比べて、Ｇの値とＷの値とは、より近くなっており、正常な状態と微小な疵のある状態との違いが小さくなり、診断に不利になっていしまっていることが見て取れる。

図２３Ａ〜図２３Ｅの各グラフを見ると、δ_１が０．０００２５、０．０００５、０．０００７５、０．００１それぞれの場合で、Ｇの値とＷの値との分離の度合に顕著な違いが表れていないことが見て取れる。

発明者らは、図２２Ａ〜図２２Ｅ、図２３Ａ〜図２３Ｅの各グラフから、δ_１が０．０００２５、０．０００５、０．０００７５の場合において、軸受１０６の診断により有利であることを見出した。これにより、発明者らは、δ_１の適切な値は、０．０００２５以上０．０００７５以下の範囲の値であるとの知見を得た。

次に、δ_１、ε_１、ε_２それぞれを０として、δ_２の値を変更しつつ、求まった実験結果について説明する。δ_２は、コスト関数Ｊ^ｄにおける抽出重み行列Ｗのスパース制約の項に係る係数である。即ち、δ_２の値は、Ｙ_ｄの分解におけるＷのスパース制約の度合を示す係数である。

図２４Ａ〜図２４Ｄに、外輪に微小な疵のある状態の軸受１０６についての実験結果を示す。また、図２５Ａ〜図２５Ｄに、正常な状態の軸受１０６についての実験結果を示す。図２４Ａ〜図２４Ｄ、図２５Ａ〜図２５Ｄの各グラフは、データ行列Ｙ_ｄの分解の結果求まったＧとＷとを示すグラフである。横軸は、ＧとＷとの各列のインデックスを示す。縦軸は、ＧとＷとの対応する列の要素の値を示す。

図２４Ａ、図２５Ａのグラフは、δ_２が０の場合の実験結果を示すグラフである。図２４Ｂ、図２５Ｂのグラフは、δ_２が０．００００５の場合の実験結果を示すグラフである。図２４Ｃ、図２５Ｃのグラフは、δ_２が０．０００２５の場合の実験結果を示すグラフである。図２４Ｄ、図２５Ｄのグラフは、δ_２が０．０００５の場合の実験結果を示すグラフである。

図２４Ａ〜図２４Ｄの各グラフを見ると、δ_２が０．００００５の場合、δ_２が０の場合と同等に、Ｇの値とＷの値とが分離していることが見て取れる。また、δ_２が０．０００２５、０．０００５の場合、δ_２が０、０．００００５の場合に比べて、Ｇの値とＷの値とがより近くなっており、正常な状態と微小な疵のある状態との違いが小さくなり、診断に不利になっていしまっていることが見て取れる。

図２５Ａ〜図２５Ｄの各グラフを見ると、δ_２が０、０．００００５それぞれの場合で、Ｇの値とＷの値との分離の度合に顕著な違いが表れていないことが見て取れる。また、δ_２が０．０００２５、０．０００５それぞれの場合で、δ_２が０、０．００００５の場合に比べて、Ｇの値とＷの値とがより分離してしまっており、正常な状態同士で違いが大きくなり、診断に不利になっていしまっていることが見て取れる。

発明者らは、図２４Ａ〜図２４Ｂ、図２５Ａ〜図２５Ｂの各グラフから、δ_２が０．０００２５、０．０００５それぞれの場合において、軸受１０６の診断により不利であることを見出した。これにより、発明者らは、δ_２の適切な値は、０以上０．００００５以下（好ましくは０以上０．００００５未満）の範囲の値であるとの知見を得た。

次に、δ_１、δ_２、ε_１それぞれを０として、ε_２の値を変更しつつ、求まった実験結果について説明する。ε_２は、コスト関数Ｊ_ｄにおける教師基底と抽出基底との直交性の向上に関する制約の項に係る係数である。以下では、直交性の向上に関する制約を、直交制約とする。即ち、ε_２の値は、Ｙ_ｄの分解におけるＦとＱとの直交制約の度合を示す係数である。

図２６Ａ〜図２６Ｄに、外輪に微小な疵のある状態の軸受１０６についての実験結果を示す。また、図２７Ａ〜図２７Ｄに、正常な状態の軸受１０６についての実験結果を示す。図２６Ａ〜図２６Ｄ、図２７Ａ〜図２７Ｄの各グラフは、データ行列Ｙ_ｄの分解の結果求まったＧとＷとを示すグラフである。横軸は、ＧとＷとの各列のインデックスを示す。縦軸は、ＧとＷとの対応する列の要素の値を示す。

図２６Ａ、図２７Ａのグラフは、ε_２が０の場合の実験結果を示すグラフである。図２６Ｂ、図２７Ｂのグラフは、ε_２が０．０２５の場合の実験結果を示すグラフである。図２６Ｃ、図２６Ｃのグラフは、ε_２が０．０５の場合の実験結果を示すグラフである。図２６Ｄ、図２７Ｄのグラフは、ε_２が０．１の場合の実験結果を示すグラフである。

図２６Ａ〜図２６Ｄの各グラフを見ると、ε_２が０．０５、０．１の場合、ε_２が０、０．０２５の場合に比べて、Ｇの値とＷの値とがより近くなっており、正常な状態と微小な疵のある状態との違いが小さくなり、診断に不利になっていしまっていることが見て取れる。

図２７Ａ〜図２７Ｄの各グラフを見ると、ε_２が０．１の場合、ε_２が０、０．０２５、０．０５の場合に比べて、Ｇの値とＷの値とがより分離してしまっており、正常な状態同士で違いが大きくなり、診断に不利になっていしまっていることが見て取れる。

発明者らは、図２６Ａ〜図２６Ｄ、図２７Ａ〜図２７Ｄの各グラフから、ε_２が０．０５、０．１それぞれの場合において、軸受１０６の診断により不利であることを見出した。これにより、発明者らは、ε_２の適切な値は、０以上０．０２５以下の範囲の値であるとの知見を得た。

（診断処理）
本実施形態の診断システムのシステム構成は、実施形態１と同様である。また、情報処理装置１８００のハードウェ構成及び機能構成についても、実施形態１と同様である。

本実施形態の処理は、診断システムがδ_１、δ_２、ε_１、ε_２の全てが０のコスト関数Ｊ_ｄの代わりにδ_１、δ_２、ε_１、ε_２のうちの少なくとも１つが０でないコスト関数Ｊ_ｄを用いること以外は、図２０、図２１で説明した実施形態１の処理と同様である。

本実施形態では、δ_１の値は、実験により見出された適切な値の範囲（０．０００２５以上０．０００７５以下の範囲）から予め定められた値（例えば、０．０００２５、０．０００５、０．０００７５等）である。また、δ_２の値は、実験により見出された適切な値の範囲（０以上０．００００５以下の範囲）から予め定められた値（例えば、０、０．００００２５、０．００００５等）である。また、ε_２の値は、実験により見出された適切な値の範囲（０以上０．０２５以下の範囲）から予め定められた値（例えば、０、０．０１、０．０２５等）である。また、ε_１の値は、予め定められた値（例えば、０等）である。

δ_１、δ_２、ε_１、ε_２のうちの少なくとも１つが０でないコスト関数Ｊ_ｄは、δ_１、δ_２、ε_１、ε_２のうちの０でない係数に対応する制約を付加した関数となる。そのため、δ_１、δ_２、ε_１、ε_２のうちの少なくとも１つが０でないコスト関数Ｊ_ｄは、教師重み行列Ｇについてのスパース制約、抽出重み行列についてのスパース制約、抽出基底の正規直交制約、および教師基底と抽出基底との間の直交制約のうちの少なくとも１つを示す。

（まとめ）
本実施形態の処理により、診断システムは、診断対象の軸受１０６が予め定められた状態であるか否かをより精度よく診断できる。尚、理論上は、式３８に示す分解の自由度を低減させるために、コスト関数に抽出基底の正規直交制約と教師基底と抽出基底との間の直交制約とを付加するのが好ましい。本実施形態においても、実施形態１で説明した変形例を採用することができる。

＜その他の実施形態＞
本実施形態では、情報処理装置１８００は、補助記憶装置１９０２に記憶されたプログラムを実行することで、図２０、図２１の処理を実現することとした。しかし、情報処理装置１８００は、外付けの記憶媒体や外部の記憶サーバ等に記憶されたプログラムを実行することで、図２０、図２１の処理を実現することとしてもよい。

また、例えば、上述した情報処理装置１８００の機能の一部又は全てをハードウェアとして情報処理装置１８００に実装してもよい。
また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。即ち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。

本発明は、例えば、周期的な運動を行う物体の状態を診断することに利用することができる。

Claims

周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、
前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定手段と、
前記決定手段により決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、抽出基底行列と、抽出重み行列と、に分解する分解手段と、
前記分解手段により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断手段と、
を有し、
前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、
前記決定手段は、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、
前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σ_ｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕ_ｓと、から導出される行列Ｕ_ｓΣ_ｓＵ_ｓ ^Ｔであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、
前記教師基底行列は、教師基底を格納する非負値の行列であり、
前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、
前記抽出基底行列は、抽出基底を格納する非負値の行列であり、
前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、
前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、
前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、
前記分解手段は、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解する情報処理装置。
前記ｓ個の固有値は、前記自己相関行列の固有値のうち、値が最大の固有値を含む請求項１記載の情報処理装置。
前記ｓ個の固有値は、前記自己相関行列の固有値の中から、値が大きいものから順に選ばれた固有値である請求項２記載の情報処理装置。
前記自己相関行列は、前記計測データの予測値と当該計測データの予測値に対応する時刻における前記計測データの実測値との二乗誤差を最小化する条件を示す条件式に含まれる行列である請求項１乃至３何れか１項記載の情報処理装置。
前記コスト関数は、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分を評価する項と、前記教師重み行列のスパース性の向上に関する制約と、前記抽出重み行列のスパース性の向上に関する制約と、前記抽出基底の正規直交性の向上に関する制約と、前記教師基底と前記抽出基底との直交性の向上に関する制約と、のうちの少なくとも１つの制約を示す項とを含む請求項１乃至４何れか１項記載の情報処理装置。
前記コスト関数の値の適正化は、前記コスト関数の値の最小化また最大化である請求項１乃至５何れか１項記載の情報処理装置。
前記計測データは、前記物体の振動に係る計測データである請求項１乃至６何れか１項記載の情報処理装置。
前記診断手段による診断の結果に係る情報を出力する出力手段を更に有する請求項１乃至７何れか１項記載の情報処理装置。
前記物体は、鉄道台車に利用される軸受であり、
前記診断手段は、前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体の状態が、前記予め定められた状態である疵のない正常な状態であるか否かを診断する請求項１乃至８何れか１項記載の情報処理装置。
情報処理装置が実行する情報処理方法であって、
周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、
前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定ステップと、
前記決定ステップで決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、抽出基底行列と、抽出重み行列と、に分解する分解ステップと、
前記分解ステップでの分解により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断ステップと、
を含み、
前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、
前記決定ステップでは、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、
前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σ_ｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕ_ｓと、から導出される行列Ｕ_ｓΣ_ｓＵ_ｓ ^Ｔであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、
前記教師基底行列は、教師基底を格納する非負値の行列であり、
前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、
前記抽出基底行列は、抽出基底を格納する非負値の行列であり、
前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、
前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、
前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、
前記分解ステップでは、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解する情報処理方法。
周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、
前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定ステップと、
前記決定ステップで決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、前抽出基底行列と、抽出重み行列と、に分解する分解ステップと、
前記分解ステップでの分解により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断ステップと、
を含むステップをコンピュータに実行させ、
前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、
前記決定ステップでは、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、
前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σ_ｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕ_ｓと、から導出される行列Ｕ_ｓΣ_ｓＵ_ｓ ^Ｔであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、
前記教師基底行列は、教師基底を格納する非負値の行列であり、
前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、
前記抽出基底行列は、抽出基底を格納する非負値の行列であり、
前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、
前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、
前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、
前記分解ステップでは、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解するプログラム。