(57)【要約】【課題】観測行列データを冗長さのない任意の長方形に分割することができるようにする。【解決手段】観測行列データ＾ＸＩを分割する際、長方形分割更新部２２による観測行列データ＾ＸＩの長方形分割Ｒの更新と、中間変数更新部２３による観測行列データ＾ＸＩと同じサイズの行列＾Ｉの長方形分割н及びバイナリ変数＾ｌの更新と、クラスタ更新部２４によるクラスタ＾ξ、及びクラスタ＾ηの更新とを繰り返し実施する。【選択図】図２
2015088117.tif 000002
【特許請求の範囲】
【請求項１】
入力された観測行列データ＾ＸＩと同じサイズの行列＾Ｉを少なくとも１つのブロックに分割したときの長方形分割нと、前記観測行列データ＾ＸＩを少なくとも１つのブロックに分割したときの長方形分割Ｒと、前記観測行列データ＾ＸＩにおける各行ｍが所属する各クラスタξｍと、前記観測行列データ＾ＸＩにおける各列ｎが所属する各クラスタηｎと、の各々に対して初期値を設定するパラメータ初期化部と、
前記長方形分割нに基づいて、前記長方形分割Ｒの新しい候補Ｒｎｅｗを求め、前記観測行列データ＾ＸＩ、前記長方形分割н、前記クラスタξｍ、及び前記クラスタηｎに基づいて算出される前記長方形分割Ｒの事後確率と前記新しい候補Ｒｎｅｗの事後確率とを比較して、前記新しい候補Ｒｎｅｗを採択するか否かを決定して、前記長方形分割Ｒを更新する長方形分割変更部と、
前記観測行列データ＾ＸＩ、前記長方形分割Ｒ、及び前記長方形分割нに基づいて、前記長方形分割нを更新する中間変数更新部と、
前記観測行列データ＾ＸＩ、前記長方形分割Ｒ、及び前記長方形分割нに基づいて、前記観測行列データ＾ＸＩにおける各行ｍのクラスタξｍと各列ｎのクラスタηｎとを更新するクラスタ更新部と、
予め定められた終了条件を満たすまで、前記長方形分割更新部による更新、前記中間変数更新部による更新、及び前記クラスタ更新部による更新を繰り返し行う終了判定部と、
を含む行列データ解析装置。
【請求項２】
前記パラメータ初期化部は、前記長方形分割нの各ブロックＨに対するバイナリ変数ｌＨの初期値を更に設定し、
前記長方形分割更新部は、前記長方形分割нの各ブロックＨについて、前記ブロックＨのバイナリ変数ｌＨに基づいて、行方向又は列方向の順序を優先して前記ブロックＨ内の要素の全順序を更新することにより、前記行列＾Ｉの要素の全順序を更新し、前記更新された前記行列＾Ｉの要素の全順序に基づいて、前記行列＾Ｉの各要素について、何れかの隣接要素と同じブロックに割り当てるか、または新しいブロックに割り当てるかを決定することにより、前記行列＾Ｉをブロックに分割し、前記行列＾Ｉの長方形分割を、前記長方形分割Ｒの新しい候補Ｒｎｅｗとして求め、前記長方形分割Ｒの事後確率と前記新しい候補Ｒｎｅｗの事後確率とを比較して、前記新しい候補Ｒｎｅｗを採択するか否かを決定して、前記長方形分割Ｒを更新し、
前記中間変数更新部は、前記観測行列データ＾ＸＩ、前記長方形分割Ｒ、及び前記長方形分割нに基づいて、前記長方形分割н及び前記バイナリ変数ｌＨを更新する
請求項１記載の行列データ解析装置。
【請求項３】
前記中間変数更新部は、長方形分割нに基づいて、前記長方形分割нのｋ番目のブロックを更に分割するか、または前記長方形分割нの隣接する２つのブロックを１つのブロックに統合するかを決定して、新たな長方形分割нを生成し、前記長方形分割нの事後確率と前記新たな長方形分割нの事後確率とに基づいて、前記新たな長方形分割нを採択するか否かを決定して、前記長方形分割нを更新する
請求項１又は請求項２記載の行列データ解析装置。
【請求項４】
前記中間変数更新部は、前記観測行列データ＾ＸＩ、前記長方形分割Ｒ、及び前記長方形分割нに基づいて、前記長方形分割нの各ブロックＨについて、前記ブロックＨに対するバイナリ変数ｌＨが０又は１を取る条件付き確率を計算し、計算された条件付き確率に基づいて、前記ブロックＨに対するバイナリ変数ｌＨを更新する
請求項２記載の行列データ解析装置。
【請求項５】
パラメータ初期化部によって、入力された観測行列データ＾ＸＩと同じサイズの行列＾Ｉを少なくとも１つのブロックに分割したときの長方形分割нと、前記観測行列データ＾ＸＩを少なくとも１つのブロックに分割したときの長方形分割Ｒと、前記観測行列データ＾ＸＩにおける各行ｍが所属する各クラスタξｍと、前記観測行列データ＾ＸＩにおける各列ｎが所属する各クラスタηｎと、の各々に対して初期値を設定し、
長方形分割変更部によって、前記長方形分割нに基づいて、前記長方形分割Ｒの新しい候補Ｒｎｅｗを求め、前記観測行列データ＾ＸＩ、前記長方形分割н、前記クラスタξｍ、及び前記クラスタηｎに基づいて算出される前記長方形分割Ｒの事後確率と前記新しい候補Ｒｎｅｗの事後確率とを比較して、前記新しい候補Ｒｎｅｗを採択するか否かを決定して、前記長方形分割Ｒを更新し、
中間変数更新部によって、前記観測行列データ＾ＸＩ、前記長方形分割Ｒ、及び前記長方形分割нに基づいて、前記長方形分割нを更新し、
クラスタ更新部によって、前記観測行列データ＾ＸＩ、前記長方形分割Ｒ、及び前記長方形分割нに基づいて、前記観測行列データ＾ＸＩにおける各行ｍのクラスタξｍと各列ｎのクラスタηｎとを更新し、
終了判定部によって、予め定められた終了条件を満たすまで、前記長方形分割更新部による更新、前記中間変数更新部による更新、及び前記クラスタ更新部による更新を繰り返す
行列データ解析方法。
【請求項６】
コンピュータを、請求項１〜請求項４の何れか１項に記載の行列データ解析装置の各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、行列データ解析装置、方法、及びプログラムに係り、特に行列データ、例えば、商品×購買者で商品購入個数を表す離散値行列や、人×人の友人関係ネットワークを表すバイナリ値行列等から長方形分割によるクラスタを抽出するための行列データ解析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
行列データをクラスタリングして抽出する手法の代表的なものとして、長方形分割による解析方法が知られている。
【０００３】
従来の長方形分割による行列データ解析方法は、図１（Ａ）に示すように、行列データの行と列がそれぞれクラスタリングされ、その直積によって各クラスタが特徴づけられる“regular grid” と呼ばれるクラスに限定されるものであった（非特許文献１）。
【０００４】
また、他の長方形分割による行列データ解析方法として、ポアソン過程の多次元への拡張であるモンドリアン過程と呼ばれるモデルが知られている（非特許文献２）。ここでは、モンドリアン過程に関する説明を簡単にするため２次元に限定して述べる。
【０００５】
２次元の領域上のモンドリアン過程Ｍは以下のような再帰的な手順でこの領域に対して確率的な分割を与える。はじめにbudgetと呼ばれる正の実数を予め指定する。各再帰の中でcostが生成され、budgetとの比較を行う。もしcostがbudgetを超過した場合、再帰的な反復は終了し現在の分割を返す。もしcostがbudgetを超過しなかった場合は新しい分割を追加し領域が２つの領域に分割され、budgetからcostが引かれる。そして、分割された２つの領域に各々独立なモンドリアン過程が再帰的に適用されていく。
【０００６】
この長方形分割の様子を示したものが図１（Ｂ）である。この場合の長方形分割は、行列データを水平方向もしくは垂直方向に分割していく行為を再帰的に行う“hierarchical”と呼ばれるクラスに限定される。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】Kemp, C. and Tenenbaum, J. B. and Griffiths, T. L. and Yamada,T. and Ueda, N. (2006). "Learning systems of concepts with an infinite relational model," Proceedings of the National Conference onArtificial Intelligence.
【非特許文献２】D. Roy, Y. W. Teh (2009). "The Mondrian Process," Advances in Neural Information Processing Systems.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
このように、従来では行列データを長方形にクラスタリングする技術においては、“regular grid”と“hierarchical”と呼ばれる限られたクラスに属する長方形分割しか抽出できず、全ての長方形分割を抽出しうる能力を持った手法は存在しなかった。
【０００９】
これは、本来の行列データが持つクラスタを表現するためには冗長な分割をする必要があったことを意味する。すなわち、ある長方形分割で“regular grid”と“hierarchical”に属さないもののクラスタを表現するため、余分な分割が加えられていることになる。
【００１０】
そこで、本発明は、観測行列データを冗長さのない任意の長方形に分割することができる行列データ解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記の目的を達成するために本発明に係る行列データ解析装置は、入力された観測行列データ＾ＸＩと同じサイズの行列＾Ｉを少なくとも１つのブロックに分割したときの長方形分割нと、前記観測行列データ＾ＸＩを少なくとも１つのブロックに分割したときの長方形分割Ｒと、前記観測行列データ＾ＸＩにおける各行ｍが所属する各クラスタξｍと、前記観測行列データ＾ＸＩにおける各列ｎが所属する各クラスタηｎと、の各々に対して初期値を設定するパラメータ初期化部と、前記長方形分割нに基づいて、前記長方形分割Ｒの新しい候補Ｒｎｅｗを求め、前記観測行列データ＾ＸＩ、前記長方形分割н、前記クラスタξｍ、及び前記クラスタηｎに基づいて算出される前記長方形分割Ｒの事後確率と前記新しい候補Ｒｎｅｗの事後確率とを比較して、前記新しい候補Ｒｎｅｗを採択するか否かを決定して、前記長方形分割Ｒを更新する長方形分割変更部と、前記観測行列データ＾ＸＩ、前記長方形分割Ｒ、及び前記長方形分割нに基づいて、前記長方形分割нを更新する中間変数更新部と、前記観測行列データ＾ＸＩ、前記長方形分割Ｒ、及び前記長方形分割нに基づいて、前記観測行列データ＾ＸＩにおける各行ｍのクラスタξｍと各列ｎのクラスタηｎとを更新するクラスタ更新部と、予め定められた終了条件を満たすまで、前記長方形分割更新部による更新、前記中間変数更新部による更新、及び前記クラスタ更新部による更新を繰り返し行う終了判定部と、を含んで構成されている。
【００１２】
本発明に係る行列データ解析方法は、パラメータ初期化部によって、入力された観測行列データ＾ＸＩと同じサイズの行列＾Ｉを少なくとも１つのブロックに分割したときの長方形分割нと、前記観測行列データ＾ＸＩを少なくとも１つのブロックに分割したときの長方形分割Ｒと、前記観測行列データ＾ＸＩにおける各行ｍが所属する各クラスタξｍと、前記観測行列データ＾ＸＩにおける各列ｎが所属する各クラスタηｎと、の各々に対して初期値を設定し、長方形分割変更部によって、前記長方形分割нに基づいて、前記長方形分割Ｒの新しい候補Ｒｎｅｗを求め、前記観測行列データ＾ＸＩ、前記長方形分割н、前記クラスタξｍ、及び前記クラスタηｎに基づいて算出される前記長方形分割Ｒの事後確率と前記新しい候補Ｒｎｅｗの事後確率とを比較して、前記新しい候補Ｒｎｅｗを採択するか否かを決定して、前記長方形分割Ｒを更新し、中間変数更新部によって、前記観測行列データ＾ＸＩ、前記長方形分割Ｒ、及び前記長方形分割нに基づいて、前記長方形分割нを更新し、クラスタ更新部によって、前記観測行列データ＾ＸＩ、前記長方形分割Ｒ、及び前記長方形分割нに基づいて、前記観測行列データ＾ＸＩにおける各行ｍのクラスタξｍと各列ｎのクラスタηｎとを更新し、終了判定部によって、予め定められた終了条件を満たすまで、前記長方形分割更新部による更新、前記中間変数更新部による更新、及び前記クラスタ更新部による更新を繰り返す。
【００１３】
本発明に係るプログラムは、上記の行列データ解析装置の各部としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００１４】
以上説明したように、本発明の行列データ解析装置、方法、及びプログラムによれば、観測行列データを冗長さのない任意の長方形に分割することができる、という効果が得られる。
【図面の簡単な説明】
【００１５】
【図１】（Ａ）“regular grid”クラスの長方形分割の例、（Ｂ）“hierarchical”クラスの長方形分割の例、（Ｃ）“arbitrary”クラスの長方形分割の例を示す図である。
【図２】本発明の実施の形態に係る行列データ解析装置の構成を示す概略図である。
【図３】本発明の実施の形態に係る行列データ解析装置における行列データ解析処理ルーチンの内容を示すフローチャートである。
【図４】本発明の実施の形態に係る行列に含まれる全要素の全順序を確率的に生成するアルゴリズムの一例を示すフローチャートである。
【図５】本発明の実施の形態に係る長方形分割過程のアルゴリズムの一例を示すフローチャートである。
【図６】（Ａ）観測行列データを示す図、（Ｂ）本発明の実施の形態に係る行列データ解析処理ルーチンを実行して生成したサンプルの一例を示す図、（Ｃ）（Ｂ）のサンプルに対する対数事後確率の変化を示した図である。
【発明を実施するための形態】
【００１６】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１７】
＜発明の原理＞
【００１８】
まず、本発明の原理について説明する。
【００１９】
本発明はベイズ推論の枠組みに基づいている。ベイズ推論は一般に事前分布と尤度関数と呼ばれる二つの確率モデルを用いて対象をモデル化する。そのモデルを用い、対象とする行列データに確率モデルをフィッティングすることで解析を行う。
【００２０】
まず、本発明の実施の形態に係る事前分布について説明する。
【００２１】
今、行列＾Ｉの全ての長方形分割を＾ＴＩと表す。ただし、ベクトル、行列、又は集合を表す記号には、記号の前に「＾」を付して表現する。事前分布μＩとしては任意の長方形分割Ｒ∈＾ＴＩに対して確率μＩ（Ｒ）が計算でき、（１）式を満たすものであれば何でもよい。
【００２２】
【数１】
2015088117.tif 000003
【００２３】
このような事前分布を与える生成モデルとして、例えばギルバート模様を考えることができる。しかし、ギルバート模様は統計的な振る舞いの解析が極めて難しいことが知られており、実用上取扱いが簡単ではない。そこで、ここでは後述する長方形分割過程を用いた方法を示すものとする。
【００２４】
次に、本発明の実施の形態に係る尤度関数について説明する。
【００２５】
行列＾Ｉに対してある長方形分割Ｒが与えられたとき、尤度関数ｆＩを設定する必要があるが、これは観測行列データ＾ＸＩに応じて制限される。例えば、観測行列データ＾ＸＩがバイナリ変数であれば、ベータ分布−二項分布の階層モデルを用いることができる。
【００２６】
このように、観測行列データ＾ＸＩに応じて尤度関数ｆＩと事前分布μＩを設定することで解析を行うことができる。また、観測行列データ＾ＸＩが複数枚あった場合においても同種の方法が適用できる。ここでは一例として、一枚の観測行列データ＾ＸＩに自然数、尤度モデルとしてディリクレ分布−多項分布の階層モデル、事前分布として以下のような長方形分割過程を用いた場合の例を示す。
【００２７】
＜長方形分割過程＞
【００２８】
以下に、行列＾Ｉに対して潜在的に全ての長方形分割（図１（Ｃ）に示すような“arbitrary”と呼ばれる何の制約も持たない長方形分割）を生成しうる確率的生成モデルの乱択アルゴリズムを示す。
【００２９】
このアルゴリズムは実数ｐ∈（０、１）、及び実数Ｂ≧０を入力として、ある事前分布μＩに従って１つの長方形分割Ｒ∈＾ＴＩのサンプルを出力する乱択アルゴリズムに相当していると共に、ある長方形分割Ｒを用意したとき、当該長方形分割Ｒが出力され得る確率を計算するためのアルゴリズムとしても使用することができる。
【００３０】
アルゴリズムは、大きく分けて２つの段階（それぞれアルゴリズム１とアルゴリズム２という）から成る。最初の段階は行列＾Ｉの全ての要素の全順序θＩ：＝（＾Ｉ；＜~）を離散モンドリアン過程に基づいて確率的に生成する。ここで“＜~”及び“＞~”は順序を含めた大小関係を示す記号である。
【００３１】
具体的には、＾ΘＩで、行列＾Ｉの全ての要素の昇順の全順序（任意の１≦ｉ＜ｊ≦ (ｍ’−ｍ＋１)(ｎ'−ｎ＋１)に対してθＩ（ｉ）＜~θＩ（ｊ））を表現し、アルゴリズム１で全順序θＩ∈＾ΘＩを確率的に生成する。そして、アルゴリズム１により全順序θＩが得られた後、この全順序θＩに基づいて、行列＾Ｉの長方形分割の１つである長方形分割Ｒを生成する。
【００３２】
また、入力データの確率的生成モデルは以下のように表せる。まず、行列＾Ｉの各行ｍ＝１,２,・・，Ｍには［０,１］の一様分布から一つの実数が割り当てられており、［０,１］上の強度をμとしたポアソン過程により点によって分けられた［０,１］上の領域において、同一の領域に入った行は同一のクラスタとみなすことでクラスタリングを行う。このように確率的に生成されるｍ（ｍ＝１,２,・・,Ｍ）番目の行のクラスタξｍを（２）式のように表現する。
【００３３】
【数２】
2015088117.tif 000004
【００３４】
ここでは簡単のため、以下＾ξ＝（ξ１,ξ２,・・,ξＭ）と表すことにする。列に関しても同様に、各列ｎ＝１,２,・・,Ｎのクラスタηｎは（３）式のように表される。
【００３５】
【数３】
2015088117.tif 000005
【００３６】
なお、行と同様に＾η＝（η１,η２,・・,ηＮ）と表わすことにする。中間変数としての｛１,２,・・,Ｍ｝×｛１,２,・・,Ｎ｝の長方形分割нは（４）式及び（５）式のようにモンドリアン過程により生成されたものと表現される。
【００３７】
【数４】
2015088117.tif 000006
【００３８】
そして、モンドリアン過程、すなわちアルゴリズム１により生成された長方形分割нに条件付けを実施することで、長方形分割Ｒが生成される。（６）式はこの状況を定式化したものである。
【００３９】
【数５】
2015088117.tif 000007
【００４０】
そして、長方形分割Ｒができた後、長方形分割Ｒの各行列要素Ｘｍ、ｎは（７）式及び（８）式のように、離散値を取るディリクレ分布−離散分布の階層モデルにより生成されたとみなされる。
【００４１】
【数６】
2015088117.tif 000008
【００４２】
ここで＾αは非負の実数を持つ任意の次元数のベクトルであり、例えば＾αは観測行列データ＾ＸＩの要素の最大値に設定すればよい。この場合、ディリクレ分布−離散分布の関係は共役になっているためφｒ（ｒ∈Ｒ）は周辺化し消去することができる。
【００４３】
そして、（９）式の事後確率に対して、マルコフ連鎖モンテカルロ法により反復的に推論を行う。
【００４４】
【数７】
2015088117.tif 000009
【００４５】
（９）式の右辺の各項は、それぞれアルゴリズム２で示した長方形分割過程、アルゴリズム１で示したモンドリアン過程、ベルヌーイ分布、ディリクレ離散分布のディリクレ変数を周辺化したもの、ポアソン過程、ポアソン過程の確率密度関数を表している。
【００４６】
以下では長方形分割Ｒ、長方形分割н、バイナリ変数＾ｌ、クラスタ＾ξ、及びクラスタ＾ηを逐次的に反復演算することで（９）式の事後確率にフィットした長方形分割Ｒ、長方形分割н、バイナリ変数＾ｌ、クラスタ＾ξ、及びクラスタ＾ηを１組としたサンプルの生成を実施する。
【００４７】
＜システム構成＞
【００４８】
次に、観測行列データを長方形分割することによりクラスタを抽出する行列データ解析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【００４９】
図２に示すように、本発明の実施の形態に係る行列データ解析装置は、ＣＰＵと、ＲＡＭと、後述する行列データ解析処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【００５０】
入力部１０により、解析対象である観測行列データ＾ＸＩが入力される。記憶部３０は、入力部１０により入力された観測行列データ＾ＸＩを記憶する。また、記憶部３０は、後述する各処理での結果を記憶すると共に、本処理ルーチンで用いる各パラメータの初期値を記憶している。
【００５１】
演算部２０は、パラメータ初期化部２１と、長方形分割更新部２２と、中間変数更新部２３と、クラスタ更新部２４と、終了判定部２５とを備えている。
【００５２】
パラメータ初期化部２１は、例えば予め定めたベクトル＾α、実数ｐ∈（０、１）、及び観測行列データ＾ＸＩを入力として、後述する処理で用いる各パラメータである長方形分割н、バイナリ変数＾ｌ、長方形分割Ｒ、クラスタ＾ξ、及びクラスタ＾ηの各初期値を設定する。
【００５３】
具体的には、長方形分割нの初期値として、観測行列データ＾ＸＩと同じサイズの行列＾Ｉの“hierarchical”クラスに属する長方形を設定する。例えば、分割が一つもなく元の行列の全ての要素が同一クラスタに入った長方形分割を初期値として設定する。
【００５４】
また、バイナリ変数＾ｌの初期値として、長方形分割нのクラスタ数と同じ次元のバイナリベクトル＾ｌの初期値を設定する（例えば、＾ｌ＝０）。
【００５５】
また、長方形分割Ｒの初期値として、観測行列データ＾ＸＩと同じサイズの行列の長方形分割を設定する。例えば、分割が一つもなく元の行列の全ての要素が同一クラスタに入った長方形分割を初期値として設定する。
【００５６】
更に、クラスタ＾ξの初期値として、全ての行ｍが同一クラスタに属するように設定し、クラスタ＾ηの初期値として、全ての列ｎが同一クラスタに属するように設定する。
【００５７】
長方形分割更新部２２は、例えば予め定めたベクトル＾α、実数ｐ、及び観測行列データ＾ＸＩ、長方形分割н、バイナリ変数＾ｌ、長方形分割Ｒ、クラスタ＾ξ、及びクラスタ＾ηを入力として、メトロポリスヘイスティングス法により長方形分割Ｒを更新して、記憶部３０に記憶する。なお、長方形分割Ｒ以外の各入力パラメータは、長方形分割更新部２２において更新されることはない。
【００５８】
具体的には、記憶部３０に記憶されている現在の長方形分割н及び長方形分割нの各ブロックＨに対するバイナリ変数＾ｌＨを用いて、長方形分割過程に従って、新しい長方形分割候補Ｒｎｅｗを求める。なお、ここでブロックとは、長方形分割された各々の長方形を指すものとする。
【００５９】
また、メトロポリスヘイスティングス法の際に使用する採択確率は、例えば（１０）式に従うものとする。
【００６０】
【数８】
2015088117.tif 000010
【００６１】
なお、長方形分割過程以外の提案分布から、新しい長方形分割候補Ｒｎｅｗを求めてもよい。また、長方形分割нに対して１行だけ変更したり、１列だけ変更したりするようなものであってもよい。
【００６２】
中間変数更新部２３は、例えば予め定めたベクトル＾α、実数ｐ、及び観測行列データ＾ＸＩ、長方形分割н、バイナリ変数＾ｌ、長方形分割Ｒ、クラスタ＾ξ、及びクラスタ＾ηを入力として、長方形分割Ｒを生成する際の中間変数である長方形分割н及びバイナリ変数＾ｌを更新して、記憶部３０に記憶する。なお、長方形分割н及びバイナリ変数＾ｌ以外の各入力パラメータは、中間変数更新部２３において更新されることはない。
【００６３】
モンドリアン過程で示される長方形分割нの更新は、各反復毎に１つの分割を加えるか、１つの分割を取り除くかをメトロポリスヘイスティングスアルゴリズムを用いて決定するリヴァーシブルマルコフ連鎖モンテカルロ法に基づいて実施され、バイナリ変数＾ｌの更新は、ギブス更新により実施される。
【００６４】
なお、長方形分割нの更新の際、モンドリアンはcostＥｋと分割位置のχｋで構成されているため、これらを同時に更新するようにする。
【００６５】
今、ｋ番目のブロックがｋ’とｋ”とに分割されたとすると、モンドリアンに関係するパラメータζは＜ζ，Ｅｋ，χｋ＞となるとみなすことができる。一方で、分割が取り除かれる場合は、Ｅｋ及びχｋが＜ζ，Ｅｋ，χｋ＞から取り除かれζになる。リヴァーシブルマルコフ連鎖モンテカルロ法を適用するためには、現在のパラメータと提案のパラメータとの次元をそろえる必要があるため、分割を加える場合はζｉと補助変数とを加えた＜ζｉ，ｕ１，ｕ２＞を用いζｉ＋１＝＜ζｉ，Ｅｋ，χｋ＞とする。このとき、提案分布Ｑ（нｉ＋１；нｉ）は、χｋを一様分布から、Ｅｋを1/τｋ（τｋは非負の実数）を平均パラメータとする指数分布から更新する。
【００６６】
ｋ番目のブロックに分割を加えるか、分割を取り除くかの決定では、分割候補となるブロック、取り除かれる分割の候補から、一様分布に従って決定される。
【００６７】
ｋ番目のブロックから分割を取り除く場合は、Ｑ（нｉ＋１；нｉ）は、新しいパラメータζｉ＋１を確定的に決めることができる。
【００６８】
提案分布Ｑが分割を取り除こうとしているとき、採択確率ａをUniform（０，１）から生成し、（１１）式に示す条件を満たすときだけ採択する。
【００６９】
【数９】
2015088117.tif 000011
【００７０】
ただし、｜∂ζｉ＋１／∂ζｉ｜はヤコビアンを表している。一方、分割を加えようとしている場合は、（１２）式に示す条件を満たす場合だけ採択する。
【００７１】
【数１０】
2015088117.tif 000012
【００７２】
ただし、ｗｋはｋ番目のブロックが正の実数であるbudgetλを使い果たす確率であり、下記のように表わせる。
【００７３】
【数１１】
2015088117.tif 000013
【００７４】
分割を加える場合は下記のように表せる。
【００７５】
【数１２】
2015088117.tif 000014
【００７６】
一方、バイナリ変数＾ｌの更新に関しては、分割された各ブロックＨ（∈н）に関して、バイナリ変数＾ｌＨが０または１を取る条件付き確率を計算し、０を取る条件付き確率と１を取る条件付き確率との割合に従って、バイナリ変数ｌＨを更新する。
【００７７】
すなわち、（９）式をバイナリ変数ｌＨ＝１として計算した確率と、（９）式をバイナリ変数ｌＨ＝０として計算した確率との割合に従って、バイナリ変数ｌＨを更新する。
【００７８】
クラスタ更新部２４は、例えば予め定めたベクトル＾α、実数ｐ、及び観測行列データ＾ＸＩ、長方形分割н、バイナリ変数＾ｌ、長方形分割Ｒ、クラスタ＾ξ、及びクラスタ＾ηを入力として、クラスタ＾ξ及びクラスタ＾ηを更新して、記憶部３０に記憶する。なお、クラスタ＾ξ及びクラスタ＾η以外の各入力パラメータは、クラスタ更新部２４において更新されることはない。
【００７９】
クラスタ＾ξ及びクラスタ＾ηの更新は、バイナリ変数＾ｌの更新と同様にギブス更新により実施される。すなわち、各ｍ＝１,２,・・,Ｍに関してξｍが取り得る各値に対する条件付き確率を計算し、計算された各値に対する条件付き確率の割合に従ってξｍを更新する。同様に、各ｎ＝１,２,・・,Ｎに関してηｎが取り得る各値に対する条件付き確率を計算し、計算された各値に対する条件付き確率の割合に従ってηｎを更新する。
【００８０】
すなわち、（９）式を、ξｍが取り得る各値に対して計算した確率の割合に従って、ξｍを更新し、（９）式を、ηｎが取り得る各値に対して計算した確率の割合に従って、ηｎを更新する。
【００８１】
終了判定部２５は、予め定められた終了条件を満足するか否かを判定し、終了条件を満足していない場合には、長方形分割更新部２２〜クラスタ更新部２４の各処理を繰り返す。
【００８２】
本実施の形態では、生成するサンプル（長方形分割Ｒ、長方形分割н、バイナリ変数＾ｌ、クラスタ＾ξ、及びクラスタ＾η）の個数の上限を予め指定しておく。計算時間の許す範囲内で、サンプルの個数の上限にできるだけ大きな値を指定するほど、生成されたサンプルは事後分布のより良い近似となるが、実用的にはサンプルの個数の上限として数千個を指定するのが一般的である。
【００８３】
この際、生成されたサンプルが指定の個数に満たない場合は反復を繰り返す。一方、生成されたサンプルが指定の個数に達した場合には反復を終了し、今までに生成した全てのサンプルを記憶部３０に記憶する。
【００８４】
なお、終了条件として、繰り返し回数が予め定められた上限回数に到達したことを用いてもよい。
【００８５】
このように更新を繰り返すことで、そこから生成されたサンプルは事後分布をよく近似したものとなる。
【００８６】
また、生成された各サンプルのうち、（９）式に適用して得られる事後確率が最大となるサンプルが最適なサンプルとなるため、実際の解析結果を用いる際には、この事後確率が最大となるサンプルのみを記憶部３０に記憶するようにしてもよい。
【００８７】
＜行列データ解析装置の作用＞
【００８８】
次に、本実施の形態に係る行列データ解析装置１００の作用について説明する。まず、解析対象のデータとして、観測行列データ＾ＸＩが行列データ解析装置１００に入力され、記憶部３０に格納される。そして、行列データ解析装置１００において、図３に示す行列データ解析処理ルーチンが実行される。
【００８９】
まず、ステップＳ３０１において、記憶部３０から観測行列データ＾ＸＩを読み込む。
【００９０】
そして、ステップＳ３０２において、予め定められたベクトル＾α、実数ｐ、及びステップＳ３０１で読み込んだ観測行列データ＾ＸＩを用いて、長方形分割н、バイナリ変数＾ｌ、長方形分割Ｒ、クラスタ＾ξ、及びクラスタ＾ηの各初期値を設定し、記憶部３０に記憶する。また、非負の実数を値に持つ変数ｔに、初期値Ｔを設定する。
【００９１】
例えば、長方形分割нの初期値として、分割が一つもなく元の行列の全ての要素が同一クラスタに入った長方形分割を設定する。また、バイナリ変数＾ｌの初期値として、例えば、＾ｌ＝０に設定する。また、長方形分割Ｒの初期値として、例えば、分割が一つもなく元の行列の全ての要素が同一クラスタに入った長方形分割を設定する。更に、クラスタ＾ξの初期値として、全ての行ｍが同一クラスタに属するように設定し、クラスタ＾ηの初期値として、全ての列ｎが同一クラスタに属するように設定する。
【００９２】
次にステップＳ３０３では、観測行列データ＾ＸＩ、非負の実数を値に持つ変数ｔ、及びバイナリ変数＾ｌに基づいて、長方形分割нのブロックＨ毎に、図４に示したアルゴリズム１を実行することにより、観測行列データ＾ＸＩと同じサイズの行列＾Ｉの全ての要素の全順序θＩを生成して、記憶部３０に記憶する。
【００９３】
図４は、全順序θＩ∈＾ΘＩを確率的に生成するアルゴリズム１の一例を示すフローチャートである。
【００９４】
なお、アルゴリズム１で全順序θＩを確率的に生成する際、全ての自然数ｕ、ｖ＜ｖ’に対して、（ｕ,ｖ）＜~（ｕ,ｖ’）、及び、全ての自然数ｖ、ｕ＜ｕ’に対して（ｕ,ｖ）＜~（ｕ’,ｖ）を満たすような全順序θＩが、全順序θＩの初期値として設定されているものとする。
【００９５】
アルゴリズム１は、長方形分割нのブロックＨに対応する行列＾Ｘ、非負の変数ｔ、及び全順序θＩを入力とする関数SampleＭ(＾Ｘ,ｔ, θＩ)として表される。また、アルゴリズム１内部の説明では、入力された行列を行列＾Ｘ、入力された非負の変数を変数ｔとして表すものとする。なお、入力された全順序θＩはそのまま全順序θＩとして表わす。また、初期化後のアルゴリズム１は、関数SampleＭ(＾Ｉ,Ｔ, θＩ)として表される。
【００９６】
まず、ステップＳ１０１では、ξ〜Exponential(ｃ({＾Ｘ}))として、指数分布に従う乱数ξを生成する。ただし、ｃ(・)は行列＾Ｘの行数及び列数の総数を表す演算子である。
【００９７】
ステップＳ１０２では、変数ｔとステップＳ１０１で生成した乱数ξの値を比較し、変数ｔが乱数ξ未満か否かを判定する。否定判定の場合はステップＳ１０７へ移行し、肯定判定の場合はステップＳ１０３へ移行する。
【００９８】
ステップＳ１０３では、行列＾Ｘに対応するブロックＨに対する現在のバイナリ変数＾ｌHが０か否かを判定し、肯定判定の場合はステップＳ１０４へ移行し、否定判定の場合はステップＳ１０５へ移行する。
【００９９】
ステップＳ１０４では、ｕ,ｖを自然数とした場合、｛(ｕ,ｖ),(ｕ',ｖ')｜(ｕ,ｖ),(ｕ',ｖ')∈＾Ｘ,ｖ＜ｖ'｝に対して(ｕ,ｖ)＜~(ｕ',ｖ')となるように、全順序θＩを更新する。一方、ステップＳ１０５では、｛(ｕ,ｖ),(ｕ',ｖ')｜(ｕ,ｖ),(ｕ',ｖ')∈＾Ｘ,ｕ＜ｕ'｝に対して(ｕ,ｖ)＜~(ｕ',ｖ')となるように、全順序θＩを更新する。
【０１００】
ステップＳ１０６では、ステップＳ１０４又はステップＳ１０５で更新した全順序θＩ及び行列＾Ｘをアルゴリズム１の読み出し元へ戻して、本アルゴリズムを終了する。
【０１０１】
一方、ステップＳ１０２で変数ｔが乱数ξ以上であると判定された場合、ステップＳ１０７へ移行し、一様分布に従って行列＾Ｘを水平方向又は垂直方向に分割して、行列＾Ｘ０及び行列＾Ｘ１を生成する。この際、行列＾Ｘ０及び行列＾Ｘ１に対して、（１９）式で示されるバイナリ変数ｌｘを割り当てる。
【０１０２】
【数１３】
2015088117.tif 000015
【０１０３】
ステップＳ１０８では、ｘ０＜~ｘ１を満たすｘ０∈＾Ｘ０、ｘ１∈＾Ｘ１が存在するか否かを判定する。そして、肯定判定の場合にはステップＳ１０９へ移行し、否定判定の場合にはステップＳ１１０へ移行する。
【０１０４】
ステップＳ１０９では、全てのｘ０∈＾Ｘ０、ｘ１∈＾Ｘ１に対してｘ０＜~ｘ１となるように、全順序θＩを更新する。一方、ステップＳ１１０では、全てのｘ０∈＾Ｘ０、ｘ１∈＾Ｘ１に対してｘ０＞~ｘ１となるように、全順序θＩを更新する。
【０１０５】
ステップＳ１１１では、ステップＳ１０７で生成した行列＾Ｘ０及び行列＾Ｘ１の各々について、行列＾Ｘ０又は行列＾Ｘ１を新たな行列＾Ｘ、変数ｔをステップＳ１０１で生成した乱数ξで減算した値（ｔ-ξ）を新たな変数ｔ、及びステップＳ１０９又はステップＳ１１０で更新した全順序θＩを新たな全順序θＩとして、SampleＭ(＾Ｘｉ,t-ξ, θＩ)を呼び出すことで、再帰的に、本アルゴリズムを実行し、本アルゴリズムを終了する。
【０１０６】
なお、この再帰的に呼び出されたSampleＭ(＾Ｘｉ,t-ξ, θＩ)の処理では、ステップＳ１０７で行列＾Ｘ０及び行列＾Ｘ１に割り当てたベクトル変数ｌXに従って、ステップＳ１０３の判定を実施するものとする。
【０１０７】
このように、全順序θＩはアルゴリズム１に示すSampleＭ(＾Ｘ,ｔ,θＩ)に従って再帰的に更新される。
【０１０８】
次のステップＳ３０４では、観測行列データ＾ＸＩと同じサイズの行列＾Ｉ、実数ｐ、及びステップＳ３０３で生成した全順序θＩに基づいて、図５に示したアルゴリズム２を実行して、行列＾Ｉの長方形分割を長方形分割候補Ｒｎｅｗとして生成し、記憶部３０に記憶する。
【０１０９】
図５は、行列＾Ｉから１つの長方形分割を生成するアルゴリズム２：SampleＲ(＾Ｉ,θＩ,ｐ)の一例を示すフローチャートである。なお、図５の説明において、θＩ（ｉ）は、全順序θＩの行列＾Ｉのｉ番目の要素を表すものとし、θＩ（１）が表す要素を出発点とする。
【０１１０】
まず、ステップＳ２０１では、ｉを２に設定する。次にステップＳ２０２では、２つの隣接した要素θＩ（ｉ）＋（−１,０）及び要素θＩ（ｉ）＋（０,−１）（以下隣接要素という）が行列＾Ｉに存在し、且つ、２つの隣接要素が既にブロックに割り当てられているか否か（条件１）を判定する。そして、否定判定の場合にはステップＳ２０８へ移行し、肯定判定の場合にはステップＳ２０３へ移行する。
【０１１１】
ステップＳ２０３では、２つの隣接要素が同一ブロックに割り当てられているか否かを判定し、肯定判定の場合にはステップＳ２０４へ移行し、否定判定の場合にはステップＳ２０５へ移行する。
【０１１２】
ステップＳ２０４では、２つの隣接要素が同一ブロックに割り当てられていることから、要素θＩ（ｉ）も２つの隣接要素と同一のブロックに割り当てるようにする（ブロック割り当て処理１）。
【０１１３】
一方、ステップＳ２０３で否定判定となった際の移行先であるステップＳ２０５では、更に要素θＩ（ｉ）の全順序に関して、θＩ（ｉ）＋（−１,０）＜~θＩ（ｉ）＋（１,０）を満たすか否かを判定する。そして、肯定判定の場合にはステップＳ２０６へ移行し、否定判定の場合にはステップＳ２０７へ移行する。
【０１１４】
ステップＳ２０６では、確率ｐで要素θＩ（ｉ）をθＩ（ｉ）＋（−１,０）の要素と同一ブロックに割り当て、確率（１−ｐ）ｑで要素θＩ（ｉ）をθＩ（ｉ）＋（０,−１）の要素と同一ブロックに割り当てる（ブロック割り当て処理２）。ここでｑは、（２０）式によって示される、実数ｐに基づいて決定される値である。
【０１１５】
【数１４】
2015088117.tif 000016
【０１１６】
また、ステップＳ２０７では、確率ｐで要素θＩ（ｉ）をθＩ（ｉ）＋（０,−１）の要素と同一ブロックに割り当て、確率（１−ｐ）ｑで要素θＩ（ｉ）をθＩ（ｉ）＋（−１,０）の要素と同一ブロックに割り当て、確率（１−ｐ）・（１−ｑ）で要素θＩ（ｉ）を新しいブロックに割り当てる（ブロック割り当て処理３）。
【０１１７】
また、ステップＳ２０８では、確率ｐで要素θＩ（ｉ）を隣接要素と同一ブロックへ割り当て、確率（１−ｐ）で要素θＩ（ｉ）を新しいブロックに割り当てる（ブロック割り当て処理４）。
【０１１８】
このようにして要素θＩ（ｉ）をブロックに割り当てた後、ステップＳ２０９で要素番号ｉを１増加し、ステップＳ２１０で、行列＾Ｉに含まれる全ての要素θＩ（ｉ）をブロックに割り当てたか否かを判定する。具体的には、要素番号ｉが (ｍ’−ｍ＋１)(ｎ'−ｎ＋１)以下か否かを判定する。ここで、ｍ及びｍ'は行列＾Ｉの行を示す値であり、ｎ及びｍ'は行列＾Ｉの列を示す値である。
【０１１９】
そして、否定判定の場合には、上記のように行列＾Ｉに含まれる要素θＩ（ｉ）をブロックに割り当てることで生成した最終的な長方形分割を、アルゴリズム２の読み出し元へ戻して本アルゴリズムを終了する。一方、肯定判定の場合には、まだ何れかのブロックに割り当てていない要素θＩ（ｉ）が存在することからステップＳ２０２へ移行し、次の要素θＩ（ｉ）を何れかのブロックに割り当てる処理を繰り返す。
【０１２０】
ステップＳ３０５では、中間変数である長方形分割нとバイナリ変数＾ｌ、クラスタ＾ξ、クラスタ＾η、行列＾Ｘ、実数ｐ、ベクトル＾α、及び実数λに基づいて長方形分割Ｒの事後確率を算出すると共に、ステップＳ３０４で生成した長方形分割候補Ｒｎｅｗの事後確率を算出する。そして、（１０）式に示す確率で、長方形分割候補Ｒｎｅｗを採択するか否かを決定し、長方形分割候補Ｒｎｅｗを採択すると決定した場合、長方形分割候補Ｒｎｅｗを新たな長方形分割Ｒとして更新し、記憶部３０に記憶する。
【０１２１】
ステップＳ３０６では、長方形分割нのｋ番目のブロックを更に分割するか、又は長方形分割нの隣接する２つのブロックを１つのブロックに統合するかを決定して分割候補нｎｅｗを生成すると共に、長方形分割нの事後確率及び分割候補нｎｅｗの事後確率を算出する。
【０１２２】
そして、長方形分割нのブロックを統合しようとする場合、すなわち分割を取り除こうとする場合には（１１）式を満足するか否かを判定し、（１１）式を満足する場合に分割候補нｎｅｗを新たな長方形分割нとして更新し、記憶部３０に記憶する。
【０１２３】
一方、長方形分割нを更に分割しようとする場合には、（１２）式を満足するか否かを判定し、（１２）式を満足する場合に分割候補нｎｅｗを新たな長方形分割нとして更新し、記憶部３０に記憶する。なお、（１１）式、（１２）式におけるнｉ＋１が分割候補нｎｅｗ、нｉが更新前の長方形分割нである。
【０１２４】
ステップＳ３０７では、長方形分割Ｒ、長方形分割н、バイナリ変数＾ｌ、クラスタ＾ξ、クラスタ＾η、行列＾Ｘに基づいて、上記（９）式に従って、長方形分割нの各ブロックＨに対するバイナリ変数であるｌＨが０または１を取る条件付き確率を計算し、０を取る条件付き確率と１を取る条件付き確率との割合に従ってバイナリ変数ｌＨを更新し、記憶部３０に記憶する。
【０１２５】
ステップＳ３０８では、長方形分割Ｒ、長方形分割н、バイナリ変数＾ｌ、クラスタ＾ξ、クラスタ＾η、行列＾Ｘに基づいて、上記（９）式に従って、各行ｍ＝１,２,・・,Ｍに関してξｍが取り得る各値に対する条件付き確率を計算し、各値に対する条件付き確率の割合に従ってξｍを更新し、記憶部３０に記憶する。
【０１２６】
ステップＳ３０９では、長方形分割Ｒ、長方形分割н、バイナリ変数＾ｌ、クラスタ＾ξ、クラスタ＾η、行列＾Ｘに基づいて、上記（９）式に従って、各列ｎ＝１,２,・・,Ｎに関してηｎが取り得る各値に対する条件付き確率を計算し、各値に対する条件付き確率の割合に従ってηｎを更新し、記憶部３０に記憶する。
【０１２７】
なお、上記ステップＳ３０５で更新された長方形分割Ｒと、上記ステップＳ３０６で更新された長方形分割нと、上記ステップＳ３０７で更新されたバイナリ変数＾ｌと、上記ステップＳ３０８で更新されたクラスタ＾ξと、上記ステップＳ３０９で更新されたクラスタ＾ηとは、行列データ解析処理ルーチンの実行により生成されたサンプルとして、互いに関連付けられた形で記憶部３０に記憶される。
【０１２８】
次のステップＳ３１０では、上記ステップＳ３０５〜ステップＳ３０９によって生成されたサンプル数が、予め記憶部３０に記憶されているサンプルの個数の上限値よりも小さいか否かを判定し、現在の生成サンプル数が上限値未満の場合には、終了条件を満足していないと判断して、記憶部３０に記憶されている現在の生成サンプル数を１増加させると共に、上記ステップＳ３０３へ戻り、上記ステップＳ３０３〜ステップＳ３０９の処理を繰り返し、次のサンプルを生成する。
【０１２９】
一方、現在の生成サンプル数が上限値以上の場合には、終了条件を満足したと判断して、サンプルの生成を終了する。
【０１３０】
そして、ステップＳ３１１では、生成されたサンプルの各々を出力部４０から出力して、行列データ解析処理ルーチンを終了する。
【０１３１】
＜実施結果＞
【０１３２】
次に、本実施の形態に係る手法の有用性を示す目的で、Jester Collaborative Filtering Dataset(http://goldberg.berkeley.edu/jester-data/)の行列データに対して、図３に示した行列データ解析処理ルーチンを実行し、長方形分割によるクラスタの抽出を行った評価実験の結果について説明する。
【０１３３】
図６は、提案法を用いた場合の評価実験結果を示した図である。図６（Ａ）は入力した観測行列データを示す。また、図６（Ｂ）は、マルコフ連鎖モンテカルロ法により１０００個のサンプルを生成した際に、事後確率が最大となるサンプルを示したものであり、図６（Ｃ）は、マルコフ連鎖モンテカルロ法における対数事後確率の変化を示したものである。
【０１３４】
図６から、観測行列データに対して、“regular grid”クラスにも“hierarchical”クラスにも属さない、“arbitrary”と呼ばれる何の制約も持たない長方形分割が実現されていることがわかる。
【０１３５】
以上説明したように、本発明の実施の形態に係る行列データ解析装置によれば、観測行列データに対して長方形分割Ｒと、長方形分割нと、バイナリ変数＾ｌと、クラスタ＾ξと、クラスタ＾ηとを繰り返し更新することにより、観測行列データが制約のない長方形に分割され、観測行列データに対するクラスタ抽出精度を向上させることができる。
【０１３６】
また、従来の制限のある長方形分割方法に、中間変数としての長方形分割нとバイナリ変数＾ｌとの更新処理を付加し、長方形分割を行うことにより、“regular grid”に限定されていない“arbitrary”と呼ばれる全ての長方形分割を潜在的に生成することができる。
【０１３７】
また、従来の制限のある長方形分割で用いられる確率分布に、逐次的な確率計算を可能とする長方形分割過程を用いることにより、計算速度の高速化を可能にする。
【０１３８】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１３９】
例えば、上述の行列データ解析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０１４０】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【０１４１】
１０  入力部
２０  演算部
２１  パラメータ初期化部
２２  長方形分割更新部
２３  中間変数更新部
２４  クラスタ更新部
２５  終了判定部
３０  記憶部
４０  出力部
１００行列データ解析装置
【図１】
2015088117.tif 000017
【図２】
2015088117.tif 000018
【図３】
2015088117.tif 000019
【図４】
2015088117.tif 000020
【図５】
2015088117.tif 000021
【図６】
2015088117.tif 000022