(57)【要約】【課題】行列の行と列の数に影響を受けずに、行列データの潜在的な構造を抽出することができる。【解決手段】データ解析部３２が、観測行列（Ｘl,m）L,Mの各要素Ｘl,mに対応するパラメータＷl,mに基づいて、各列ｌの［０，１］の範囲における仮想位置Ｕlと、各行ｍの［０，１］の範囲における仮想位置Ｕ'mとに対するランダム関数Θ（Ｕl，Ｕ'm）がパラメータＷl,mと等しくなるときに、［０，１］上に非負の値を持つ関数である基底Ｄと、［０，１］上に非負の値を持つ関数であるアクティベーションＨとの直積により、ランダム関数Θ（Ｕl，Ｕ'm）が生成されるように、基底Ｄに対応する潜在基底ｄと、アクティベーションＨに対応する潜在アクティベーションｈと、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、ランダム関数Θ（Ｕl，Ｕ'm）とを推定する。【選択図】図３
2018073234.tif 000002
【特許請求の範囲】
【請求項１】
観測行列（Ｘl,m）L,Mに潜在する構造を解析する行列データ解析装置であって、
前記観測行列（Ｘl,m）L,Mの各要素Ｘl,mに対応するパラメータＷl,mに基づいて、各列ｌの［０，１］の範囲における仮想位置Ｕlと、各行ｍの［０，１］の範囲における仮想位置Ｕ'mとに対するランダム関数Θ（Ｕl，Ｕ'm）がパラメータＷl,mと等しくなるときに、［０，１］上に非負の値を持つ関数である基底Ｄと、［０，１］上に非負の値を持つ関数であるアクティベーションＨとの直積により、前記ランダム関数Θ（Ｕl，Ｕ'm）が生成されるように、前記基底Ｄに対応する潜在基底ｄと、前記アクティベーションＨに対応する潜在アクティベーションｈと、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、前記ランダム関数Θ（Ｕl，Ｕ'm）とを推定するデータ解析部
を含む行列データ解析装置。
【請求項２】
前記データ解析部は、
各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、前記潜在基底ｄに対する補助変数εと、前記潜在アクティベーションｈに対する補助変数ηと、前記ランダム関数Θ（Ｕl，Ｕ'm）とを初期化するパラメータ初期化部と、
前記ランダム関数Θ（Ｕl，Ｕ'm）に基づいて、予め定められた補助変数ε及び補助変数ηの尤度関数を最適化するように、前記補助変数εと、前記補助変数ηとを更新する補助変数更新部と、
各列ｌの仮想位置Ｕlに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕlを更新し、
各行ｍの仮想位置Ｕ'mに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕ'mを更新する仮想位置更新部と、
各列ｌの仮想位置Ｕl、各行ｍの仮想位置Ｕ'm、及び前記パラメータＷl,mに基づいて、前記ランダム関数Θ（Ｕl，Ｕ'm）を更新する関数更新部と、
予め定められた終了条件を満たすまで、前記補助変数更新部、前記仮想位置更新部、及び前記関数更新部による更新を繰り返し行う終了判定部と、
を含む請求項１に記載の行列データ解析装置。
【請求項３】
前記補助変数更新部は、以下（１）式に従って、前記補助変数εを更新し、以下（２）式に従って、前記補助変数ηを更新する請求項２に記載の行列データ解析装置。
【数１】
2018073234.tif 000043
・・・（１）
2018073234.tif 000044
・・・（２）
ただし、ρは乱数を表し、
【数２】
2018073234.tif 000045
は、前記ランダム関数Θ（Ｕl，Ｕ'm）が与えられたもとでの前記補助変数ε及び補助変数ηの尤度関数である。
【請求項４】
観測行列（Ｘl,m）L,Mに潜在する構造を解析する行列データ解析装置における行列データ解析方法であって、
データ解析部が、前記観測行列（Ｘl,m）L,Mの各要素Ｘl,mに対応するパラメータＷl,mに基づいて、各列ｌの［０，１］の範囲における仮想位置Ｕlと、各行ｍの［０，１］の範囲における仮想位置Ｕ'mとに対するランダム関数Θ（Ｕl，Ｕ'm）がパラメータＷl,mと等しくなるときに、［０，１］上に非負の値を持つ関数である基底Ｄと、［０，１］上に非負の値を持つ関数であるアクティベーションＨとの直積により、前記ランダム関数Θ（Ｕl，Ｕ'm）が生成されるように、前記基底Ｄに対応する潜在基底ｄと、前記アクティベーションＨに対応する潜在アクティベーションｈと、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、前記ランダム関数Θ（Ｕl，Ｕ'm）とを推定するステップ
を含む行列データ解析方法。
【請求項５】
前記データ解析部が推定するステップは、
パラメータ初期化部が、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、前記潜在基底ｄに対する補助変数εと、前記潜在アクティベーションｈに対する補助変数ηと、前記ランダム関数Θ（Ｕl，Ｕ'm）とを初期化するステップと、
補助変数更新部が、前記ランダム関数Θ（Ｕl，Ｕ'm）に基づいて、予め定められた補助変数ε及び補助変数ηの尤度関数を最適化するように、前記補助変数εと、前記補助変数ηとを更新するステップと、
仮想位置更新部が、各列ｌの仮想位置Ｕlに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕlを更新するステップと、
仮想位置更新部が、各行ｍの仮想位置Ｕ' mに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕ'mを更新するステップと、
関数更新部とが、各列ｌの仮想位置Ｕl、各行ｍの仮想位置Ｕ'm、及び前記パラメータＷl,mに基づいて、前記ランダム関数Θ（Ｕl，Ｕ'm）を更新するステップ、
終了判定部が、予め定められた終了条件を満たすまで、前記補助変数更新部、前記仮想位置更新部、及び前記関数更新部による更新を繰り返し行うステップと、
を含む請求項４に記載の行列データ解析方法。
【請求項６】
前記補助変数更新部が更新するステップは、以下（３）式に従って、前記補助変数εを更新し、以下（４）式に従って、前記補助変数ηを更新する請求項５に記載の行列データ解析方法。
【数３】
2018073234.tif 000046
・・・（３）
2018073234.tif 000047
・・・（４）
ただし、ρは乱数を表し、
【数４】
2018073234.tif 000048
は、前記ランダム関数Θ（Ｕl，Ｕ'm）が与えられたもとでの前記補助変数ε及び補助変数ηの尤度関数である。
【請求項７】
コンピュータを、請求項１〜請求項３のいずれか１項に記載の行列データ解析装置の各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、行列データ解析装置、方法、及びプログラムに係り、特に、観測行列に潜在する構造を解析するための行列データ解析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
従来より、行列データ（ネットワークデータ、関係データを含む）の解析技術が知られている。入力された行列データの潜在的な構造を抽出することにより、欠損値の推定や予測などの各種行列データ解析のアプリケーションに利用することが出来る。
【０００３】
行列データの解析技術は、技術分野としては、広くは機械学習におけるノンパラメトリックベイズ分野に属し、関係データ解析・ネットワーク解析・行列分解・関数分解の要素技術と深く関連している。
【０００４】
行列データの解析技術に、ガウス過程による無限交換可能なグラフ表現がある（非特許文献１参照）。これは、行と列のインデックスの並び替えによって確率が変わらないような行列の確率モデルを表現するものである。行と列の各要素を［０，１］上の一様乱数に写像し、［０，１］２上のランダム関数を定義することによってモデルを構成する。
【０００５】
また、ガウス過程因子を用いた非負値行列因子分解がある（非特許文献２参照）。これは、入力行列を少数の主要構成因子に分解抽出する手法である。分解される因子にガウス過程事前モデルが導入されていることにより、所望の形に分解結果を誘導することが出来る。
【０００６】
また、交換可能性の表現定理がある。まず、二次元の配列
【０００７】
【数１】
2018073234.tif 000003
【０００８】
を考える。配列Ｚがseparately exchangeableであることは、
【０００９】
【数２】
2018073234.tif 000004
【００１０】
の任意の置換π，π'に対して
【００１１】
【数３】
2018073234.tif 000005
【００１２】
が成り立つことである。このような性質をもつ配列は以下のような形に表現出来ることが知られている。配列Ｚがseparately exchangeableであるとは以下と同値である。つまり、ある可測関数Ｇ：[０，１]３→Ｚがあり、Uniform[０，１]に従う一様変数Ｕirow,Ｕjcolumnand Ｕi,jに対して以下（１）式が成り立つ。
【００１３】
【数４】
2018073234.tif 000006
・・・（１）
【先行技術文献】
【非特許文献】
【００１４】
【非特許文献１】J. Lloyd, P. Orbanz, Z. Ghahramani, and D. M. Roy. Random function priors for exchangeable arrays with applications to graphs and relational data. In Advances in Neural Information Processing Systems, 2012.
【非特許文献２】M. N. Schmidt and H. Laurberg. Non-negative matrix factorization with Gaussian process priors. Computational Intelligence and Neuroscience, 2008.
【非特許文献３】D. J. Aldous. Representations for partially exchangeable arrays of random variables. Journal of Multivariate Analysis, 11:581-598, 1981.
【非特許文献４】D. N. Hoover. Relations on probability spaces and arrays of random variables.Technical report, Institute of Advanced Study, Princeton, 1979.
【発明の概要】
【発明が解決しようとする課題】
【００１５】
行列への解析結果はその入力行列の行と列の数に強く依存するため、データのサイズを変える度にそれに応じたアルゴリズムのパラメータチューニングを本質的に逃れることが出来ない。これは解析に用いる行列への確率的生成モデルが行と列のインデックスに対して交換可能性（インデックスの並び替えによって確率が変わらないこと）を満たしていないことに起因する。
【００１６】
本発明は、上記問題点を解決するために成されたものであり、行列の行と列の数に影響を受けずに、行列データの潜在的な構造を抽出することができる行列データ解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１７】
上記目的を達成するために、第１の発明に係る行列データ解析装置は、観測行列（Ｘl,m）L,Mに潜在する構造を解析する行列データ解析装置であって、前記観測行列（Ｘl,m）L,Mの各要素Ｘl,mに対応するパラメータＷl,mに基づいて、各列ｌの［０，１］の範囲における仮想位置Ｕlと、各行ｍの［０，１］の範囲における仮想位置Ｕ'mとに対するランダム関数Θ（Ｕl，Ｕ'm）がパラメータＷl,mと等しくなるときに、［０，１］上に非負の値を持つ関数である基底Ｄと、［０，１］上に非負の値を持つ関数であるアクティベーションＨとの直積により、前記ランダム関数Θ（Ｕl，Ｕ'm）が生成されるように、前記基底Ｄに対応する潜在基底ｄと、前記アクティベーションＨに対応する潜在アクティベーションｈと、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、前記ランダム関数Θ（Ｕl，Ｕ'm）とを推定するデータ解析部、を含んで構成されている。
【００１８】
また、第１の発明に係る行列データ解析装置において、前記データ解析部は、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、前記潜在基底ｄに対する補助変数εと、前記潜在アクティベーションｈに対する補助変数ηと、前記ランダム関数Θ（Ｕl，Ｕ'm）とを初期化するパラメータ初期化部と、前記ランダム関数Θ（Ｕl，Ｕ'm）に基づいて、予め定められた補助変数ε及び補助変数ηの尤度関数を最適化するように、前記補助変数εと、前記補助変数ηとを更新する補助変数更新部と、各列ｌの仮想位置Ｕlに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕlを更新し、各行ｍの仮想位置Ｕ'mに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕ'mを更新する仮想位置更新部と、各列ｌの仮想位置Ｕl、各行ｍの仮想位置Ｕ'm、及び前記パラメータＷl,mに基づいて、前記ランダム関数Θ（Ｕl，Ｕ'm）を更新する関数更新部と、予め定められた終了条件を満たすまで、前記補助変数更新部、前記仮想位置更新部、及び前記関数更新部による更新を繰り返し行う終了判定部と、を含むようにしてもよい。
【００１９】
また、第１の発明に係る行列データ解析装置において、前記補助変数更新部は、以下（２９）式に従って、前記補助変数εを更新し、以下（３０）式に従って、前記補助変数ηを更新するようにしてもよい。
【数５】
2018073234.tif 000007
・・・（２９）
2018073234.tif 000008
・・・（３０）
ただし、ρは乱数を表し、
【数６】
2018073234.tif 000009
は、前記ランダム関数Θ（Ｕl，Ｕ'm）が与えられたもとでの前記補助変数ε及び補助変数ηの尤度関数である。
【００２０】
第２の発明に係る行列データ解析方法は、観測行列（Ｘl,m）L,Mに潜在する構造を解析する行列データ解析装置における行列データ解析方法であって、データ解析部が、前記観測行列（Ｘl,m）L,Mの各要素Ｘl,mに対応するパラメータＷl,mに基づいて、各列ｌの［０，１］の範囲における仮想位置Ｕlと、各行ｍの［０，１］の範囲における仮想位置Ｕ'mとに対するランダム関数Θ（Ｕl，Ｕ'm）がパラメータＷl,mと等しくなるときに、［０，１］上に非負の値を持つ関数である基底Ｄと、［０，１］上に非負の値を持つ関数であるアクティベーションＨとの直積により、前記ランダム関数Θ（Ｕl，Ｕ'm）が生成されるように、前記基底Ｄに対応する潜在基底ｄと、前記アクティベーションＨに対応する潜在アクティベーションｈと、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、前記ランダム関数Θ（Ｕl，Ｕ'm）とを推定するステップを含んで実行することを特徴とする。
【００２１】
また、第２の発明に係る行列データ解析方法において、前記データ解析部が推定するステップは、パラメータ初期化部が、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、前記潜在基底ｄに対する補助変数εと、前記潜在アクティベーションｈに対する補助変数ηと、前記ランダム関数Θ（Ｕl，Ｕ'm）とを初期化するステップと、補助変数更新部が、前記ランダム関数Θ（Ｕl，Ｕ'm）に基づいて、予め定められた補助変数ε及び補助変数ηの尤度関数を最適化するように、前記補助変数εと、前記補助変数ηとを更新するステップと、仮想位置更新部が、各列ｌの仮想位置Ｕlに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕlを更新するステップと、仮想位置更新部が、各行ｍの仮想位置Ｕ'mに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕ'mを更新するステップと、関数更新部とが、各列ｌの仮想位置Ｕl、各行ｍの仮想位置Ｕ'm、及び前記パラメータＷl,mに基づいて、前記ランダム関数Θ（Ｕl，Ｕ'm）を更新するステップ、終了判定部が、予め定められた終了条件を満たすまで、前記補助変数更新部、前記仮想位置更新部、及び前記関数更新部による更新を繰り返し行うステップと、を含むようにしてもよい。
【００２２】
また、第２の発明に係る行列データ解析方法において、前記補助変数更新部が更新するステップは、以下（３１）式に従って、前記補助変数εを更新し、以下（３２）式に従って、前記補助変数ηを更新するようにしてもよい。
【数７】
2018073234.tif 000010
・・・（３１）
2018073234.tif 000011
・・・（３２）
ただし、ρは乱数を表し、
【数８】
2018073234.tif 000012
は、前記ランダム関数Θ（Ｕl，Ｕ'm）が与えられたもとでの前記補助変数ε及び補助変数ηの尤度関数である。
【００２３】
また、第３の発明に係るプログラムは、コンピュータを、第１の発明に係る行列データ解析装置の各部として機能させるためのプログラムである。
【発明の効果】
【００２４】
本発明の行列データ解析装置、方法、及びプログラムによれば、観測行列（Ｘl,m）L,Mの各要素Ｘl,mに対応するパラメータＷl,mに基づいて、各列ｌの［０，１］の範囲における仮想位置Ｕlと、各行ｍの［０，１］の範囲における仮想位置Ｕ'mとに対するランダム関数Θ（Ｕl，Ｕ'm）がパラメータＷl,mと等しくなるときに、［０，１］上に非負の値を持つ関数である基底Ｄと、［０，１］上に非負の値を持つ関数であるアクティベーションＨとの直積により、ランダム関数Θ（Ｕl，Ｕ'm）が生成されるように、基底Ｄに対応する潜在基底ｄと、アクティベーションＨに対応する潜在アクティベーションｈと、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、ランダム関数Θ（Ｕl，Ｕ'm）とを推定することにより、行列の行と列の数に影響を受けずに、行列データの潜在的な構造を抽出することができる、という効果が得られる。
【図面の簡単な説明】
【００２５】
【図１】ＵＬ∈[０,１],ＵＭ∈[０,１]の直観的なイメージを表した図である。
【図２】Θの生成モデルの直観的なイメージを表した図である。
【図３】本発明の実施の形態に係る行列データ解析装置の構成を示すブロック図である。
【図４】Animal feature データセットを用いて解析したΘの一例を示す図である。
【図５】Animal feature データセットを用いて解析したΘの２つの因子（Ｄ，Ｈ）の一例を示す図である。
【図６】本発明の実施の形態に係る行列データ解析装置における行列データ解析処理ルーチンを示すフローチャートである。
【発明を実施するための形態】
【００２６】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００２７】
＜本発明の実施の形態に係る原理＞
【００２８】
以下、本発明の実施の形態における原理を説明する。
【００２９】
まず、Aldous-Hoover定理により、交換可能行列（Ｘｌ,ｍ）Ｌ,Ｍは[０，１]２上のランダム関数Ｇのモデル化について説明する。
【００３０】
まず観測行列（Ｘｌ,ｍ）Ｌ,Ｍの尤度関数の設計から始める。Aldous-Hoover定理に従い、補助変数Ｕ１,Ｕ２,...,ＵＬ, Ｕ１,Ｕ２,...,ＵＭを導入して、以下（２）式及び（３）式のように独立な[０，１]上の一様乱数とする。
【００３１】
【数９】
2018073234.tif 000013
・・・（２）
2018073234.tif 000014
・・・（３）
【００３２】
この補助変数の導入により、観測行列は[０，１]２の長方形を介して表現することが出来るようになる。ここで[０，１]２上のランダム関数Ｇを次のような階層モデルによって表現することで観測行列の生成モデルとする。
【００３３】
【数１０】
2018073234.tif 000015
・・・（４）
2018073234.tif 000016
・・・（５）
【００３４】
ここで、Ｐ（・｜・）は任意の尤度関数を用いることができる。Θ:[０，１]２→[０，１]は尤度関数に用いるパラメータへの事前分布とする。例えば、観測行列の各要素がバイナリ（０または１）である場合は、尤度関数にベルヌーイ分布、事前分布にはベータ分布を採用するのが標準的である。
【００３５】
次にΘへの確率的生成モデルを考える。ここで、Θは１ランク非負値因子の重ね合わせで表現出来ると仮定する。それぞれの因子は[０，１]上に非負の値を持つ基底と、同じく[０，１]上に非負の値を持つアクティベーションの積で表すことが出来る。このような基底とアクティベーションを以下のように設計する。まず[０，Ｋ]上に実数の値を持つｈを導入し、これが０を平均としｈを共分散とするガウス過程から生成されたものとみなす。
【００３６】
【数１１】
2018073234.tif 000017
・・・（６）
【００３７】
より正確にはｈは次の確率分布に従う。
【００３８】
【数１２】
2018073234.tif 000018
・・・（７）
【００３９】
次にこのｈを非負値因子Ｈ(Ｋ×[０，１])に変換することを考える。非負値
【００４０】
【数１３】
2018073234.tif 000019
【００４１】
は以下のリンク関数を介してｈ(ｋ＋ｕ)から得ることが出来る。
【００４２】
【数１４】
2018073234.tif 000020
・・・（８）
【００４３】
ただし、vec()は[０，１]上のＫ個の実数値の集合を[０，Ｋ]上の実数値の集合へ変換する写像とする。具体的にはｈ(ｋ＋ｕ)＝ ｆｈ(Ｈｋ(ｕ))を意味する。リンク関数としては例えば次のものを用いることが出来る。
【００４４】
【数１５】
2018073234.tif 000021
・・・（９）
【００４５】
ここでerfは誤差関数、σiはｈのｉ番目の対角要素、λは重みの定数を表している。以上のアクティベーションと同様に、基底Ｄも次のように構成することが出来る。
【００４６】
【数１６】
2018073234.tif 000022
・・・（１０）
2018073234.tif 000023
・・・（１１）
【００４７】
これらを用いてΘは次のように表すことが出来る。
【００４８】
【数１７】
2018073234.tif 000024
・・・（１２）
【００４９】
ただし、
【００５０】
【数１８】
2018073234.tif 000025
【００５１】
は(Ｕｌ,Ｕ'０)∈[０,１]２上のＤＨの値を表している。
【００５２】
図１にＵＬ∈[０,１],ＵＭ∈[０,１]の直観的なイメージを表した図を示す。観測行列の各要素Ｘl,mはΘ（Ｕｌ,Ｕ'０）を与えられたもとでの尤度関数から生成される。
【００５３】
また、図２にΘの生成モデルの直観的なイメージを表した図を示す。基底Ｄ（[０,１]×Ｋ）とアクティベーションＨ（Ｋ×[０,１]）の直積により構成される。
【００５４】
[アルゴリズムの説明]
【００５５】
ここでは提案モデルを観測データにフィッティングさせる方法としてマルコフ連鎖モンテカルロ法による推論手法の一例を示す。この手法は反復的なアルゴリズムに基づいており、潜在アクティベーションｈ、潜在基底ｄ、仮想位置Ｕｌ,Ｕ'０、ランダム関数Θを逐次的に更新していることに相当する。
【００５６】
リンク関数の逆関数（ｆｈ−１）'の一階微分は以下のように表すことが出来る。
【００５７】
【数１９】
2018073234.tif 000026
・・・（１３）
【００５８】
ここで補助変数εとηを導入し以下のような変形を考える。
【００５９】
【数２０】
2018073234.tif 000027
・・・（１４）
2018073234.tif 000028
・・・（１５）
【００６０】
ただし、ｆｄ−１及びｆｈ−１はリンク関数ｆｄ、ｆｈの逆関数として与えられ、ＣｄとＣｈはそれぞれ共分散行列Σｄ,Σｈの行列平方根であり、Σｄ,Σｈのチョムスキー分解を表している。このとき補助変数の確率分布は以下に従う。
【００６１】
【数２１】
2018073234.tif 000029
・・・（１６）
2018073234.tif 000030
・・・（１７）
【００６２】
これにより、Θが与えられたもとで、εとηの尤度関数は以下のように表すことが出来る。
【００６３】
【数２２】
2018073234.tif 000031
・・・（１８）
【００６４】
この一次微分はηについては、
【００６５】
【数２３】
2018073234.tif 000032
・・・（１９）
2018073234.tif 000033
・・・（２０）
2018073234.tif 000034
・・・（２１）
【００６６】
となる。同様にεについては、
【００６７】
【数２４】
2018073234.tif 000035
・・・（２２）
2018073234.tif 000036
・・・（２３）
2018073234.tif 000037
・・・（２４）
【００６８】
となる。
【００６９】
以上のような変数変換により実質的なアルゴリズムは、補助変数ε、補助変数η、仮想位置Ｕｌ(ｌ＝１,２,...,Ｌ),Ｕ'０(ｍ＝１,２,...,Ｍ)となる。
【００７０】
以下の表１に各パラメータを示す。
【００７１】
【表１】
2018073234.tif 000038
【００７２】
＜本発明の実施の形態に係る行列データ解析装置の構成＞
【００７３】
次に、本発明の実施の形態に係る行列データ解析装置の構成について説明する。図３に示すように、本発明の実施の形態に係る行列データ解析装置１００は、ＣＰＵと、ＲＡＭと、後述する行列データ解析処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この行列データ解析装置１００は、機能的には図３に示すように入力部１０と、演算部２０と、出力部５０とを備えている。
【００７４】
入力部１０は、観測行列（Ｘl,m）L,Mを受け付ける。
【００７５】
演算部２０は、データ解析部３２を含んで構成されている。
【００７６】
データ解析部３２は、以下に説明する各部の処理により、観測行列（Ｘl,m）L,Mの各要素Ｘl,mに対応するパラメータＷl,mに基づいて、各列ｌの［０，１］の範囲における仮想位置Ｕlと、各行ｍの［０，１］の範囲における仮想位置Ｕ'mとに対するランダム関数Θ（Ｕl，Ｕ'm）がパラメータＷl,mと等しくなるときに、［０，１］上に非負の値を持つ関数である基底Ｄと、［０，１］上に非負の値を持つ関数であるアクティベーションＨとの直積により、ランダム関数Θ（Ｕl，Ｕ'm）が生成されるように、基底Ｄに対応する潜在基底ｄと、アクティベーションＨに対応する潜在アクティベーションｈと、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、ランダム関数Θ（Ｕl，Ｕ'm）とを推定する。
【００７７】
データ解析部３２は、パラメータ初期化部３４と、補助変数更新部３６と、仮想位置更新部３８と、関数更新部４０と、終了判定部４２とを含んで構成されている。
【００７８】
パラメータ初期化部３４は、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、潜在基底ｄに対する補助変数εと、潜在アクティベーションｈに対する補助変数ηと、ランダム関数Θ（Ｕl，Ｕ'm）とを初期化する。また、入力部１０で受け付けた観測行列（Ｘl,m）L,Mの各要素Ｘl,mに対応するパラメータＷl,mを設定する。具体的には、観測行列（Ｘl,m）L,Mの各要素Ｘl,mの値を、そのままパラメータＷl,mに設定する。
【００７９】
補助変数ε及び補助変数ηは、下記の確率分布からサンプリングして初期化する。
【００８０】
【数２５】
2018073234.tif 000039
・・・（２５）
2018073234.tif 000040
・・・（２６）
【００８１】
仮想位置Ｕl(ｌ＝１,２,...,Ｌ),Ｕ0(ｍ＝１,２,...,Ｍ)は、［０,１］上の一様分布からサンプリングして初期化する。
【００８２】
ランダム関数Θ（Ｕl，Ｕ'm）は、仮想位置Ｕlと、仮想位置Ｕ'mと、パラメータＷl,mとに基づいて初期化する。
【００８３】
補助変数更新部３６は、ランダム関数Θ（Ｕl，Ｕ'm）に基づいて、予め定められた補助変数ε及び補助変数ηの尤度関数を最適化するように、補助変数εと、補助変数ηとを更新する。
【００８４】
具体的には、補助変数εについて、現時点のサンプルとランダム関数Θ（Ｕl，Ｕ'm）を使って、次式に含まれる尤度関数を最適化するように新しいサンプルを生成する。
【００８５】
【数２６】
2018073234.tif 000041
・・・（２７）
【００８６】
補助変数ηについては、現時点のサンプルとランダム関数Θ（Ｕl，Ｕ'm）を使って、次式に含まれる尤度関数を最適化するように新しいサンプルを生成する。
【００８７】
【数２７】
2018073234.tif 000042
・・・（２８）
【００８８】
ただし、ρは乱数を表す。
【００８９】
ここで、補助変数εと、補助変数ηとが求まることによって、上記（１４）式及び（１５）式に示すように、基底Ｄ及びアクティベーションＨが直ちに求められる。なお、最初に補助変数ε及び補助変数ηを更新したときには、Ｃｄ及びＣｈは考慮せずに、上記（１４）式及び（１５）式を計算して基底Ｄ及びアクティベーションＨを求める。また、基底Ｄ及びアクティベーションＨを求める度に、基底Ｄ及びアクティベーションＨに基づいてΣｄ及びΣｈが得られ、Ｃｄ及びＣｈが求まるため、次の更新において用いるようにすればよい。
【００９０】
仮想位置更新部３８は、各列ｌの仮想位置Ｕlに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕlを更新する。また、各行ｍの仮想位置Ｕ'mに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕ'mを更新する。
【００９１】
仮想位置Ｕlについては、具体的には、まずｌ番目の行ごとに、［０，１］上のインテンシティ１のポアソン過程からサンプル点を生成する。これは［０，１］の数直線を分割していることに相当する。そして、ｌ番目の行が［０，１］の分割の各クラスに属した場合の事後確率を計算し、Ｕｌが属するクラスをこの事後確率から生成する。決定したクラス上の数直線の一様分布からＵｌを生成する。
【００９２】
仮想位置Ｕ'mについては、具体的には、まずｍ番目の列ごとに、［０，１］上のインテンシティ１のポアソン過程からサンプル点を生成する。これは［０，１］の数直線を分割していることに相当する。そして、ｍ番目の列が［０，１］の分割の各クラスに属した場合の事後確率を計算し、Ｕ'mが属するクラスをこの事後確率から生成する。決定したクラス上の数直線の一様分布からＵ'mを生成する。
【００９３】
関数更新部４０は、仮想位置Ｕlと、仮想位置Ｕ'mと、パラメータＷl,mとに基づいてランダム関数Θ（Ｕl，Ｕ'm）を更新する。上記（５）式の尤度関数は観測行列Ｘに応じて選ぶことができ、例えば観測行列の各要素がバイナリ（０または１）である場合は、尤度関数にベルヌーイ分布、事前分布にはベータ分布を採用するのが標準的である。尤度関数の選ぶ方によらず非特許文献１の5.2節または5.3節の方法により、新しいランダム関数Θ（Ｕl，Ｕ'm）のサンプルを一つ生成する。直感的には離散行列をなめらかな長方形に変換するイメージで、新しいランダム関数Θ（Ｕl，Ｕ'm）のサンプルを一つ生成する。
【００９４】
終了判定部４２は、予め定められた終了条件を満たすまで、補助変数更新部３６、仮想位置更新部３８、及び関数更新部４０による更新を繰り返し行う。本実施の形態では、終了条件として、事前に指定した既定の回数の反復が終了したかを判定する。終了時の仮想位置Ｕlと、仮想位置Ｕ'mと、補助変数εと、補助変数ηと、ランダム関数Θ（Ｕl，Ｕ'm）とを出力部５０に出力する。
【００９５】
図４に、Animal feature データセットを用いて解析したΘの一例を示す。
【００９６】
また図５に、Animal feature データセットを用いて解析したΘの２つの因子（Ｄ，Ｈ）の一例を示す。
【００９７】
＜本発明の実施の形態に係る行列データ解析装置の作用＞
【００９８】
次に、本発明の実施の形態に係る行列データ解析装置１００の作用について説明する。入力部１０において観測行列（Ｘl,m）L,Mを受け付けると、行列データ解析装置１００は、図６に示す行列データ解析処理ルーチンを実行する。
【００９９】
まず、ステップＳ１００では、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、潜在基底ｄに対する補助変数εと、潜在アクティベーションｈに対する補助変数ηと、ランダム関数Θ（Ｕl，Ｕ'm）とを初期化する。また、入力部１０で受け付けた観測行列（Ｘl,m）L,Mの各要素Ｘl,mに対応するパラメータＷl,mを設定する。
【０１００】
ステップＳ１０２では、ランダム関数Θ（Ｕl，Ｕ'm）に基づいて、上記（２７）式に従って、補助変数εの尤度関数を最適化するように、補助変数εを更新する。
【０１０１】
ステップＳ１０４では、ランダム関数Θ（Ｕl，Ｕ'm）に基づいて、上記（２８）式に従って、補助変数ηの尤度関数を最適化するように、補助変数ηとを更新する。
【０１０２】
ステップＳ１０６では、各列ｌの仮想位置Ｕlに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕlを更新する。
【０１０３】
ステップＳ１０８では、各行ｍの仮想位置Ｕ'mに対し、ポアソン過程から生成されたサンプル点で、［０，１］の範囲を分割したときの各クラスに属する事後確率の分布から、仮想位置Ｕlが属するクラスを決定し、仮想位置Ｕ'mを更新する。
【０１０４】
ステップＳ１１０では、仮想位置Ｕlと、仮想位置Ｕ'mと、パラメータＷl,mとに基づいてランダム関数Θ（Ｕl，Ｕ'm）を更新する。
【０１０５】
ステップＳ１１２では、事前に指定した既定の回数の反復が終了したかを判定し、終了していれば、終了時の仮想位置Ｕlと、仮想位置Ｕ'mと、補助変数εと、補助変数ηと、ランダム関数Θ（Ｕl，Ｕ'm）とを出力部５０に出力して処理を終了し、終了していなければステップＳ１０２に戻って処理を繰り返す。
【０１０６】
以上説明したように、本発明の実施の形態に係る行列データ解析装置によれば、観測行列（Ｘl,m）L,Mの各要素Ｘl,mに対応するパラメータＷl,mに基づいて、各列ｌの［０，１］の範囲における仮想位置Ｕlと、各行ｍの［０，１］の範囲における仮想位置Ｕ'mとに対するランダム関数Θ（Ｕl，Ｕ'm）がパラメータＷl,mと等しくなるときに、［０，１］上に非負の値を持つ関数である基底Ｄと、［０，１］上に非負の値を持つ関数であるアクティベーションＨとの直積により、ランダム関数Θ（Ｕl，Ｕ'm）が生成されるように、基底Ｄに対応する潜在基底ｄと、アクティベーションＨに対応する潜在アクティベーションｈと、各列ｌの仮想位置Ｕlと、各行ｍの仮想位置Ｕ'mと、ランダム関数Θ（Ｕl，Ｕ'm）とを推定することにより、行列の行と列の数に影響を受けずに、行列データの潜在的な構造を抽出することができる。
【０１０７】
また、行列データの中から主要要素を見つける因子分解において、従来非特許文献２型の観測データの行列分解モデルを用いられることがあったが、これは限られた有限の観測データに過学習してしまう問題を本質的に逃れることが出来ない。本発明の実施の形態の手法によれば、Aldous-Hooverの交換可能行列の表現定理（非特許文献３）を利用して有限の観測データの背後にある潜在的な[０,１]２の上の写像に対し、Ｄ及びＨの関数分解を行うことによってこの問題を解決している。
【０１０８】
一方、交換可能行列に対して潜在的な[０,１]２上の写像へのモデルを介して表現する方法が非特許文献１で提案されているが、その写像の自由度は極めて高く、何らかの性質を埋め込むことが容易ではなかった。本発明の実施の形態の手法では、その写像に対し少数の因子で構成されているという仮定のもと、ガウス過程の直積型のモデルを導入したことで過学習を防ぐことを提案した。
【０１０９】
なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１１０】
例えば、上述した実施の形態では、観測行列（Ｘl,m）L,Mの各要素Ｘl,mの値を、そのままパラメータＷl,mに設定していたが、これに限定されるものではなく、観測行列（Ｘl,m）L,Mの各要素Ｘl,mに基づいて、上記（５）式の尤度関数からパラメータＷl,mを求めるようにしてもよい。
【符号の説明】
【０１１１】
１０   入力部
２０   演算部
３２   データ解析部
３４   パラメータ初期化部
３６   補助変数更新部
３８   仮想位置更新部
４０   関数更新部
４２   終了判定部
５０   出力部
１００ 行列データ解析装置
【図１】
2018073234.tif 000049
【図２】
2018073234.tif 000050
【図３】
2018073234.tif 000051
【図４】
2018073234.tif 000052
【図５】
2018073234.tif 000053
【図６】
2018073234.tif 000054