(57)【要約】      （修正有）【課題】画像や音響信号のスペクトログラムのような非負値（０以上の値）からなるデータから、ポアソン分布とガウス分布の混合分布により定義される関数を尤度関数とし、適応的に観測データの統計分布の分布形を推定し、ＮＭＦ（非負値行列分解）を行う。【解決手段】ポアソン分布とガウス分布の基底行列、混合比および重み行列の初期値を生成し、観測データに基づいたガウス分布の分散より生成されたガウス分布率行列で混合比率を更新し、基底行列および重み行列の更新を収束するまで繰り返す。【選択図】図１
2010039723.tif 000002
【特許請求の範囲】
【請求項１】
非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算する非負値行列分解の数値計算方法において、
初期パラメータ生成部が、
前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成し、
前記ガウス分布率行列生成部が、
前記初期パラメータ生成部から入力された前記パラメータの初期値、あるいは、パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散の値に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成し、
前記パラメータ更新部の分布混合比更新部が、
前記ガウス分布率行列に基づき、前記分布混合比を更新し、
前記パラメータ更新部のガウス分布分散更新部が、
前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新し、
前記パラメータ更新部の基底行列更新部が、
前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新し、
前記パラメータ更新部の重み行列更新部が、
前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新し、
収束判定部は、
前記パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かを判定し、
前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成部に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記パラメータ更新部により更新されたパラメータの更新を再度判定し、
分解行列出力部は、
前記収束判定部により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、更新された前記基底行列および前記重み行列を出力する
ことを特徴とする非負値行列分解の数値計算方法。
【請求項２】
前記収束判定部は、前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、前記ポアソン分布とガウス分布の混合分布により定義された目的関数を最大化する所定の基準を満たしているか否かを判断することを特徴とする請求項１に記載の非負値行列分解の数値計算方法。
【請求項３】
非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算する非負値行列分解の数値計算装置において、
前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成する初期パラメータ生成部と、
前記初期パラメータ生成部から入力された前記パラメータの初期値、あるいは、パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成するガウス分布率行列生成部と、
前記ガウス分布率行列生成部から入力された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散を、それぞれ更新する基底行列更新部、重み行列更新部、分布混合比更新部、およびガウス分布分散更新部を有するパラメータ更新部と、
前記パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かを判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成部に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記パラメータ更新部により更新されたパラメータの更新を再度判定する収束判定部と、
前記収束判定部により、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、前記基底行列および前記重み行列を出力する分解行列出力部とを備え、
前記パラメータ更新部は、
前記ガウス分布率行列に基づき、前記分布混合比を更新する分布混合比更新部と、
前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新するガウス分布分散更新部と、
前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新する基底行列更新部と、
前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新する重み行列更新部を有することを特徴とする非負値行列分解の数値計算装置。
【請求項４】
前記収束判定部は、前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、前記ポアソン分布とガウス分布の混合分布により定義された目的関数を最大化する所定の基準を満たしているか否かを判断することを特徴とする請求項３に記載の非負値行列分解の数値計算装置。
【請求項５】
非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算するコンピュータに、
前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成する初期パラメータ生成手段と、
前記初期パラメータ生成手段から入力された前記パラメータの初期値、あるいは、更新手段により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成するガウス分布率行列生成手段と、
前記ガウス分布率行列に基づき、前記分布混合比を更新する更新手段と、
前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新する更新手段と、
前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新するガウス分布分散更新手段と、
前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新する基底行列更新手段と、
前記更新手段により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かの前記パラメータの更新を判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成手段に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記更新手段により更新された前記パラメータの更新を再度判定する収束判定手段と、
前記収束判定手段により、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、前記基底行列および前記重み行列を出力する分解行列出力手段と、
を実行させるためのプログラム。
【請求項６】
非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算するコンピュータに、
前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成する初期パラメータ生成手段と、
前記初期パラメータ生成手段から入力された前記パラメータの初期値、あるいは、更新手段により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成するガウス分布率行列生成手段と、
前記ガウス分布率行列に基づき、前記分布混合比を更新する更新手段と、
前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新する更新手段と、
前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新する更新手段と、
前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新する更新手段と、
前記更新手段により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かの前記パラメータの更新を判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成手段に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記更新手段により更新された前記パラメータの更新を再度判定する収束判定手段と、
前記収束判定手段により、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、前記基底行列および前記重み行列を出力する分解行列出力手段とを実行させるためのプログラムを記憶した記憶媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、主に、画像の特徴抽出や音源分離で用いられる技術に関するものであって、特に、画像や音響信号のスペクトログラムのような非負値からなるデータに基づき、ＮＭＦ（非負値行列分解）モデルのパラメータを計算する非負値行列分解の数値計算方法、非負値行列分解の数値計算装置、プログラムおよび記憶媒体に関する。
【背景技術】
【０００２】
非負値行列分解（Ｎｏｎ−ｎｅｇａｔｉｖｅ  Ｍａｔｒｉｘ  Ｆａｃｔｏｒｉｚａｔｉｏｎ；ＮＭＦ）は、画像や音響信号のスペクトログラムなどの非負値データから特徴的なパターンを抽出する方法として有効であり、画像特徴抽出や音源分離など多岐にわたって応用されている。この非負値行列とは、全成分が非負（０以上）の値で構成される行列であって、非負値行列分解とは、観測データ（非負値ベクトル）のすべてを非負値基底ベクトルの非負結合でよく近似できるような適当な基底系を構成する理論である。
言い換えると、非負値行列分解は、すべての観測データ（Ｊ個のＫ次元非負値ベクトル）
【数１】
2010039723.tif 000003
を、Ｉ個のＫ次元非負値基底ベクトル
【数２】
2010039723.tif 000004
の非負結合
【数３】
2010039723.tif 000005
で良く近似できるような適当な基底系を構成する理論である。
【０００３】
すなわち、Ｊ個のＫ次元非負値ベクトルを並べた非負値行列Ｘ＝（ｘｋ，ｊ）Ｋ×Ｊ＝（ｘ１，ｘ２，・・・ｘＪ）、Ｉ個のＫ次元非負値基底ベクトルを並べた基底行列Ｗ＝（ｗｋ，ｉ）Ｋ×Ｉ＝（ｗ１，ｗ２，・・・ｗＩ）、非負結合の係数を要素とした重み行列Ｈ＝（ｈｉ，ｊ）Ｉ×Ｊとした場合、
【数４】
2010039723.tif 000006
となるような非負値行列である基底行列Ｗ，重み行列Ｈを求めるものが非負値行列分解（ＮＭＦ）である。
【０００４】
ここで、これらの近さの尺度を、フロベニウス（Ｆｒｏｂｅｎｉｕｓ）ノルム
【数５】
2010039723.tif 000007
と、Ｉダイバージェンス（またはＫＬダイバージェンス）
【数６】
2010039723.tif 000008
で定義した場合の最適な基底行列Ｗおよび重み行列Ｈを数値計算するための効率的な反復アルゴリズムが、開示されている（非特許文献１）。これらは乗法更新（Ｍｕｌｔｉｐｌｉｃａｔｉｖｅ  Ｕｐｄａｔｅ）アルゴリズムと呼ばれている。なお、式６において、［・］ｋ，ｊは、「・」の行列の（ｋ，ｊ）成分を表す。
【０００５】
ＷＴをＷの転置行列とすると、フロベニウスノルムで定義された式５に示すＦ（Ｘ，ＷＨ）を小さくする更新式は、それぞれ
【数７】
2010039723.tif 000009
【数８】
2010039723.tif 000010
で与えられ、Ｉダイバージェンスで定義された式６に示すＩ（Ｘ，ＷＨ）を小さくする更新式はそれぞれ
【数９】
2010039723.tif 000011
【数１０】
2010039723.tif 000012
で与えられる。
【０００６】
いずれの更新式も、基底行列Ｗと重み行列Ｈの初期値が非負値行列であれば必ず非負値行列に更新されることが保証されている点に特徴があり、目的関数の定義によって異なる更新式が導入されるため、当然のことながら目的関数の定義によって異なる値に収束する。
したがって、基底行列Ｗおよび重み行列Ｈの行列同士の近さをフロベニウスノルムで測った場合で最適解を得ることは、目的関数を式５に示したフロベニウスノルムとして定義し、尤度関数を次式に示すガウス（Ｇａｕｓｓ）分布と仮定して最尤推定していることに相当する。
【数１１】
2010039723.tif 000013
【０００７】
また、基底行列Ｗおよび重み行列Ｈの行列同士の近さをＩダイバージェンスで測った場合で最適解を得ることは、目的関数を式６に示したＩダイバージェンスとして定義し、尤度関数を次式に示すポアソン（Ｐｏｉｓｓｏｎ）分布と仮定して最尤推定していることに相当する。
【数１２】
2010039723.tif 000014
ただし、観測データｘｋ，ｊは量子化された整数値データ
【数１３】
2010039723.tif 000015
とする。
【０００８】
すなわち、観測データが実際にどのような統計分布に従って生成されたかに従って、式５に示したフロベニウスノルムあるいは式６に示したＩダイバージェンスとしての目的関数のうち、最小にするべき目的関数がいずれであるかが決まるわけである。
しかしながら、一般に、観測データがどのような統計分布に従って生成されたかを事前に知ることはできず、かつ、必ずしもガウス分布やポアソン分布に理想的に従っているとは限らない。このため、従来、状況に応じてどちらかの目的関数を採用することが好ましいかを経験的に選択していた。
【非特許文献１】Ｄ．Ｄ．Ｌｅｅ  ａｎｄ  Ｈ．Ｓ．Ｓｅｕｎｇ，“Ｌｅａｒｎｉｎｇ  ｔｈｅ  Ｐａｒｔｓ  ｏｆ  Ｏｂｊｅｃｔｓ  ｂｙ  Ｎｏｎ−ｎｅｇａｔｉｖｅ  Ｍａｔｒｉｘ  Ｆａｃｔｏｒｉｚａｔｉｏｎ，”Ｎａｔｕｒａ  Ｖｏｌ．４０１，ｐｐ．７８８−７９１，１９９９．
【発明の開示】
【発明が解決しようとする課題】
【０００９】
しかしながら、後で詳細に説明するが、図３に示す分布図からわかるように、ポアソン分布の尤度関数の場合、尤度関数の分布の裾は、観測値ｘが大きいほど広くなる。このため、観測値ｘ（観測行列の要素の値）が小さいほど、ＮＭＦモデル（基底行列Ｗと重み行列Ｈの積をとった行列）の要素値と観測値ｘとのずれに対して感度が鋭くなるようなコスト（尤度関数）を設けていることに相当する。一方、図４に示す分布図からわかるように、ガウス分布の尤度関数の場合、尤度関数の分布の裾は、観測値ｘの大小にかかわらず一定であるため、観測値ｘが大きいほど、ＮＭＦモデルの要素値と観測値ｘとのずれに対して感度が鋭くなるようなコストを設けていることに相当する。
【００１０】
つまり、尤度関数をポアソン分布と仮定すると、小さな観測値ｘの微妙な誤差に影響を受けてＮＭＦモデルのパラメータの推定精度が低下しやすくなる。逆に、尤度関数をガウス分布と仮定すると、観測値ｘに含まれる外れ値や飛び値の影響を受けて、ＮＭＦモデルのパラメータの推定精度が低下しやすくなる。
また、音響信号のスペクトログラムのようなスパースなデータ（鋭いピーク成分がまばらにあり、それ以外は極端に小さい成分からなるデータ）に関しては、上記のいずれの影響も陽に受けることになる。したがって、尤度関数としてポアソン分布と仮定した場合は、小さい成分の影響を大きく受け、尤度関数としてガウス分布と仮定した場合は、鋭いピーク成分（外れ値や飛び値）の影響を大きく受けることになる。このため、どちらの尤度関数と仮定したとしても、ＮＭＦモデルのパラメータの推定精度が低下してしまうという問題がある。
【００１１】
本発明は、このような事情を考慮し、上記の問題を解決すべくなされたものであって、その目的は、ポアソン分布とガウス分布の混合分布により定義される関数を尤度関数とし、観測データに応じてＮＭＦパラメータだけでなく、観測データの統計分布の分布形も適応的に推定するための非負値行列分解の数値計算方法、非負値行列分解の数値計算装置、プログラムおよび記憶媒体を提供することにある。
【課題を解決するための手段】
【００１２】
上記問題を解決するために、本発明は、非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算する非負値行列分解の数値計算方法において、初期パラメータ生成部が、前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成し、前記ガウス分布率行列生成部が、前記初期パラメータ生成部から入力された前記パラメータの初期値、あるいは、パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散の値に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成し、前記パラメータ更新部の分布混合比更新部が、前記ガウス分布率行列に基づき、前記分布混合比を更新し、前記パラメータ更新部のガウス分布分散更新部が、
前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新し、前記パラメータ更新部の基底行列更新部が、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新し、前記パラメータ更新部の重み行列更新部が、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新し、収束判定部は、前記パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かを判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成部に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記パラメータ更新部により更新されたパラメータの更新を再度判定し、分解行列出力部は、前記収束判定部により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、更新された前記基底行列および前記重み行列を出力することを特徴とする非負値行列分解の数値計算方法である。
【００１３】
また、本発明にかかる収束判定部は、前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、前記ポアソン分布とガウス分布の混合分布により定義された目的関数を最大化する所定の基準を満たしているか否かを判断することを特徴とする非負値行列分解の数値計算方法である。
【００１４】
また、本発明は、非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算する非負値行列分解の数値計算装置において、前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成する初期パラメータ生成部と、前記初期パラメータ生成部から入力された前記パラメータの初期値、あるいは、パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成するガウス分布率行列生成部と、前記ガウス分布率行列生成部から入力された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散を、それぞれ更新する基底行列更新部、重み行列更新部、分布混合比更新部、およびガウス分布分散更新部を有するパラメータ更新部と、前記パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かを判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成部に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記パラメータ更新部により更新されたパラメータの更新を再度判定する収束判定部と、前記収束判定部により、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、前記基底行列および前記重み行列を出力する分解行列出力部とを備え、前記パラメータ更新部は、前記ガウス分布率行列に基づき、前記分布混合比を更新する分布混合比更新部と、前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新するガウス分布分散更新部と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新する基底行列更新部と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新する重み行列更新部を有することを特徴とする非負値行列分解の数値計算装置である。
【００１５】
また、本発明にかかる収束判定部は、前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、前記ポアソン分布とガウス分布の混合分布により定義された目的関数を最大化する所定の基準を満たしているか否かを判断することを特徴とする非負値行列分解の数値計算装置である。
【００１６】
また、本発明は、非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算するコンピュータに、前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成する初期パラメータ生成手段と、前記初期パラメータ生成手段から入力された前記パラメータの初期値、あるいは、更新手段により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成するガウス分布率行列生成手段と、前記ガウス分布率行列に基づき、前記分布混合比を更新する更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新する更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新するガウス分布分散更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新する基底行列更新手段と、前記更新手段により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かの前記パラメータの更新を判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成手段に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記更新手段により更新された前記パラメータの更新を再度判定する収束判定手段と、前記収束判定手段により、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、前記基底行列および前記重み行列を出力する分解行列出力手段と、を実行させるためのプログラムである。
【００１７】
また、本発明は、非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算するコンピュータに、前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成する初期パラメータ生成手段と、前記初期パラメータ生成手段から入力された前記パラメータの初期値、あるいは、更新手段により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成するガウス分布率行列生成手段と、前記ガウス分布率行列に基づき、前記分布混合比を更新する更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新する更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新する更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新する更新手段と、前記更新手段により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かの前記パラメータの更新を判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成手段に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記更新手段により更新された前記パラメータの更新を再度判定する収束判定手段と、前記収束判定手段により、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、前記基底行列および前記重み行列を出力する分解行列出力手段とを実行させるためのプログラムを記憶した記憶媒体である。
【発明の効果】
【００１８】
本発明によれば、観測データが、ポアソン分布ないしガウス分布のいずれかに従っている場合、どちらの分布に従っているかが推定されるため、そのもとで適切なＮＭＦパラメータを推定することができる。
また、本発明によれば、観測データがポアソン分布ないしガウス分布のいずれかに従っていない場合であっても、観測データの実際の統計分布にポアソン分布とガウス分布の混合分布で表現可能な範囲内で柔軟に尤度関数が適応されるため、そのもとで最尤のＮＭＦパラメータを得ることができる。
さらに、本発明によれば、効率的にＮＭＦパラメータを算出可能な乗法更新アルゴリズムを実現することができる。
また、本発明によれば、音響信号のスペクトログラムのようなスパースなデータに対しても、高精度なＮＭＦパラメータを推定することができる。
【発明を実施するための最良の形態】
【００１９】
以下、図面を参照して、本発明の一実施形態について説明する。図１は、本発明の実施の形態の構成を示すブロック図であり、図２は、その動作を説明するフローチャートである。
図１に示すとおり、本実施の形態の非負値行列分解の計算装置は、初期パラメータ生成部１、ガウス分布率行列生成部２、パラメータ更新部３、収束判定部４および分解行列出力部５を有する。
【００２０】
初期パラメータ生成部１は、非負値行列分解（ＮＭＦ）モデルのパラメータである基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２の初期値を生成し、ガウス分布率行列生成部２に出力する。例えば、初期パラメータ生成部１は、観測データｘに基づき、所定の値が設定されたパラメータの初期値を記憶し、また、設定変更・修正が可能な構成を有する。
ただし、基底行列Ｗと重み行列Ｈは、それぞれＫ×Ｉ行列、Ｉ×Ｊ行列である。なお、ここで言うＫ，Ｊは、観測データ行列Ｘの行数および列数であり、Ｉは、基底ベクトルが備える非負値基底ベクトルの個数であって予め定められている任意の定数である。
それぞれの初期値の選択可能な範囲は、基底行列Ｗ＝（ｗｋ，ｉ）Ｋ×Ｉに関しては、すべてのｋとｉに対しての範囲として、
【数１４】
2010039723.tif 000016
また、重み行列Ｈ＝（ｈｉ，ｊ）Ｉ×Ｊに関しては、すべてのｉとｊに対しての範囲として、
【数１５】
2010039723.tif 000017
分布混合比λに関しては、
【数１６】
2010039723.tif 000018
ガウス分布の分散σ２に関しては、σ２＞０とする。
なお、基底行列Ｗ＝（ｗｋ，ｉ）Ｋ×Ｉは、Ｋ×Ｉ行列で（ｋ，ｉ）成分がｗｋ，ｉである行列（ｗ１，ｗ２，・・・ｗＩ）を表す。また、分布混合比λは、ガウス分布とポアソン分布の混合比率を表す。
【００２１】
ガウス分布率行列生成部２は、入力信号として入力された観測データ行列Ｘに基づき、初期パラメータ生成部１から入力された、あるいはパラメータ更新部３から読み出した基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２を用いて、次式に基づきｋ，ｊすべてについてガウス分布率γｋ，ｉを算出する。
【数１７】
2010039723.tif 000019
なお、観測データ行列Ｘは、所定の観測データの出力部（図示せず）により、非負値データを格納した非負値ベクトルが並べられた観測データ行列である。なお、観測データｘｋ，ｊは、観測データ行列Ｘのｋ行ｊ行列目の成分、［・］ｋ，ｊは行列「・」のｋ行ｊ列目の成分を表す。
ガウス分布率行列生成部２は、上述の式１７を利用して計算されたガウス分布率γｋ，ｊを要素とした、次式に示すガウス分布率行列Γを生成し、パラメータ更新部３に出力する。
【数１８】
2010039723.tif 000020
なお、ガウス分布率行列Γとは、観測データ行列Ｘの各要素ｘがガウス分布とポアソン分布のうち、ガウス分布から生成された確率を表したものである。
また、ガウス分布率行列生成部２は、入力された観測データ行列Ｘ、および初期パラメータ生成部１あるいは分解行列出力部５から入力された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２を、パラメータ更新部３に出力する。
【００２２】
パラメータ更新部３は、分布混合比更新部３１、ガウス分布分散更新部３２、基底行列更新部３３、重み行列更新部３４を含み、これらにより更新された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２を収束判定部４に出力する。
【００２３】
分布混合比更新部３１は、ガウス分布率行列生成部２から入力されたガウス分布率行列Γの要素であるガウス分布率γｋ，ｊをガウス分布率γ´ｋ，ｊとし、これを用いて次式に基づき分布混合比λを更新する。
【数１９】
2010039723.tif 000021
分布混合比更新部３１は、更新した分布混合比λを収束判定部４に出力する。
【００２４】
ガウス分布分散更新部３２は、ガウス分布率行列生成部２から入力されたガウス分布の分散σ２を次式に基づき更新する。つまり、ガウス分布分散更新部３２は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布率γｋ，ｊを、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布率γ´ｋ，ｊとし、これらと観測データ行列Ｘ＝（ｘｋ，ｊ）Ｋ×Ｊを用いて、次式に基づきガウス分布の分散σ２を更新する。
【数２０】
2010039723.tif 000022
ガウス分布分散更新部３２は、更新したガウス分布の分散σ２を収束判定部４に出力する。
【００２５】
基底行列更新部３３は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布行列Γ、およびガウス分布の分散σ２を、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布行列Γ´、ガウス分布の分散σ´２とし、これらと観測データ行列Ｘ＝（ｘｋ，ｊ）Ｋ×Ｊを用いて、次式に基づき基底行列Ｗを更新する。
【数２１】
2010039723.tif 000023
ただし、行列Ａ，Ｂ，Ｃ（太字）は、それぞれ次式と定義される。
【数２２】
2010039723.tif 000024
行列の除法は、行列要素同士で商をとる演算とし、○中に中黒「・」で示す記号
【数２３】
2010039723.tif 000025
は、Ｈａｄａｍａｒｄ積（行列の要素同士の積をとる演算）を、√内に「・」で示す表記は、「・」が行列のとき、行列の要素ごとの平方根をとる演算を表す。また、「１（太字）」は、要素がすべて１のＫ×Ｊの行列を表す。
基底行列更新部３３は、更新した基底行列Ｗを収束判定部４に出力する。
【００２６】
重み行列更新部３４は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布行列Γ、およびガウス分布の分散σ２を、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布行列Γ´、ガウス分布の分散σ´２とし、これらと観測データ行列Ｘ＝（ｘｋ，ｊ）Ｋ×Ｊを用いて、次式に基づき重み行列Ｈを更新する。
【数２４】
2010039723.tif 000026
ただし、行列Ｄ，Ｅ，Ｆ（太字）は、それぞれ次式で定義される。
【数２５】
2010039723.tif 000027
重み行列更新部３４は、更新した重み行列Ｈを収束判定部４に出力する。
【００２７】
収束判定部４は、パラメータ更新部３により入力された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２が、それぞれ所定の基準を満たしているか否かのパラメータの更新を判定する。収束判定部４は、基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２が、それぞれ所定の基準を満たしていない場合、ガウス分布率行列生成部２にガウス分布率行列の生成を指示する信号を出力し、それぞれのパラメータが所定の基準を満たすまで、この生成されたガウス分布率行列Γに基づきパラメータ更新部３により更新されたパラメータの更新の判定を繰り返す。
本実施の形態において、収束判定部４は、パラメータ更新部３により、基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２が、下に示す式２６のポアソン分布とガウス分布の混合分布により定義された目的関数を最大化するよう更新され、この目的関数の値の変化率が所定値以下となったか否かを判定する。
【数２６】
2010039723.tif 000028
収束判定部４は、式２６に示す目的関数の変化率が所定値以下となった場合、基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２が、所定の基準を満たしていると判定し、更新された最新の基底行列Ｗおよび重み行列Ｈを出力する指示の信号を、分解行列出力部５に出力する。
【００２８】
なお、収束判定部４は、上述の判定方法に限られず、例えば、パラメータの更新計算を反復して行った回数が、所定の回数に達したか否かを検出し、所定に回数に達した場合、分解行列出力部５に指示の信号を出力する構成であってもよい。また、収束判定部４は、パラメータの反復更新計算において、それぞれのパラメータの更新の変化率が所定の値以下であるか否かを検出し、変化率が所定値以下になった場合、分解行列出力部５に指示の信号を出力する構成であってもよい。
【００２９】
分解行列出力部５は、収束判定部４から出力の指示信号を受けて、更新された最新の基底行列Ｗおよび重み行列Ｈを出力する。
【００３０】
この構成により、本発明は、ポアソン分布とガウス分布の混合分布により定義される尤度関数を、観測データｘに応じて、ＮＭＦパラメータだけでなく、観測データｘの統計分布の分布形も適応的に推定することができる。
つまり、観測データｘがポアソン分布ないしガウス分布のいずれかに従っている場合、どちらの分布に従っているか収束判定部４により判定されるまで、パラメータ更新部３による更新が繰り返される。このため、観測データｘがどちらの分布に従っているかが推定され、同時にそのもとで適切なＮＭＦパラメータが推定される。
一方、観測データｘがポアソン分布ないしガウス分布のいずれかに従っていない場合であっても、観測データｘの実際の統計分布にポアソン分布とガウス分布の混合分布で表現可能な範囲内で柔軟に尤度関数が適応されるため、そのもとで最尤のＮＭＦパラメータを得ることができる。
したがって、本発明は、観測データｘが、いずれの統計分布の分布形に従っているか不明な場合であっても、観測データｘに応じた最尤のＮＭＦパラメータを得ることができる。
【００３１】
さらに、本発明によれば、従来経験的に観測データｘの統計分布の分布形を判断していたところ、自動的に最適な統計分布に応じた最尤のＮＭＦパラメータを計算することができるため、効率的にＮＭＦパラメータを計算可能な乗法更新アルゴリズムを実現することができる。
また、本発明によれば、音響信号のスペクトログラムのようなスパースなデータに対しても、高精度なＮＭＦパラメータを推定することができる。
【００３２】
次に、図２を用いて、本発明にかかる非負値行列分解の数値計算方法について説明する。図２は、本実施の形態における非負値行列分解の数値計算方法の一例を示すフローチャートである。
【００３３】
図２に示すとおり、まず、初期パラメータ生成部１により、基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２の初期値が生成される（Ｓ１）。生成された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２の初期値は、ガウス分布率行列生成部２に出力される。
ガウス分布率行列生成部２は、入力された観測データ行列Ｘに基づき、これら基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２の初期値を用いて、式１７によりガウス分布率γｋ，ｉを計算し、このガウス分布率γｋ，ｉよりガウス分布率行列Γを生成する（Ｓ２）。ガウス分布率行列生成部２は、生成したガウス分布率行列Γと、観測データ行列Ｘ、基底行列Ｗ、重み行列Ｈ、分布混合比λおよびガウス分布の分散σ２の初期値を、パラメータ更新部３に出力する。
【００３４】
パラメータ更新部３は、入力された情報に基づき、ＮＭＦパラメータの更新を行う。
すなわち、分布混合比更新部３１は、ガウス分布率行列生成部２から入力されたガウス分布率行列Γの要素であるガウス分布率γｋ，ｊをガウス分布率γ´ｋ，ｊとし、これを用いて式１９に基づき分布混合比λを更新する（Ｓ３）。
ガウス分布分散更新部３２は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布率γｋ，ｊを、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布率γ´ｋ，ｊとし、これらと観測データ行列Ｘ＝（ｘｋ，ｊ）Ｋ×Ｊを用いて、式２０に基づきガウス分布の分散σ２を更新する（Ｓ４）。
基底行列更新部３３は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布行列Γ、およびガウス分布の分散σ２を、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布行列Γ´、ガウス分布の分散σ´２とし、これらと観測データ行列Ｘ＝（ｘｋ，ｊ）Ｋ×Ｊを用いて、式２１に基づき基底行列Ｗを更新する（Ｓ５）。
重み行列更新部３４は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布行列Γ、およびガウス分布の分散σ２を、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布行列Γ´、ガウス分布の分散σ´２とし、これらと観測データ行列Ｘ＝（ｘｋ，ｊ）Ｋ×Ｊを用いて、式２４に基づき重み行列Ｈを更新する（Ｓ６）。
【００３５】
収束判定部４は、入力された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２が、それぞれ所定の基準を満たしているか否かを判定する（Ｓ７）。
基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２が、それぞれ所定の基準を満たしている場合（Ｓ７−ＹＥＳ）、収束判定部４は、更新された最新の基底行列Ｗおよび重み行列Ｈを出力する指示の信号を、分解行列出力５に出力する。
分解行列出力部５は、収束判定部４から出力の指示信号を受けて、更新された最新の基底行列Ｗおよび重み行列Ｈを出力する（Ｓ８）。
【００３６】
一方、基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２が、それぞれ所定の基準を満たしていない場合（Ｓ７−ＮＯ）、収束判定部４は、ガウス分布率行列生成部２にガウス分布率行列の生成を指示する信号を出力する。
ステップＳ２に戻って、指示信号を受けたガウス分布率行列生成部２は、パラメータ更新部３から更新された最新の基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２を読み出し、観測データｘに基づきガウス分布率行列Γを生成し、パラメータ更新部３に出力する。パラメータ更新部３は、ガウス分布率行列生成部２から入力されたガウス分布率行列をおよび観測データｘに基づき基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２を更新し、収束判定部４に出力する。これら観測データｘおよびパラメータ等が入力された収束判定部４は、パラメータの更新の判定を行う。
すなわち、収束判定部４は、それぞれのパラメータが所定の基準を満たすまで、この生成されたガウス分布率行列Γに基づきパラメータ更新部３により更新されたパラメータの更新の判定を繰り返す。
【００３７】
なお、パラメータ更新部３によるこれら更新の順番は任意であり、例えば、基底行列Ｗおよび重み行列Ｈの更新を先に実行し、その後に分布混合比λ、ガウス分布の分散σ２の更新を実行する構成であってもよく、全ての更新を同時に実行する構成であってもよい。
パラメータ更新部３は、更新された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２を収束判定部４に出力する。
【００３８】
［観測データに基づく統計分布設定の問題］
ここで、ガウス分布とポアソン分布の尤度関数について説明する。
図３は、ポアソン分布の尤度関数の概形を示す図であり、図４は、ガウス分布の尤度関数の概形を示す図であり、図５は、ポアソン分布とガウス分布の混合分布の概形を示す図である。図３〜５において、凸型の複数の波形ａ〜ｏは、それぞれ異なる観測値ｘの尤度関数を示している。
【００３９】
図３に示すとおり、ポアソン分布の尤度関数である複数の波形ａ〜ｅは、波形ａの観測データｘが最も小さく、アルファベット順に観測データｘが大きなり、それぞれ観測データｘ＝１，１０，３０，６０，９０であるときのポアソン分布の尤度関数を示す。図３に示すポアソン分布の概形は、横軸に観測データｘ、縦軸にポアソン分布の尤度関数（式８８参照）を規定する。
図３からわかるように、尤度関数がポアソン分布の場合、観測データｘが大きいほど尤度関数の分布の裾が広くなる。このため、本実施の形態では、観測データｘの値が大きいほど、ＮＭＦモデル（基底行列Ｗと重み行列Ｈの積をとった行列）［ＷＨ］ｋ，ｊの要素値と観測データｘｋ，ｊとのずれに対して感度が鈍くなるようなコスト（尤度関数）を設けていることに相当する。
【００４０】
一方、図４に示すとおり、ガウス分布の尤度関数である複数の波形ｆ〜ｊは、波形ｆの観測データｘが最も小さく、アルファベット順に観測データｘが大きくなり、それぞれ観測データｘ＝１，１０，３０，６０，９０であるときのガウス分布の尤度関数を示す。図４に示すガウス分布の概形は、横軸に観測値、縦軸にガウス分布の尤度関数（式８９参照）を規定する。
図４からわかるように、尤度関数がガウス分布の場合、ポアソン分布の場合とは逆に、尤度関数の分布の裾は観測データｘの大きさによらず一定である。このため、本実施の形態では、ポアソン分布のときに比べて観測値ｘｋ，ｊの値が大きいほどＮＭＦモデル［ＷＨ］ｋ，ｊの要素値と観測値ｘとのずれに対して感度が鋭くなるようなコストを設けていることに相当する。
【００４１】
したがって、尤度関数がポアソン分布の場合のほうが、小さい観測データｘｋ，ｊにおける微量な誤差によって、推定結果全体が影響をうけやすく、一方、尤度関数がガウス分布の場合のほうが、観測データｘの外れ値や飛び値がある場合の影響を受けやすい。
特に音響信号のスペクトログラムのようにスパースなデータ（鋭いピーク成分がまばらにあり、それ以外は極端に小さい成分からなるデータ）に関しては、上記のいずれの影響も陽に受けることになるため、いずれか一方の尤度関数を仮定した場合、必然的に仮定した尤度関数の影響を受けることとなってしまう。
【００４２】
つまり、ＮＭＦの性能を向上させるためには、観測データが従う統計分布により即した目的関数をもとに、ＮＭＦを行うことが望ましい。そのためには、観測データに応じた統計分布を推定しながら、ＮＭＦパラメータである基底行列Ｗ，重み行列Ｈを推定する必要がある。
本発明は、ポアソン分布とガウス分布の混合分布により尤度関数を定義し、観測データに応じてＮＭＦパラメータだけでなく、観測データの統計分布の分布形も適応的に推定することができる。よって本発明は、これら統計分布と尤度関数とのバランスを観測データから適応的に調節し、適切な感度のコストを設けることが可能であり、いかなる応用場面でも最適化基準の選択に迫られることなく、統一的に利用可能なＮＭＦの枠組みを確率する可能性を発揮することができる。
【００４３】
［ハイブリッド尤度関数の定義］
次に、図５を参照して、ポアソン分布とガウス分布の混合分布について説明する。
図５に示すポアソン分布とガウス分布の混合分布は、
【数２７】
2010039723.tif 000029
と定義される。ただし、ガウス分布は、
【数２８】
2010039723.tif 000030
ポアソン分布は、
【数２９】
2010039723.tif 000031
とする。
この式２７の混合分布の尤度関数は、式５に示すフロベニウスノルムおよび式６に示すＩダイバージェンスで定義される目的関数の両特性を含むハイブリッド尤度関数である。
ただし、ガウス分布とポアソン分布の分布混合比λは、０≦λ≦１である。すなわち、分布混合比λが大きいほど、混合分布はガウス分布に近似し、分布混合比λが小さいほど、混合分布はポアソン分布に近似する。
【００４４】
したがって、対数尤度関数Ｌは、
【数３０】
2010039723.tif 000032
となり、この対数尤度関数Ｌが最大となるような基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ２を算出することが好ましい。なお、この分布混合比λとガウス分布の分散σ２を自由パラメータとして扱うことで、観測データｘが従う確率分布をできるだけよくあらわすことができるようになる。つまり、統計分布と尤度関数とのバランスを観測データから適応的に調節し、適切な感度のコストを設けることが可能となる。
なお、本実施の形態において、パラメータ更新部３は、上述の対数尤度関数Ｌを目的関数として、この目的関数を最大化するように各パラメータを更新し、収束判定部４は、この目的関数の値の変化率が所定値以下となったか否かを判定している。
【００４５】
［ＮＭＦパラメータと分布パラメータの最尤推定］
次に、本実施の形態における乗法更新アルゴリズムを導くための、最適な補助関数の計算について説明する。
まず、対数関数の凸性に基づき、
【数３１】
2010039723.tif 000033
という条件もとで、次のような不等式が成り立つ。
【数３２】
2010039723.tif 000034
ただし、ガウス分布率行列Γは次式のとおりである。
【数３３】
2010039723.tif 000035
また、
【数３４】
2010039723.tif 000036
かつ
【数３５】
2010039723.tif 000037
という条件のもとでは、次のような不等式が成り立つ。
【数３６】
2010039723.tif 000038
さらに、
【数３７】
2010039723.tif 000039
かつ
【数３８】
2010039723.tif 000040
という条件のもとでは、次のような不等式が成り立つ。
【数３９】
2010039723.tif 000041
【００４６】
そこで、上述の不等式（式３２，３６，３９）において、いずれも
【数４０】
2010039723.tif 000042
である場合に等号が成立するため、
【数４１】
2010039723.tif 000043
が成り立つ。
だたし、３階のテンソルＭは
【数４２】
2010039723.tif 000044
とする。
【００４７】
以上より、次のような不等式が成り立つ。
【数４３】
2010039723.tif 000045
【数４４】
2010039723.tif 000046
従って、ガウス分布率γｋ，ｊが次のような場合、上述の不等式（式４３）が成り立つ。
【数４５】
2010039723.tif 000047
また、ｍi，k，jが次のような場合、上述の不等式（式４４）が成り立つ。
【数４６】
2010039723.tif 000048
【００４８】
このように、式４１，４４に基づき算出された関数Ｇ（Ｗ，Ｈ，λ，σ２，Γ，Ｍ）は、対数尤度関数Ｌ（Ｗ，Ｈ，λ，σ２）を効率的に局所最大化するための補助関数であり、式４５，４６による更新と、基底行列Ｗ、重み行列Ｈ、分布混合比λ、およびガウス分布の分散σ２の更新を繰り返すことで、対数尤度関数Ｌ（Ｗ，Ｈ，λ，σ２）を単調増加させることができる。
本実施の形態においては、上述の理論に基づき、ガウス分布率行列生成部２によるガウス分布率行列Γの更新、およびパラメータ更新部３による基底行列Ｗ、重み行列Ｈ、分布混合比λ、およびガウス分布の分散σ２の更新を繰り返すことで、本発明にかかる非負値行列分解の数値計算装置は、対数尤度関数Ｌを単調増加させている。
【００４９】
［乗法更新式の導入］
式４１により算出された補助関数Ｇに基づき算出される各パラメータの更新式について説明する。
まず、ガウス分布率行列Γ＝Γ´という条件のもとで、補助関数Ｇ（Ｗ，Ｈ，λ，σ２，Γ，Ｍ）が最大となる分布混合比λを求める。
【数４７】
2010039723.tif 000049
を解くと、
【数４８】
2010039723.tif 000050
となる。
なお、ガウス分布率γｋ，ｉの式４５より
【数４９】
2010039723.tif 000051
であるから
【数５０】
2010039723.tif 000052
となる。
【００５０】
次に、次式に基づき、基底行列Ｗ＝Ｗ´、重み行列Ｈ＝Ｈ´、ガウス分布率行列Γ＝Γ´、３階のテンソルＭ＝Ｍ´という条件のもとで、補助関数Ｇ（Ｗ，Ｈ，λ，σ２，Γ，Ｍ）が最大となるガウス分布の分散σ２を求める。
【数５１】
2010039723.tif 000053
この式５１を、補助関数Ｇ（Ｗ，Ｈ，λ，σ２，Γ，Ｍ）に代入して、
【数５２】
2010039723.tif 000054
を解くと、
【数５３】
2010039723.tif 000055
となる。
【００５１】
次に、次式に基づき、重み行列Ｈ＝Ｈ´、ガウス分布率行列Γ＝Γ´、３階のテンソルＭ＝Ｍ´、ガウス分布の分散σ２＝σ´２という条件のもとで、補助関数Ｇ（Ｗ，Ｈ，λ，σ２，Γ，Ｍ）が最大となる基底行列Ｗを求める。
【数５４】
2010039723.tif 000056
この式５４に
【数５５】
2010039723.tif 000057
を代入すると、
【数５６】
2010039723.tif 000058
となる。ただし、行列の除法は、要素同士で商をとる演算とし、○中に中黒「・」で示す記号（数２３参照）は、Ｈａｄａｍａｒｄ積（行列の要素同士の積をとる演算）を、√内に「・」で示す表記は、「・」が行列のとき、行列の要素ごとの平方根をとる演算とする。なお、式５６において、「１（太字）」は要素がすべて１のＫ×Ｊの行列である。
この式５６の両辺に、
【数５７】
2010039723.tif 000059
を乗じると、次に示す２次方程式となる。
【数５８】
2010039723.tif 000060
ここで、
【数５９】
2010039723.tif 000061
【数６０】
2010039723.tif 000062
【数６１】
2010039723.tif 000063
とすると、式５８の方程式の解は、
【数６２】
2010039723.tif 000064
となる。
ただし、
【数６３】
2010039723.tif 000065
のとき、すなわち、式４３の√の中が行列のとき、行列の要素ごとの平方根をとる演算子を表す。
以上により、基底行列Ｗの更新式は、
【数６４】
2010039723.tif 000066
となる。式６２の解は、非負値行列であるため、式６４の基底行列Ｗ´が非負値行列であれば基底行列Ｗは必ず非負値行列に更新される。
【００５２】
次に、次式に基づき、基底行列Ｗ＝Ｗ´、ガウス分布率行列Γ＝Γ´、３階のテンソルＭ＝Ｍ´、ガウス分布の分散σ２＝σ´２という条件のもとで、補助関数Ｇ（Ｗ，Ｈ，λ，σ２，Γ，Ｍ）が最大となる重み行列Ｈを求める。
【数６５】
2010039723.tif 000067
この式６５に、式５５を代入すると、
【数６６】
2010039723.tif 000068
となる。
この式６６の両辺に、
【数６７】
2010039723.tif 000069
を乗じると、次に示す２次方程式となる。
【数６８】
2010039723.tif 000070
【００５３】
ここで、
【数６９】
2010039723.tif 000071
【数７０】
2010039723.tif 000072
【数７１】
2010039723.tif 000073
とすると、式６８の方程式の解は、
【数７２】
2010039723.tif 000074
となる。
以上により、重み行列Ｈの更新式は、
【数７３】
2010039723.tif 000075
となる。基底行列Ｗと同様に、式７２の解は非負値行列であるため、式７３の重み行列Ｈ´が非負値行列であれば重み行列Ｈは必ず非負値行列に更新される。
本実施の形態において、パラメータ更新部３は、上述のとおり算出された式４８，５３，６４，７３に基づき、基底行列Ｗ、重み行列Ｈ、分布混合比λおよびガウス分布の分散σ２を更新し、ガウス分布率行列生成部２により生成されたガウス分布率行列Γのもとで補助関数Ｇを最大化させている。つまり、この補助関数Ｇを最大化させるように更新を繰り返すことによって、収束判定部４は、目的関数Ｌを効率的に局所最大化させるための基底行列Ｗ、重み行列Ｈ、分布混合比λおよびガウス分布の分散σ２を得て、パラメータの更新判定を実行することができる。
【００５４】
［更新式の性質］
ここで、下に示す式７４および式７５のときの、
【数７４】
2010039723.tif 000076
【数７５】
2010039723.tif 000077
前記基底行列Ｗおよび重み行列Ｈの更新式（式６４，７３参照）の性質について説明する。
式７４のとき、ガウス分布率γｋ，ｊの式４５より
【数７６】
2010039723.tif 000078
すなわち、ガウス分布率行列Γは、
【数７７】
2010039723.tif 000079
である。
同様に、式７５のとき、ガウス分布率行列Γは、
【数７８】
2010039723.tif 000080
である。
【００５５】
まず、ガウス分布率行列Γ´＝「１（太字）」すなわち、ガウス分布率行列Γ´がすべて１のＫ×Ｊの行列であるとき、
【数７９】
2010039723.tif 000081
となり、これらを基底行列Ｗおよび重み行列Ｈの更新式（式６４，７３参照）に代入すると、基底行列Ｗおよび重み行列Ｈは、それぞれ
【数８０】
2010039723.tif 000082
【数８１】
2010039723.tif 000083
となり、この結果は、フロベニウスノルムで定義されたＦ（Ｘ，ＷＨ）を小さくする基底行列Ｗおよび重み行列Ｈの更新式（式７，８参照）とそれぞれ一致する。
【００５６】
一方、
【数８２】
2010039723.tif 000084
という条件のもとで、次式が成り立つ。
【数８３】
2010039723.tif 000085
ただし、｜・｜は、「・」が行列のとき、行列の要素ごとの絶対値をとる演算子を表す。
【００５７】
同様にして
【数８４】
2010039723.tif 000086
では、
【数８５】
2010039723.tif 000087
となり、これらを基底行列Ｗおよび重み行列Ｈの更新式（式６４，７３参照）に代入すると、基底行列Ｗおよび重み行列Ｈは、それぞれ
【数８６】
2010039723.tif 000088
【数８７】
2010039723.tif 000089
となり、この結果は、Ｉダイバージェンスで定義されたＩ（Ｘ，ＷＨ）を小さくする基底行列Ｗおよび重み行列Ｈの更新式（式９，１０参照）とそれぞれ一致する。
【００５８】
したがって、基底行列Ｗの更新式（式６４）および重み行列Ｈの更新式（式７３）は、２種類の乗法更新式を内包するハイブリッドな乗法更新式である。
本発明は、パラメータ更新部３により、このハイブリッドな乗法更新式を用いて基底行列Ｗおよび重み行列Ｈを更新している。
このため、本発明は、ポアソン分布とガウス分布の混合分布により尤度関数を定義し、観測データに応じてＮＭＦパラメータだけでなく、観測データの統計分布の分布形も適用的に推定することができる。
【００５９】
なお、上述の図３に示したポアソン分布の尤度関数の概形は、次式に示すポアソン分布の尤度関数の概形である。
【数８８】
2010039723.tif 000090
また、上述の図４に示したガウス分布の尤度関数の概形は、次式に示すガウス分布の尤度関数の概形である。
【数８９】
2010039723.tif 000091
【００６０】
ここで、以下上述の数式についての記号の定義について説明する。
非負値行列とは、全成分が非負（０以上）の値で構成される行列である。
次式に示す記号は、「全てのｋ，ｉに対して」という意味を表している。
【数９０】
2010039723.tif 000092
次式に示す記号は、「行列Ｈ´の転置行列」という意味を表している。
【数９１】
2010039723.tif 000093
次式に示す記号は、「左辺の式を右辺の式で定義する」という意味を表している。
【数９２】
2010039723.tif 000094
次式に示す記号は、「Ｋ次元非負値ベクトルの空間」という意味を表している。
【数９３】
2010039723.tif 000095
次式に示す記号は、「ベクトルｗの全ての成分が非負（０以上）である」という意味を表している。
【数９４】
2010039723.tif 000096
次式に示す記号は、「左辺の式を右変の式で近似する」という意味を表している。
【数９５】
2010039723.tif 000097
式９６に示す記号は、式９７に示す「フロベニウスノルム」を表している。
【数９６】
2010039723.tif 000098
【数９７】
2010039723.tif 000099
次式に示す記号は、「右辺を左辺に代入する（左辺の式を右辺の式で置き換える）」という意味を表している。
【数９８】
2010039723.tif 000100
次式に示す記号は、「整数全体からなる集合」という意味を表している。
【数９９】
2010039723.tif 000101
次式に示す記号は、「０から１の間（区間）に含まれる（０＜ガウス分布率γｋ，ｊ＜１）」という意味を表している。
【数１００】
2010039723.tif 000102
【００６１】
なお、上述の初期パラメータ生成部１、ガウス分布率行列生成部２、パラメータ更新部３、収束判定部４、分解行列出力部５における、図２を用いて説明したそれぞれの動作の過程は、コンピュータに実行させるためのプログラムや、このプログラムとしてコンピュータ読み取り可能な記録媒体として利用可能であり、コンピュータシステムが読み出して実行することによって、上記処理が行われる。なお、ここでいう「コンピュータシステム」とは、ＣＰＵ及び各種メモリやＯＳ、周辺機器等のハードウェアを含むものである。
また、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、フラッシュメモリ等の書き込み可能な不揮発性メモリ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。
【００６２】
さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムが送信された場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリ（例えばＤＲＡＭ（Dynamic Random Access Memory））のように、一定時間プログラムを保持しているものも含むものとする。
また、上記プログラムは、このプログラムを記憶装置等に格納したコンピュータシステムから、伝送媒体を介して、あるいは、伝送媒体中の伝送波により他のコンピュータシステムに伝送されてもよい。ここで、プログラムを伝送する「伝送媒体」は、インターネット等のネットワーク（通信網）や電話回線等の通信回線（通信線）のように情報を伝送する機能を有する媒体のことをいう。
また、上記プログラムは、前述した機能の一部を実現するためのものであっても良い。さらに、前述した機能をコンピュータシステムに既に記録されているプログラムとの組合せで実現できるもの、いわゆる差分ファイル（差分プログラム）であっても良い。
【図面の簡単な説明】
【００６３】
【図１】本実施形態による非負値行列分解の数値計算装置の一例を示す概略ブロック図である。
【図２】本実施形態による非負値行列分解の数値計算方法の一例を示すフローチャートである。
【図３】ポアソン分布の尤度関数の概形である。
【図４】ガウス分布の尤度関数の概形である。
【図５】ポアソン分布とガウス分布の混合分布の概形である。
【符号の説明】
【００６４】
１  初期パラメータ生成部
２  ガウス分布率行列生成部
３  パラメータ更新部
４  収束判定部
５  分解行列出力部
３１  分布混合比更新部
３２  ガウス分布分散更新部
３３  基底行列更新部
３４  重み行列更新部
【図１】
2010039723.tif 000103
【図２】
2010039723.tif 000104
【図３】
2010039723.tif 000105
【図４】
2010039723.tif 000106
【図５】
2010039723.tif 000107