JP5068228B2

JP5068228B2 - Non-negative matrix decomposition numerical calculation method, non-negative matrix decomposition numerical calculation apparatus, program, and storage medium

Info

Publication number: JP5068228B2
Application number: JP2008201274A
Authority: JP
Inventors: 弘和亀岡; 邦夫柏野
Original assignee: Nippon Telegraph and Telephone Corp
Current assignee: Nippon Telegraph and Telephone Corp
Priority date: 2008-08-04
Filing date: 2008-08-04
Publication date: 2012-11-07
Anticipated expiration: 2028-08-04
Also published as: JP2010039723A

Description

本発明は、主に、画像の特徴抽出や音源分離で用いられる技術に関するものであって、特に、画像や音響信号のスペクトログラムのような非負値からなるデータに基づき、ＮＭＦ（非負値行列分解）モデルのパラメータを計算する非負値行列分解の数値計算方法、非負値行列分解の数値計算装置、プログラムおよび記憶媒体に関する。 The present invention mainly relates to techniques used in image feature extraction and sound source separation, and in particular, NMF (non-negative matrix decomposition) based on data consisting of non-negative values such as spectrograms of images and sound signals. The present invention relates to a numerical calculation method for non-negative matrix decomposition, a numerical calculation apparatus for non-negative matrix decomposition, a program, and a storage medium.

非負値行列分解（Ｎｏｎ−ｎｅｇａｔｉｖｅＭａｔｒｉｘＦａｃｔｏｒｉｚａｔｉｏｎ；ＮＭＦ）は、画像や音響信号のスペクトログラムなどの非負値データから特徴的なパターンを抽出する方法として有効であり、画像特徴抽出や音源分離など多岐にわたって応用されている。この非負値行列とは、全成分が非負（０以上）の値で構成される行列であって、非負値行列分解とは、観測データ（非負値ベクトル）のすべてを非負値基底ベクトルの非負結合でよく近似できるような適当な基底系を構成する理論である。
言い換えると、非負値行列分解は、すべての観測データ（Ｊ個のＫ次元非負値ベクトル）

を、Ｉ個のＫ次元非負値基底ベクトル

の非負結合

で良く近似できるような適当な基底系を構成する理論である。 Non-negative matrix decomposition (NMF) is effective as a method for extracting characteristic patterns from non-negative data such as spectrograms of images and sound signals, and has a wide range of applications such as image feature extraction and sound source separation. Has been. This non-negative matrix is a matrix composed of all non-negative (0 or more) values. Non-negative matrix decomposition is a non-negative combination of all non-negative basis vectors for all observed data (non-negative vector). Is a theory that constructs an appropriate basis set that can be approximated by
In other words, non-negative matrix decomposition is performed on all observed data (J K-dimensional non-negative vectors).

, I K-dimensional non-negative basis vectors

Non-negative coupling of

This is a theory that constructs an appropriate basis set that can be approximated by

すなわち、Ｊ個のＫ次元非負値ベクトルを並べた非負値行列Ｘ＝（ｘ_ｋ，_ｊ）_Ｋ×Ｊ＝（ｘ_１，ｘ_２，・・・ｘ_Ｊ）、Ｉ個のＫ次元非負値基底ベクトルを並べた基底行列Ｗ＝（ｗ_ｋ，_ｉ）_Ｋ×Ｉ＝（ｗ_１，ｗ_２，・・・ｗ_Ｉ）、非負結合の係数を要素とした重み行列Ｈ＝（ｈ_ｉ，_ｊ）_Ｉ×Ｊとした場合、

となるような非負値行列である基底行列Ｗ，重み行列Ｈを求めるものが非負値行列分解（ＮＭＦ）である。 That is, a non-negative matrix X = (x _k , _j ) _{K × J} = (x ₁ , x ₂ ,... X _J ) in which _J K-dimensional non-negative vectors are arranged, I K-dimensional non-negative bases Vector-based basis matrix W = (w _k , _i ) _{K × I} = (w ₁ , w ₂ ,... W _I ), weight matrix H = (h _i , _j ) with non-negative coupling coefficients as elements _{If I × J} ,

Non-negative matrix decomposition (NMF) is a method for obtaining a base matrix W and a weight matrix H, which are non-negative matrixes such that

ここで、これらの近さの尺度を、フロベニウス（Ｆｒｏｂｅｎｉｕｓ）ノルム

と、Ｉダイバージェンス（またはＫＬダイバージェンス）

で定義した場合の最適な基底行列Ｗおよび重み行列Ｈを数値計算するための効率的な反復アルゴリズムが、開示されている（非特許文献１）。これらは乗法更新（ＭｕｌｔｉｐｌｉｃａｔｉｖｅＵｐｄａｔｅ）アルゴリズムと呼ばれている。なお、式６において、［・］_ｋ，_ｊは、「・」の行列の（ｋ，ｊ）成分を表す。 Here, the measure of these closenesses is expressed by the Frobenius norm.

And I divergence (or KL divergence)

An efficient iterative algorithm for numerically calculating the optimal base matrix W and weight matrix H when defined in (1) is disclosed (Non-Patent Document 1). These are called multiplicative update algorithms. In Equation 6, [•] _k , _j represents the (k, j) component of the matrix “•”.

Ｗ^ＴをＷの転置行列とすると、フロベニウスノルムで定義された式５に示すＦ（Ｘ，ＷＨ）を小さくする更新式は、それぞれ

で与えられ、Ｉダイバージェンスで定義された式６に示すＩ（Ｘ，ＷＨ）を小さくする更新式はそれぞれ

で与えられる。 When the W ^T and transposed matrix of W, update expression for reducing the F (X, WH) shown in Equation 5 defined in the Frobenius norm, respectively

The update equations for reducing I (X, WH) shown in Equation 6 defined by I divergence are

Given in.

いずれの更新式も、基底行列Ｗと重み行列Ｈの初期値が非負値行列であれば必ず非負値行列に更新されることが保証されている点に特徴があり、目的関数の定義によって異なる更新式が導入されるため、当然のことながら目的関数の定義によって異なる値に収束する。
したがって、基底行列Ｗおよび重み行列Ｈの行列同士の近さをフロベニウスノルムで測った場合で最適解を得ることは、目的関数を式５に示したフロベニウスノルムとして定義し、尤度関数を次式に示すガウス（Ｇａｕｓｓ）分布と仮定して最尤推定していることに相当する。

Each update formula is characterized in that it is guaranteed to be updated to a non-negative matrix whenever the initial values of the base matrix W and the weight matrix H are non-negative matrices, and the update varies depending on the definition of the objective function. Since formulas are introduced, it naturally converges to different values depending on the definition of the objective function.
Therefore, to obtain an optimal solution when the proximity of the base matrix W and the weight matrix H is measured by the Frobenius norm, the objective function is defined as the Frobenius norm shown in Equation 5, and the likelihood function is expressed by the following equation: This corresponds to the maximum likelihood estimation assuming the Gauss distribution shown in FIG.

また、基底行列Ｗおよび重み行列Ｈの行列同士の近さをＩダイバージェンスで測った場合で最適解を得ることは、目的関数を式６に示したＩダイバージェンスとして定義し、尤度関数を次式に示すポアソン（Ｐｏｉｓｓｏｎ）分布と仮定して最尤推定していることに相当する。

ただし、観測データｘ_ｋ，_ｊは量子化された整数値データ

とする。 Also, to obtain an optimal solution when the proximity of the base matrix W and the weight matrix H is measured by I divergence, the objective function is defined as I divergence shown in Equation 6, and the likelihood function is expressed by the following equation: This corresponds to the maximum likelihood estimation assuming the Poisson distribution shown in FIG.

However, observation data x _k , _j is quantized integer value data

And

すなわち、観測データが実際にどのような統計分布に従って生成されたかに従って、式５に示したフロベニウスノルムあるいは式６に示したＩダイバージェンスとしての目的関数のうち、最小にするべき目的関数がいずれであるかが決まるわけである。
しかしながら、一般に、観測データがどのような統計分布に従って生成されたかを事前に知ることはできず、かつ、必ずしもガウス分布やポアソン分布に理想的に従っているとは限らない。このため、従来、状況に応じてどちらかの目的関数を採用することが好ましいかを経験的に選択していた。
Ｄ．Ｄ．ＬｅｅａｎｄＨ．Ｓ．Ｓｅｕｎｇ，“ＬｅａｒｎｉｎｇｔｈｅＰａｒｔｓｏｆＯｂｊｅｃｔｓｂｙＮｏｎ−ｎｅｇａｔｉｖｅＭａｔｒｉｘＦａｃｔｏｒｉｚａｔｉｏｎ，”ＮａｔｕｒａＶｏｌ．４０１，ｐｐ．７８８−７９１，１９９９． That is, which of the objective functions to be minimized is the Frobenius norm shown in Equation 5 or the I divergence objective function shown in Equation 6 depending on what statistical distribution is actually generated according to the statistical distribution. It is decided.
However, in general, it is not possible to know in advance what kind of statistical distribution the observation data is generated in, and it does not always follow ideally the Gaussian distribution or the Poisson distribution. For this reason, conventionally, it has been empirically selected which one of the objective functions is preferably adopted depending on the situation.
D. D. Lee and H.C. S. Seung, “Learning the Part of Objects by Non-Negative Matrix Factorization,” Nature Vol. 401, pp. 788-791, 1999.

しかしながら、後で詳細に説明するが、図３に示す分布図からわかるように、ポアソン分布の尤度関数の場合、尤度関数の分布の裾は、観測値ｘが大きいほど広くなる。このため、観測値ｘ（観測行列の要素の値）が小さいほど、ＮＭＦモデル（基底行列Ｗと重み行列Ｈの積をとった行列）の要素値と観測値ｘとのずれに対して感度が鋭くなるようなコスト（尤度関数）を設けていることに相当する。一方、図４に示す分布図からわかるように、ガウス分布の尤度関数の場合、尤度関数の分布の裾は、観測値ｘの大小にかかわらず一定であるため、観測値ｘが大きいほど、ＮＭＦモデルの要素値と観測値ｘとのずれに対して感度が鋭くなるようなコストを設けていることに相当する。 However, as will be described later in detail, as can be seen from the distribution diagram shown in FIG. 3, in the case of the Poisson distribution likelihood function, the tail of the distribution of the likelihood function becomes wider as the observation value x is larger. For this reason, the smaller the observation value x (value of the element of the observation matrix), the more sensitive the deviation between the element value of the NMF model (matrix obtained by multiplying the base matrix W and the weight matrix H) and the observation value x is. This corresponds to providing a sharp cost (likelihood function). On the other hand, as can be seen from the distribution chart shown in FIG. 4, in the case of the likelihood function of the Gaussian distribution, the tail of the distribution of the likelihood function is constant regardless of the magnitude of the observed value x. This corresponds to providing a cost that makes the sensitivity sharp with respect to the deviation between the element value of the NMF model and the observed value x.

つまり、尤度関数をポアソン分布と仮定すると、小さな観測値ｘの微妙な誤差に影響を受けてＮＭＦモデルのパラメータの推定精度が低下しやすくなる。逆に、尤度関数をガウス分布と仮定すると、観測値ｘに含まれる外れ値や飛び値の影響を受けて、ＮＭＦモデルのパラメータの推定精度が低下しやすくなる。
また、音響信号のスペクトログラムのようなスパースなデータ（鋭いピーク成分がまばらにあり、それ以外は極端に小さい成分からなるデータ）に関しては、上記のいずれの影響も陽に受けることになる。したがって、尤度関数としてポアソン分布と仮定した場合は、小さい成分の影響を大きく受け、尤度関数としてガウス分布と仮定した場合は、鋭いピーク成分（外れ値や飛び値）の影響を大きく受けることになる。このため、どちらの尤度関数と仮定したとしても、ＮＭＦモデルのパラメータの推定精度が低下してしまうという問題がある。 That is, if the likelihood function is assumed to be Poisson distribution, the estimation accuracy of the parameters of the NMF model is likely to be deteriorated due to a subtle error of the small observation value x. On the other hand, assuming that the likelihood function is a Gaussian distribution, the estimation accuracy of the parameters of the NMF model is likely to decrease due to the influence of outliers and jump values included in the observed value x.
In addition, regarding sparse data such as a spectrogram of an acoustic signal (data having sharp peak components sparsely and otherwise consisting of extremely small components), any of the above effects is positively affected. Therefore, when assuming a Poisson distribution as a likelihood function, it is greatly affected by small components, and when assuming a Gaussian distribution as a likelihood function, it is greatly affected by sharp peak components (outliers and jump values). become. For this reason, whichever likelihood function is assumed, there is a problem that the estimation accuracy of the parameters of the NMF model is lowered.

本発明は、このような事情を考慮し、上記の問題を解決すべくなされたものであって、その目的は、ポアソン分布とガウス分布の混合分布により定義される関数を尤度関数とし、観測データに応じてＮＭＦパラメータだけでなく、観測データの統計分布の分布形も適応的に推定するための非負値行列分解の数値計算方法、非負値行列分解の数値計算装置、プログラムおよび記憶媒体を提供することにある。 The present invention has been made in consideration of such circumstances, and has been made to solve the above problem. The object of the present invention is to use a function defined by a mixture distribution of Poisson distribution and Gaussian distribution as a likelihood function, and to observe it. Provided are a non-negative matrix decomposition numerical calculation method, a non-negative matrix decomposition numerical calculation device, a program, and a storage medium for adaptively estimating not only the NMF parameters but also the distribution of the statistical distribution of observation data according to the data There is to do.

上記問題を解決するために、本発明は、非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算する非負値行列分解の数値計算方法において、初期パラメータ生成部が、前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成し、前記ガウス分布率行列生成部が、前記初期パラメータ生成部から入力された前記パラメータの初期値、あるいは、パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散の値に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成し、前記パラメータ更新部の分布混合比更新部が、前記ガウス分布率行列に基づき、前記分布混合比を更新し、前記パラメータ更新部のガウス分布分散更新部が、
前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新し、前記パラメータ更新部の基底行列更新部が、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新し、前記パラメータ更新部の重み行列更新部が、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新し、収束判定部は、前記パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かを判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成部に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記パラメータ更新部により更新されたパラメータの更新を再度判定し、分解行列出力部は、前記収束判定部により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、更新された前記基底行列および前記重み行列を出力することを特徴とする非負値行列分解の数値計算方法である。 In order to solve the above problem, the present invention is based on non-negative matrix decomposition that calculates parameters of a non-negative matrix decomposition model according to observation data based on observation data in which non-negative vectors storing non-negative data are arranged. In the numerical calculation method, the initial parameter generation unit is a parameter of the non-negative matrix decomposition model, each of a basis matrix, a weight matrix, a distribution mixture ratio representing a mixture ratio of a Gaussian distribution and a Poisson distribution, and a variance of a Gaussian distribution. The initial value of the parameter input from the initial parameter generation unit, or the base matrix updated by the parameter update unit, the weight matrix, and the distribution Based on the mixing ratio and the value of the variance of the Gaussian distribution, each element of the observation data is Gaussian and Poisson A Gaussian distribution rate matrix representing the probability generated from the Gaussian distribution of the cloth is generated, and the distribution mixture ratio update unit of the parameter update unit updates the distribution mixture ratio based on the Gaussian distribution rate matrix, The Gaussian distribution update part of the parameter update part
The variance of the Gaussian distribution is updated based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, and the observation data, and the basis matrix update unit of the parameter update unit includes the Gaussian distribution rate matrix, the basis matrix , Updating the basis matrix based on the variance of the weight matrix, the observation data, and the Gaussian distribution, and a weight matrix update unit of the parameter update unit includes the Gaussian distribution rate matrix, the basis matrix, the weight matrix, The weighting matrix is updated based on the observation data and the variance of the Gaussian distribution, and the convergence determination unit is configured to update the basis matrix, the weighting matrix, the distribution mixture ratio, and the Gaussian distribution updated by the parameter updating unit. It is determined whether each variance satisfies a predetermined criterion, and the updated base matrix, weight matrix, distribution mixture ratio, and distribution ratio are determined. And when the variance of the Gaussian distribution does not satisfy a predetermined criterion, the Gaussian distribution rate matrix generation unit outputs a signal instructing generation of the Gaussian distribution rate matrix, and the generated Gaussian distribution rate matrix Re-determining the update of the parameter updated by the parameter update unit based on the base matrix, the weight matrix, the distribution mixture ratio and the variance of the Gaussian distribution updated by the convergence determination unit Is a numerical calculation method for non-negative matrix decomposition, wherein the updated base matrix and weight matrix are output when it is determined that each satisfies a predetermined criterion.

また、本発明にかかる収束判定部は、前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、前記ポアソン分布とガウス分布の混合分布により定義された目的関数を最大化する所定の基準を満たしているか否かを判断することを特徴とする非負値行列分解の数値計算方法である。 In addition, the convergence determination unit according to the present invention maximizes an objective function in which the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution are defined by the Poisson distribution and the Gaussian distribution mixture. It is a numerical calculation method of non-negative matrix decomposition characterized by determining whether or not a predetermined criterion is satisfied.

また、本発明は、非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算する非負値行列分解の数値計算装置において、前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成する初期パラメータ生成部と、前記初期パラメータ生成部から入力された前記パラメータの初期値、あるいは、パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成するガウス分布率行列生成部と、前記ガウス分布率行列生成部から入力された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散を、それぞれ更新する基底行列更新部、重み行列更新部、分布混合比更新部、およびガウス分布分散更新部を有するパラメータ更新部と、前記パラメータ更新部により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かを判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成部に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記パラメータ更新部により更新されたパラメータの更新を再度判定する収束判定部と、前記収束判定部により、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、前記基底行列および前記重み行列を出力する分解行列出力部とを備え、前記パラメータ更新部は、前記ガウス分布率行列に基づき、前記分布混合比を更新する分布混合比更新部と、前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新するガウス分布分散更新部と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新する基底行列更新部と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新する重み行列更新部を有することを特徴とする非負値行列分解の数値計算装置である。 Further, the present invention provides a non-negative matrix decomposition numerical calculation device for calculating a parameter of a non-negative matrix decomposition model according to observation data based on observation data in which non-negative vectors storing non-negative data are arranged. An initial parameter generation unit for generating respective initial values of a basis matrix, a weight matrix, a distribution mixture ratio representing a mixture ratio of a Gaussian distribution and a Poisson distribution, and a variance of a Gaussian distribution, which are parameters of the non-negative matrix decomposition model; Based on the initial value of the parameter input from the initial parameter generation unit, or the basis matrix updated by the parameter update unit, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution, Gaussian component that represents the probability that each element is generated from Gaussian or Poisson distribution A Gaussian distribution rate matrix generation unit for generating a rate matrix, and a basis matrix for updating the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution input from the Gaussian distribution rate matrix generation unit, respectively A parameter updating unit having an updating unit, a weight matrix updating unit, a distribution mixture ratio updating unit, and a Gaussian distribution dispersion updating unit; the basis matrix updated by the parameter updating unit; the weight matrix; the distribution mixture ratio and the Gaussian It is determined whether the variance of the distribution satisfies a predetermined criterion, and the updated base matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution each satisfy a predetermined criterion. If not, the Gaussian distribution rate matrix generation unit outputs a signal instructing generation of the Gaussian distribution rate matrix, and the generated Gaussian distribution A convergence determination unit that re-determines update of the parameter updated by the parameter update unit based on the matrix, and the updated determination of the basis matrix, the weight matrix, the distribution mixture ratio, and the Gaussian distribution by the convergence determination unit. And a decomposition matrix output unit that outputs the basis matrix and the weight matrix when the variances are determined to satisfy predetermined criteria, respectively, and the parameter update unit is based on the Gaussian distribution rate matrix, A distribution mixture ratio update unit that updates a distribution mixture ratio; a Gaussian distribution dispersion update unit that updates a variance of the Gaussian distribution based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, and the observation data; and Update the basis matrix based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the variance of the Gaussian distribution And a weight matrix update unit that updates the weight matrix based on a variance of the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the Gaussian distribution. This is a numerical calculation device for non-negative matrix decomposition.

また、本発明にかかる収束判定部は、前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、前記ポアソン分布とガウス分布の混合分布により定義された目的関数を最大化する所定の基準を満たしているか否かを判断することを特徴とする非負値行列分解の数値計算装置である。 In addition, the convergence determination unit according to the present invention maximizes an objective function in which the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution are defined by the Poisson distribution and the Gaussian distribution mixture. It is a non-negative matrix decomposition numerical calculation apparatus characterized by determining whether or not a predetermined criterion is satisfied.

また、本発明は、非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算するコンピュータに、前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成する初期パラメータ生成手段と、前記初期パラメータ生成手段から入力された前記パラメータの初期値、あるいは、更新手段により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成するガウス分布率行列生成手段と、前記ガウス分布率行列に基づき、前記分布混合比を更新する更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新する更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新するガウス分布分散更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新する基底行列更新手段と、前記更新手段により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かの前記パラメータの更新を判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成手段に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記更新手段により更新された前記パラメータの更新を再度判定する収束判定手段と、前記収束判定手段により、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、前記基底行列および前記重み行列を出力する分解行列出力手段と、を実行させるためのプログラムである。 Further, the present invention provides a computer for calculating parameters of a non-negative matrix decomposition model corresponding to the observation data based on observation data in which non-negative vectors storing non-negative data are arranged, to the non-negative matrix decomposition model. Initial parameters generating means for generating initial values of parameters, ie, a base matrix, a weight matrix, a distribution mixture ratio representing a mixing ratio of Gaussian distribution and Poisson distribution, and a variance of Gaussian distribution, and input from the initial parameter generating means Each element of the observation data is Gaussian distribution and Poisson distribution on the basis of the initial value of the parameter, or the basis matrix updated by the updating means, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution A Gaussian component that generates a Gaussian distribution matrix representing the probability generated from the Gaussian distribution Based on the Gaussian distribution rate matrix, the updating unit for updating the distribution mixture ratio based on the Gaussian distribution rate matrix, the Gaussian distribution rate matrix, the basis matrix, the weight matrix, and the observation data, Updating means for updating the variance, Gaussian distribution variance updating means for updating the basis matrix based on the variance of the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the Gaussian distribution, and the Gaussian A basis matrix updating means for updating the weight matrix based on a distribution rate matrix, the basis matrix, the weight matrix, the observation data, and a variance of the Gaussian distribution; the basis matrix updated by the updating means; and the weight It is determined whether the parameter, whether the matrix, the distribution mixture ratio, and the variance of the Gaussian distribution satisfy predetermined criteria, is updated. If the updated base matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution do not satisfy predetermined criteria, the Gaussian distribution rate matrix generation means sends the Gaussian distribution rate matrix A signal for instructing generation is output, and a convergence determination unit for re-determining an update of the parameter updated by the update unit based on the generated Gaussian distribution rate matrix, and the updated by the convergence determination unit Decomposition matrix output means for outputting the basis matrix and the weight matrix when it is determined that the basis matrix, the weight matrix, the distribution mixture ratio and the variance of the Gaussian distribution satisfy predetermined criteria, respectively. Is a program for executing

また、本発明は、非負値データを格納した非負値ベクトルが並べられた観測データに基づき、この観測データに応じた非負値行列分解モデルのパラメータを計算するコンピュータに、前記非負値行列分解モデルのパラメータである、基底行列、重み行列、ガウス分布とポアソン分布の混合比率をあらわす分布混合比、およびガウス分布の分散のそれぞれの初期値を生成する初期パラメータ生成手段と、前記初期パラメータ生成手段から入力された前記パラメータの初期値、あるいは、更新手段により更新された前記基底行列、前記重み行列、前記分布混合比、および前記ガウス分布の分散に基づき、前記観測データの各要素がガウス分布およびポアソン分布のうちガウス分布から生成された確率を表したガウス分布率行列を生成するガウス分布率行列生成手段と、前記ガウス分布率行列に基づき、前記分布混合比を更新する更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、および前記観測データに基づき、前記ガウス分布の分散を更新する更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記基底行列を更新する更新手段と、前記ガウス分布率行列、前記基底行列、前記重み行列、前記観測データ、および前記ガウス分布の分散に基づき、前記重み行列を更新する更新手段と、前記更新手段により更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしているか否かの前記パラメータの更新を判定し、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていない場合、前記ガウス分布率行列生成手段に前記ガウス分布率行列の生成を指示する信号を出力し、この生成された前記ガウス分布率行列に基づき前記更新手段により更新された前記パラメータの更新を再度判定する収束判定手段と、前記収束判定手段により、前記更新された前記基底行列、前記重み行列、前記分布混合比および前記ガウス分布の分散が、それぞれ所定の基準を満たしていると判定された場合、前記基底行列および前記重み行列を出力する分解行列出力手段とを実行させるためのプログラムを記憶した記憶媒体である。 Further, the present invention provides a computer for calculating parameters of a non-negative matrix decomposition model corresponding to the observation data based on observation data in which non-negative vectors storing non-negative data are arranged, to the non-negative matrix decomposition model. Initial parameters generating means for generating initial values of parameters, ie, a base matrix, a weight matrix, a distribution mixture ratio representing a mixing ratio of Gaussian distribution and Poisson distribution, and a variance of Gaussian distribution, and input from the initial parameter generating means Each element of the observation data is Gaussian distribution and Poisson distribution on the basis of the initial value of the parameter, or the basis matrix updated by the updating means, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution A Gaussian component that generates a Gaussian distribution matrix representing the probability generated from the Gaussian distribution Based on the Gaussian distribution rate matrix, the updating unit for updating the distribution mixture ratio based on the Gaussian distribution rate matrix, the Gaussian distribution rate matrix, the basis matrix, the weight matrix, and the observation data, Updating means for updating the variance; updating means for updating the basis matrix based on the variance of the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the Gaussian distribution; and the Gaussian distribution rate matrix , The basis matrix, the weight matrix, the observation data, and update means for updating the weight matrix based on the variance of the Gaussian distribution, the basis matrix updated by the update means, the weight matrix, and the distribution mixture Determining whether the parameter and variance of the Gaussian distribution satisfy predetermined criteria, respectively, and updating the parameter; When the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution do not satisfy predetermined criteria, a signal that instructs the Gaussian distribution rate matrix generation unit to generate the Gaussian distribution rate matrix is output. And a convergence determination means for re-determining an update of the parameter updated by the update means based on the generated Gaussian distribution rate matrix, and the basis matrix and the weight matrix updated by the convergence determination means. Storing a program for executing the base matrix and the decomposition matrix output means for outputting the weight matrix when it is determined that the distribution mixture ratio and the variance of the Gaussian distribution satisfy predetermined criteria, respectively. Storage medium.

本発明によれば、観測データが、ポアソン分布ないしガウス分布のいずれかに従っている場合、どちらの分布に従っているかが推定されるため、そのもとで適切なＮＭＦパラメータを推定することができる。
また、本発明によれば、観測データがポアソン分布ないしガウス分布のいずれかに従っていない場合であっても、観測データの実際の統計分布にポアソン分布とガウス分布の混合分布で表現可能な範囲内で柔軟に尤度関数が適応されるため、そのもとで最尤のＮＭＦパラメータを得ることができる。
さらに、本発明によれば、効率的にＮＭＦパラメータを算出可能な乗法更新アルゴリズムを実現することができる。
また、本発明によれば、音響信号のスペクトログラムのようなスパースなデータに対しても、高精度なＮＭＦパラメータを推定することができる。 According to the present invention, when the observation data follows either a Poisson distribution or a Gaussian distribution, it is estimated which distribution is followed, and therefore an appropriate NMF parameter can be estimated based on the estimated distribution.
Further, according to the present invention, even if the observation data does not follow either the Poisson distribution or the Gaussian distribution, the actual statistical distribution of the observation data is within the range that can be expressed by the mixed distribution of the Poisson distribution and the Gaussian distribution. Since the likelihood function is flexibly adapted, the maximum likelihood NMF parameter can be obtained under that.
Furthermore, according to the present invention, it is possible to realize a multiplicative update algorithm that can efficiently calculate NMF parameters.
Further, according to the present invention, it is possible to estimate a highly accurate NMF parameter even for sparse data such as a spectrogram of an acoustic signal.

以下、図面を参照して、本発明の一実施形態について説明する。図１は、本発明の実施の形態の構成を示すブロック図であり、図２は、その動作を説明するフローチャートである。
図１に示すとおり、本実施の形態の非負値行列分解の計算装置は、初期パラメータ生成部１、ガウス分布率行列生成部２、パラメータ更新部３、収束判定部４および分解行列出力部５を有する。 Hereinafter, an embodiment of the present invention will be described with reference to the drawings. FIG. 1 is a block diagram showing a configuration of an embodiment of the present invention, and FIG. 2 is a flowchart for explaining the operation thereof.
As shown in FIG. 1, the non-negative matrix decomposition calculation apparatus according to the present embodiment includes an initial parameter generation unit 1, a Gaussian distribution rate matrix generation unit 2, a parameter update unit 3, a convergence determination unit 4, and a decomposition matrix output unit 5. Have.

初期パラメータ生成部１は、非負値行列分解（ＮＭＦ）モデルのパラメータである基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２の初期値を生成し、ガウス分布率行列生成部２に出力する。例えば、初期パラメータ生成部１は、観測データｘに基づき、所定の値が設定されたパラメータの初期値を記憶し、また、設定変更・修正が可能な構成を有する。
ただし、基底行列Ｗと重み行列Ｈは、それぞれＫ×Ｉ行列、Ｉ×Ｊ行列である。なお、ここで言うＫ，Ｊは、観測データ行列Ｘの行数および列数であり、Ｉは、基底ベクトルが備える非負値基底ベクトルの個数であって予め定められている任意の定数である。
それぞれの初期値の選択可能な範囲は、基底行列Ｗ＝（ｗ_ｋ，_ｉ）_Ｋ×Ｉに関しては、すべてのｋとｉに対しての範囲として、

また、重み行列Ｈ＝（ｈ_ｉ，_ｊ）_Ｉ×Ｊに関しては、すべてのｉとｊに対しての範囲として、

分布混合比λに関しては、

ガウス分布の分散σ^２に関しては、σ^２＞０とする。
なお、基底行列Ｗ＝（ｗ_ｋ，_ｉ）_Ｋ×Ｉは、Ｋ×Ｉ行列で（ｋ，ｉ）成分がｗ_ｋ，_ｉである行列（ｗ_１，ｗ_２，・・・ｗ_Ｉ）を表す。また、分布混合比λは、ガウス分布とポアソン分布の混合比率を表す。 The initial parameter generation unit 1 generates initial values of a base matrix W, a weight matrix H, a distribution mixture ratio λ, and a Gaussian distribution σ ² that are parameters of a non-negative matrix decomposition (NMF) model, and generates a Gaussian distribution rate matrix Output to part 2. For example, the initial parameter generation unit 1 stores an initial value of a parameter to which a predetermined value is set based on the observation data x, and has a configuration that allows setting change / correction.
However, the base matrix W and the weight matrix H are a K × I matrix and an I × J matrix, respectively. Here, K and J are the number of rows and the number of columns of the observation data matrix X, and I is the number of non-negative basis vectors included in the basis vector, and is a predetermined arbitrary constant.
The selectable range of each initial value is the range for all k and i with respect to the base matrix W = (w _k , _i ) _{K × I} ,

For the weight matrix H = (h _i , _j ) _{I × J} , the range for all i and j is

Regarding the distribution mixing ratio λ,

Regarding the variance σ ^{2 of the} Gaussian distribution, σ ² > 0.
Note that the base matrix W = (w _k , _i ) _{K × I} is a matrix (w ₁ , w ₂ ,... W _I ) in which the (k, i) component is w _k , _i in the K × I matrix. To express. The distribution mixture ratio λ represents the mixing ratio of the Gaussian distribution and the Poisson distribution.

ガウス分布率行列生成部２は、入力信号として入力された観測データ行列Ｘに基づき、初期パラメータ生成部１から入力された、あるいはパラメータ更新部３から読み出した基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２を用いて、次式に基づきｋ，ｊすべてについてガウス分布率γ_ｋ，_ｉを算出する。

なお、観測データ行列Ｘは、所定の観測データの出力部（図示せず）により、非負値データを格納した非負値ベクトルが並べられた観測データ行列である。なお、観測データｘ_ｋ，_ｊは、観測データ行列Ｘのｋ行ｊ行列目の成分、［・］_ｋ，_ｊは行列「・」のｋ行ｊ列目の成分を表す。
ガウス分布率行列生成部２は、上述の式１７を利用して計算されたガウス分布率γ_ｋ，_ｊを要素とした、次式に示すガウス分布率行列Γを生成し、パラメータ更新部３に出力する。

なお、ガウス分布率行列Γとは、観測データ行列Ｘの各要素ｘがガウス分布とポアソン分布のうち、ガウス分布から生成された確率を表したものである。
また、ガウス分布率行列生成部２は、入力された観測データ行列Ｘ、および初期パラメータ生成部１あるいは分解行列出力部５から入力された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２を、パラメータ更新部３に出力する。 The Gaussian distribution rate matrix generation unit 2 is based on the observation data matrix X input as an input signal, and is input from the initial parameter generation unit 1 or read out from the parameter update unit 3. Using the ratio λ and the variance σ ² of the Gaussian distribution, the Gaussian distribution ratios γ _k , _i are calculated for all k and j based on the following equation.

The observation data matrix X is an observation data matrix in which non-negative value vectors storing non-negative value data are arranged by a predetermined observation data output unit (not shown). Note that the observation data x _k , _j represents the component in the k-th row and j-th matrix of the observation data matrix X, and [•] _k , _j represents the component in the k-th row and j-th column of the matrix “•”.
The Gaussian distribution rate matrix generation unit 2 generates a Gaussian distribution rate matrix Γ shown in the following equation using the Gaussian distribution rates γ _k , _j calculated using the above Equation 17 as elements, and the parameter updating unit 3 Output.

The Gaussian distribution rate matrix Γ represents the probability that each element x of the observation data matrix X is generated from the Gaussian distribution among the Gaussian distribution and the Poisson distribution.
In addition, the Gaussian distribution rate matrix generation unit 2 includes the input observation data matrix X, the basis matrix W, the weight matrix H, the distribution mixture ratio λ, and the Gaussian distribution input from the initial parameter generation unit 1 or the decomposition matrix output unit 5. the variance sigma ^2, and outputs to the parameter updating unit 3.

パラメータ更新部３は、分布混合比更新部３１、ガウス分布分散更新部３２、基底行列更新部３３、重み行列更新部３４を含み、これらにより更新された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２を収束判定部４に出力する。 The parameter update unit 3 includes a distribution mixture ratio update unit 31, a Gaussian distribution dispersion update unit 32, a basis matrix update unit 33, and a weight matrix update unit 34, and the basis matrix W, weight matrix H, distribution mixture ratio updated thereby. The λ and the variance σ ^{2 of the} Gaussian distribution are output to the convergence determination unit 4.

分布混合比更新部３１は、ガウス分布率行列生成部２から入力されたガウス分布率行列Γの要素であるガウス分布率γ_ｋ，_ｊをガウス分布率γ´_ｋ，_ｊとし、これを用いて次式に基づき分布混合比λを更新する。

分布混合比更新部３１は、更新した分布混合比λを収束判定部４に出力する。 Distributive mixing ratio updating unit 31, the element a is the Gaussian distribution ratio gamma _k, _j of the Gaussian distribution ratio matrix Γ inputted from a Gaussian distribution index matrix generator 2 Gaussian distribution ratio gamma prime _k, and _j, by using the The distribution mixture ratio λ is updated based on the following equation.

The distribution mixture ratio update unit 31 outputs the updated distribution mixture ratio λ to the convergence determination unit 4.

ガウス分布分散更新部３２は、ガウス分布率行列生成部２から入力されたガウス分布の分散σ^２を次式に基づき更新する。つまり、ガウス分布分散更新部３２は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布率γ_ｋ，_ｊを、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布率γ´_ｋ，_ｊとし、これらと観測データ行列Ｘ＝（ｘ_ｋ，_ｊ）_Ｋ×Ｊを用いて、次式に基づきガウス分布の分散σ^２を更新する。

ガウス分布分散更新部３２は、更新したガウス分布の分散σ^２を収束判定部４に出力する。 The Gaussian distribution update unit 32 updates the variance σ ² of the Gaussian distribution input from the Gaussian distribution rate matrix generation unit 2 based on the following equation. That is, the Gaussian distribution variance updating unit 32 converts the base matrix W, the weight matrix H, and the Gaussian distribution rates γ _k and _j input from the Gaussian distribution rate matrix generation unit 2 into the base matrix W ′, the weight matrix H ′, and the Gauss, respectively. Using the distribution rate γ ′ _k , _j and these and the observed data matrix X = (x _k , _j ) _{K × J} , the variance σ ^{2 of the} Gaussian distribution is updated based on the following equation.

The Gaussian distribution variance updating unit 32 outputs the updated variance σ ² of the Gaussian distribution to the convergence determining unit 4.

基底行列更新部３３は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布行列Γ、およびガウス分布の分散σ^２を、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布行列Γ´、ガウス分布の分散σ´^２とし、これらと観測データ行列Ｘ＝（ｘ_ｋ，_ｊ）_Ｋ×Ｊを用いて、次式に基づき基底行列Ｗを更新する。

ただし、行列Ａ，Ｂ，Ｃ（太字）は、それぞれ次式と定義される。

行列の除法は、行列要素同士で商をとる演算とし、○中に中黒「・」で示す記号

は、Ｈａｄａｍａｒｄ積（行列の要素同士の積をとる演算）を、√内に「・」で示す表記は、「・」が行列のとき、行列の要素ごとの平方根をとる演算を表す。また、「１（太字）」は、要素がすべて１のＫ×Ｊの行列を表す。
基底行列更新部３３は、更新した基底行列Ｗを収束判定部４に出力する。 The base matrix update unit 33 converts the base matrix W, the weight matrix H, the Gaussian distribution matrix Γ, and the variance σ ^{2 of the} Gaussian distribution input from the Gaussian distribution rate matrix generation unit 2 into a base matrix W ′ and a weight matrix H ′, respectively. , Gaussian matrix gamma prime, as a dispersion Shiguma' ² of a Gaussian distribution, using them as observation data matrix _{_{X = (x k, j)}} K × J, updates the basis matrix W based on the following equation.

However, the matrices A, B, and C (bold) are respectively defined as the following equations.

Matrix division is an operation that takes a quotient between matrix elements, and a symbol indicated by a middle black “•” in ○

Indicates a Hadamard product (an operation that takes the product of elements of a matrix), and the notation indicated by “·” in √ represents an operation that takes the square root of each element of the matrix when “·” is a matrix. “1 (bold)” represents a K × J matrix having all elements of 1.
The base matrix update unit 33 outputs the updated base matrix W to the convergence determination unit 4.

重み行列更新部３４は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布行列Γ、およびガウス分布の分散σ^２を、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布行列Γ´、ガウス分布の分散σ´^２とし、これらと観測データ行列Ｘ＝（ｘ_ｋ，_ｊ）_Ｋ×Ｊを用いて、次式に基づき重み行列Ｈを更新する。

ただし、行列Ｄ，Ｅ，Ｆ（太字）は、それぞれ次式で定義される。

重み行列更新部３４は、更新した重み行列Ｈを収束判定部４に出力する。 The weight matrix updating unit 34 converts the basis matrix W, the weight matrix H, the Gaussian distribution matrix Γ, and the variance σ ^{2 of the} Gaussian distribution input from the Gaussian distribution rate matrix generation unit 2 into the basis matrix W ′ and the weight matrix H ′, respectively. , Gaussian matrix gamma prime, as a dispersion Shiguma' ² of a Gaussian distribution, using them as observation data matrix _{_{X = (x k, j)}} K × J, updates the weight matrix H based on the following equation.

However, the matrices D, E, and F (bold) are defined by the following equations, respectively.

The weight matrix update unit 34 outputs the updated weight matrix H to the convergence determination unit 4.

収束判定部４は、パラメータ更新部３により入力された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２が、それぞれ所定の基準を満たしているか否かのパラメータの更新を判定する。収束判定部４は、基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２が、それぞれ所定の基準を満たしていない場合、ガウス分布率行列生成部２にガウス分布率行列の生成を指示する信号を出力し、それぞれのパラメータが所定の基準を満たすまで、この生成されたガウス分布率行列Γに基づきパラメータ更新部３により更新されたパラメータの更新の判定を繰り返す。
本実施の形態において、収束判定部４は、パラメータ更新部３により、基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２が、下に示す式２６のポアソン分布とガウス分布の混合分布により定義された目的関数を最大化するよう更新され、この目的関数の値の変化率が所定値以下となったか否かを判定する。

収束判定部４は、式２６に示す目的関数の変化率が所定値以下となった場合、基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２が、所定の基準を満たしていると判定し、更新された最新の基底行列Ｗおよび重み行列Ｈを出力する指示の信号を、分解行列出力部５に出力する。 The convergence determination unit 4 updates the parameters as to whether the basis matrix W, the weight matrix H, the distribution mixture ratio λ, and the variance σ ^{2 of the} Gaussian distribution input by the parameter update unit 3 satisfy predetermined criteria, respectively. judge. When the basis matrix W, the weight matrix H, the distribution mixture ratio λ, and the variance σ ² of the Gaussian distribution do not satisfy the predetermined criteria, the convergence determination unit 4 causes the Gaussian distribution rate matrix generation unit 2 to A signal instructing generation is output, and the parameter update determination updated by the parameter update unit 3 is repeated based on the generated Gaussian distribution rate matrix Γ until each parameter satisfies a predetermined criterion.
In the present embodiment, the convergence determination unit 4 uses the parameter update unit 3 to determine that the basis matrix W, the weight matrix H, the distribution mixture ratio λ, and the variance σ ^{2 of the} Gaussian distribution are the Poisson distribution and the Gaussian distribution of Expression 26 shown below. It is updated so as to maximize the objective function defined by the mixed distribution, and it is determined whether or not the rate of change of the value of the objective function has become a predetermined value or less.

When the rate of change of the objective function shown in Expression 26 is equal to or less than a predetermined value, the convergence determination unit 4 satisfies the predetermined criterion when the base matrix W, the weight matrix H, the distribution mixture ratio λ, and the variance σ ^{2 of the} Gaussian distribution satisfy the predetermined criterion. And outputs an instruction signal for outputting the latest updated base matrix W and weight matrix H to the decomposition matrix output unit 5.

なお、収束判定部４は、上述の判定方法に限られず、例えば、パラメータの更新計算を反復して行った回数が、所定の回数に達したか否かを検出し、所定に回数に達した場合、分解行列出力部５に指示の信号を出力する構成であってもよい。また、収束判定部４は、パラメータの反復更新計算において、それぞれのパラメータの更新の変化率が所定の値以下であるか否かを検出し、変化率が所定値以下になった場合、分解行列出力部５に指示の信号を出力する構成であってもよい。 The convergence determination unit 4 is not limited to the above-described determination method. For example, the convergence determination unit 4 detects whether or not the number of repeated parameter update calculations has reached a predetermined number, and has reached the predetermined number. In this case, an instruction signal may be output to the decomposition matrix output unit 5. Further, the convergence determining unit 4 detects whether or not the change rate of each parameter update is equal to or less than a predetermined value in the parameter iterative update calculation. The configuration may be such that an instruction signal is output to the output unit 5.

分解行列出力部５は、収束判定部４から出力の指示信号を受けて、更新された最新の基底行列Ｗおよび重み行列Ｈを出力する。 The decomposition matrix output unit 5 receives the output instruction signal from the convergence determination unit 4 and outputs the updated latest base matrix W and weight matrix H.

この構成により、本発明は、ポアソン分布とガウス分布の混合分布により定義される尤度関数を、観測データｘに応じて、ＮＭＦパラメータだけでなく、観測データｘの統計分布の分布形も適応的に推定することができる。
つまり、観測データｘがポアソン分布ないしガウス分布のいずれかに従っている場合、どちらの分布に従っているか収束判定部４により判定されるまで、パラメータ更新部３による更新が繰り返される。このため、観測データｘがどちらの分布に従っているかが推定され、同時にそのもとで適切なＮＭＦパラメータが推定される。
一方、観測データｘがポアソン分布ないしガウス分布のいずれかに従っていない場合であっても、観測データｘの実際の統計分布にポアソン分布とガウス分布の混合分布で表現可能な範囲内で柔軟に尤度関数が適応されるため、そのもとで最尤のＮＭＦパラメータを得ることができる。
したがって、本発明は、観測データｘが、いずれの統計分布の分布形に従っているか不明な場合であっても、観測データｘに応じた最尤のＮＭＦパラメータを得ることができる。 With this configuration, the present invention adapts the likelihood function defined by the mixed distribution of Poisson distribution and Gaussian distribution not only according to the observation data x, but also in the distribution form of the statistical distribution of the observation data x. Can be estimated.
That is, when the observation data x follows either a Poisson distribution or a Gaussian distribution, the update by the parameter update unit 3 is repeated until the convergence determination unit 4 determines which distribution is being followed. For this reason, it is estimated which distribution the observation data x follows, and at the same time, an appropriate NMF parameter is estimated.
On the other hand, even when the observed data x does not follow either the Poisson distribution or the Gaussian distribution, the actual statistical distribution of the observed data x is flexible within a range that can be expressed by a mixed distribution of Poisson distribution and Gaussian distribution. Since the function is adapted, the most likely NMF parameter can be obtained under it.
Therefore, the present invention can obtain the maximum likelihood NMF parameter corresponding to the observation data x even if it is not clear which of the statistical distributions the observation data x follows.

さらに、本発明によれば、従来経験的に観測データｘの統計分布の分布形を判断していたところ、自動的に最適な統計分布に応じた最尤のＮＭＦパラメータを計算することができるため、効率的にＮＭＦパラメータを計算可能な乗法更新アルゴリズムを実現することができる。
また、本発明によれば、音響信号のスペクトログラムのようなスパースなデータに対しても、高精度なＮＭＦパラメータを推定することができる。 Furthermore, according to the present invention, since the distribution form of the statistical distribution of the observation data x is conventionally determined empirically, the maximum likelihood NMF parameter corresponding to the optimal statistical distribution can be automatically calculated. A multiplicative update algorithm that can efficiently calculate the NMF parameters can be realized.
Further, according to the present invention, it is possible to estimate a highly accurate NMF parameter even for sparse data such as a spectrogram of an acoustic signal.

次に、図２を用いて、本発明にかかる非負値行列分解の数値計算方法について説明する。図２は、本実施の形態における非負値行列分解の数値計算方法の一例を示すフローチャートである。 Next, a numerical calculation method for non-negative matrix decomposition according to the present invention will be described with reference to FIG. FIG. 2 is a flowchart showing an example of a numerical calculation method for non-negative matrix decomposition in the present embodiment.

図２に示すとおり、まず、初期パラメータ生成部１により、基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２の初期値が生成される（Ｓ１）。生成された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２の初期値は、ガウス分布率行列生成部２に出力される。
ガウス分布率行列生成部２は、入力された観測データ行列Ｘに基づき、これら基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２の初期値を用いて、式１７によりガウス分布率γ_ｋ，_ｉを計算し、このガウス分布率γ_ｋ，_ｉよりガウス分布率行列Γを生成する（Ｓ２）。ガウス分布率行列生成部２は、生成したガウス分布率行列Γと、観測データ行列Ｘ、基底行列Ｗ、重み行列Ｈ、分布混合比λおよびガウス分布の分散σ^２の初期値を、パラメータ更新部３に出力する。 As shown in FIG. 2, first, the initial parameter generation unit 1 generates initial values of the base matrix W, the weight matrix H, the distribution mixture ratio λ, and the variance σ ² of the Gaussian distribution (S1). Generated basis matrix W, the weighting matrix H, distributive mixing ratio lambda, the initial value of the variance sigma ² of the Gaussian distribution is outputted to a Gaussian distribution ratio matrix generating unit 2.
Based on the input observation data matrix X, the Gaussian distribution rate matrix generation unit 2 uses the initial values of the basis matrix W, the weight matrix H, the distribution mixture ratio λ, and the variance σ ² of the Gaussian distribution to obtain a Gaussian according to Expression 17. The distribution rate γ _k , _i is calculated, and a Gaussian distribution rate matrix Γ is generated from the Gaussian distribution rate γ _k , _i (S2). The Gaussian distribution rate matrix generation unit 2 uses the generated Gaussian distribution rate matrix Γ, the observed data matrix X, the base matrix W, the weight matrix H, the distribution mixture ratio λ, and the initial values of the Gaussian distribution σ ² as parameter update units. 3 is output.

パラメータ更新部３は、入力された情報に基づき、ＮＭＦパラメータの更新を行う。
すなわち、分布混合比更新部３１は、ガウス分布率行列生成部２から入力されたガウス分布率行列Γの要素であるガウス分布率γ_ｋ，_ｊをガウス分布率γ´_ｋ，_ｊとし、これを用いて式１９に基づき分布混合比λを更新する（Ｓ３）。
ガウス分布分散更新部３２は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布率γ_ｋ，_ｊを、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布率γ´_ｋ，_ｊとし、これらと観測データ行列Ｘ＝（ｘ_ｋ，_ｊ）_Ｋ×Ｊを用いて、式２０に基づきガウス分布の分散σ^２を更新する（Ｓ４）。
基底行列更新部３３は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布行列Γ、およびガウス分布の分散σ^２を、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布行列Γ´、ガウス分布の分散σ´^２とし、これらと観測データ行列Ｘ＝（ｘ_ｋ，_ｊ）_Ｋ×Ｊを用いて、式２１に基づき基底行列Ｗを更新する（Ｓ５）。
重み行列更新部３４は、ガウス分布率行列生成部２から入力された基底行列Ｗ、重み行列Ｈ、ガウス分布行列Γ、およびガウス分布の分散σ^２を、それぞれ基底行列Ｗ´、重み行列Ｈ´、ガウス分布行列Γ´、ガウス分布の分散σ´^２とし、これらと観測データ行列Ｘ＝（ｘ_ｋ，_ｊ）_Ｋ×Ｊを用いて、式２４に基づき重み行列Ｈを更新する（Ｓ６）。 The parameter update unit 3 updates the NMF parameter based on the input information.
In other words, distributive mixing ratio updating unit 31, a Gaussian distribution ratio gamma _k, _j is an element of a Gaussian distribution ratio matrix Γ inputted from a Gaussian distribution index matrix generator 2 Gaussian distribution ratio gamma prime _k, and _j, it The distribution mixture ratio λ is updated based on the equation 19 (S3).
The Gaussian distribution variance updating unit 32 converts the base matrix W, the weight matrix H, and the Gaussian distribution rates γ _k and _j input from the Gaussian distribution rate matrix generation unit 2 into a base matrix W ′, a weight matrix H ′, and a Gaussian distribution rate, respectively. Using γ ′ _k , _j and these and the observed data matrix X = (x _k , _j ) _{K × J} , the variance σ ^{2 of the} Gaussian distribution is updated based on Expression 20 (S4).
The base matrix update unit 33 converts the base matrix W, the weight matrix H, the Gaussian distribution matrix Γ, and the variance σ ^{2 of the} Gaussian distribution input from the Gaussian distribution rate matrix generation unit 2 into a base matrix W ′ and a weight matrix H ′, respectively. , Gaussian matrix gamma prime, as a dispersion Shiguma' ² of a Gaussian distribution, using them as observation data matrix _{_{X = (x k, j)}} K × J, updates the basis matrix W based on the equation 21 (S5).
The weight matrix updating unit 34 converts the basis matrix W, the weight matrix H, the Gaussian distribution matrix Γ, and the variance σ ^{2 of the} Gaussian distribution input from the Gaussian distribution rate matrix generation unit 2 into the basis matrix W ′ and the weight matrix H ′, respectively. The Gaussian distribution matrix Γ ′ and the Gaussian distribution σ ′ ² are used, and the observed data matrix X = (x _k , _j ) _{K × J} is used to update the weight matrix H based on Expression 24 (S6).

収束判定部４は、入力された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２が、それぞれ所定の基準を満たしているか否かを判定する（Ｓ７）。
基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２が、それぞれ所定の基準を満たしている場合（Ｓ７−ＹＥＳ）、収束判定部４は、更新された最新の基底行列Ｗおよび重み行列Ｈを出力する指示の信号を、分解行列出力５に出力する。
分解行列出力部５は、収束判定部４から出力の指示信号を受けて、更新された最新の基底行列Ｗおよび重み行列Ｈを出力する（Ｓ８）。 The convergence determination unit 4 determines whether or not the input base matrix W, weight matrix H, distribution mixture ratio λ, and Gaussian distribution variance σ ² satisfy predetermined criteria, respectively (S7).
When the base matrix W, the weight matrix H, the distribution mixture ratio λ, and the variance σ ^{2 of the} Gaussian distribution satisfy predetermined criteria (S7-YES), the convergence determination unit 4 determines that the updated latest base matrix W is updated. And an instruction signal for outputting the weight matrix H is output to the decomposition matrix output 5.
The decomposition matrix output unit 5 receives the output instruction signal from the convergence determination unit 4, and outputs the updated latest base matrix W and weight matrix H (S8).

一方、基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２が、それぞれ所定の基準を満たしていない場合（Ｓ７−ＮＯ）、収束判定部４は、ガウス分布率行列生成部２にガウス分布率行列の生成を指示する信号を出力する。
ステップＳ２に戻って、指示信号を受けたガウス分布率行列生成部２は、パラメータ更新部３から更新された最新の基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２を読み出し、観測データｘに基づきガウス分布率行列Γを生成し、パラメータ更新部３に出力する。パラメータ更新部３は、ガウス分布率行列生成部２から入力されたガウス分布率行列をおよび観測データｘに基づき基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２を更新し、収束判定部４に出力する。これら観測データｘおよびパラメータ等が入力された収束判定部４は、パラメータの更新の判定を行う。
すなわち、収束判定部４は、それぞれのパラメータが所定の基準を満たすまで、この生成されたガウス分布率行列Γに基づきパラメータ更新部３により更新されたパラメータの更新の判定を繰り返す。 On the other hand, when the basis matrix W, the weight matrix H, the distribution mixture ratio λ, and the variance σ ^{2 of the} Gaussian distribution do not satisfy the predetermined criteria (S7-NO), the convergence determining unit 4 is a Gaussian distribution rate matrix generating unit. 2 outputs a signal instructing generation of a Gaussian distribution rate matrix.
Returning to step S 2, the Gaussian distribution rate matrix generation unit 2 that has received the instruction signal obtains the latest base matrix W, weight matrix H, distribution mixture ratio λ, and Gaussian distribution variance σ ² updated from the parameter update unit 3. Read and generate a Gaussian distribution rate matrix Γ based on the observation data x and output it to the parameter update unit 3. The parameter update unit 3 updates the base matrix W, the weight matrix H, the distribution mixture ratio λ, and the Gaussian distribution σ ² based on the Gaussian distribution rate matrix input from the Gaussian distribution rate matrix generation unit 2 and the observation data x. And output to the convergence determination unit 4. The convergence determination unit 4 to which the observation data x and the parameters are input determines whether to update the parameters.
That is, the convergence determination unit 4 repeats the determination of updating the parameters updated by the parameter update unit 3 based on the generated Gaussian distribution rate matrix Γ until each parameter satisfies a predetermined criterion.

なお、パラメータ更新部３によるこれら更新の順番は任意であり、例えば、基底行列Ｗおよび重み行列Ｈの更新を先に実行し、その後に分布混合比λ、ガウス分布の分散σ^２の更新を実行する構成であってもよく、全ての更新を同時に実行する構成であってもよい。
パラメータ更新部３は、更新された基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２を収束判定部４に出力する。 The order of these updates by the parameter update unit 3 is arbitrary. For example, the base matrix W and the weight matrix H are updated first, and then the distribution mixture ratio λ and the Gaussian distribution σ ² are updated. The structure which performs may be sufficient and the structure which performs all the updates simultaneously may be sufficient.
The parameter update unit 3 outputs the updated base matrix W, weight matrix H, distribution mixture ratio λ, and Gaussian distribution variance σ ² to the convergence determination unit 4.

［観測データに基づく統計分布設定の問題］
ここで、ガウス分布とポアソン分布の尤度関数について説明する。
図３は、ポアソン分布の尤度関数の概形を示す図であり、図４は、ガウス分布の尤度関数の概形を示す図であり、図５は、ポアソン分布とガウス分布の混合分布の概形を示す図である。図３〜５において、凸型の複数の波形ａ〜ｏは、それぞれ異なる観測値ｘの尤度関数を示している。 [Problem of setting statistical distribution based on observation data]
Here, the likelihood function of the Gaussian distribution and the Poisson distribution will be described.
3 is a diagram showing an outline of the likelihood function of the Poisson distribution, FIG. 4 is a diagram showing an outline of the likelihood function of the Gaussian distribution, and FIG. 5 is a mixed distribution of the Poisson distribution and the Gaussian distribution. FIG. 3 to 5, a plurality of convex waveforms a to o indicate likelihood functions of different observation values x.

図３に示すとおり、ポアソン分布の尤度関数である複数の波形ａ〜ｅは、波形ａの観測データｘが最も小さく、アルファベット順に観測データｘが大きなり、それぞれ観測データｘ＝１，１０，３０，６０，９０であるときのポアソン分布の尤度関数を示す。図３に示すポアソン分布の概形は、横軸に観測データｘ、縦軸にポアソン分布の尤度関数（式８８参照）を規定する。
図３からわかるように、尤度関数がポアソン分布の場合、観測データｘが大きいほど尤度関数の分布の裾が広くなる。このため、本実施の形態では、観測データｘの値が大きいほど、ＮＭＦモデル（基底行列Ｗと重み行列Ｈの積をとった行列）［ＷＨ］_ｋ，_ｊの要素値と観測データｘ_ｋ，_ｊとのずれに対して感度が鈍くなるようなコスト（尤度関数）を設けていることに相当する。 As shown in FIG. 3, the plurality of waveforms a to e which are Poisson distribution likelihood functions have the smallest observed data x of the waveform a and the larger observed data x in alphabetical order, and the observed data x = 1, 10, The likelihood function of Poisson distribution when it is 30, 60, 90 is shown. The Poisson distribution outline shown in FIG. 3 defines the observation data x on the horizontal axis and the likelihood function (see Equation 88) of the Poisson distribution on the vertical axis.
As can be seen from FIG. 3, when the likelihood function is a Poisson distribution, the larger the observation data x, the wider the tail of the likelihood function distribution. Therefore, in the present embodiment, as the value of the observation data x is larger, the element value of the NMF model (matrix obtained by multiplying the base matrix W and the weight matrix H) [WH] _k , _j and the observation data x _k , _This corresponds to providing a cost (likelihood function) that makes the sensitivity less sensitive to a deviation from _j .

一方、図４に示すとおり、ガウス分布の尤度関数である複数の波形ｆ〜ｊは、波形ｆの観測データｘが最も小さく、アルファベット順に観測データｘが大きくなり、それぞれ観測データｘ＝１，１０，３０，６０，９０であるときのガウス分布の尤度関数を示す。図４に示すガウス分布の概形は、横軸に観測値、縦軸にガウス分布の尤度関数（式８９参照）を規定する。
図４からわかるように、尤度関数がガウス分布の場合、ポアソン分布の場合とは逆に、尤度関数の分布の裾は観測データｘの大きさによらず一定である。このため、本実施の形態では、ポアソン分布のときに比べて観測値ｘ_ｋ，_ｊの値が大きいほどＮＭＦモデル［ＷＨ］_ｋ，_ｊの要素値と観測値ｘとのずれに対して感度が鋭くなるようなコストを設けていることに相当する。 On the other hand, as shown in FIG. 4, in the plurality of waveforms f to j which are likelihood functions of the Gaussian distribution, the observation data x of the waveform f is the smallest, the observation data x is increased in alphabetical order, and the observation data x = 1, The likelihood function of the Gaussian distribution when 10, 30, 60, 90 is shown. The outline of the Gaussian distribution shown in FIG. 4 defines the observed value on the horizontal axis and the likelihood function (see Equation 89) of the Gaussian distribution on the vertical axis.
As can be seen from FIG. 4, when the likelihood function is a Gaussian distribution, the tail of the distribution of the likelihood function is constant regardless of the size of the observation data x, contrary to the case of the Poisson distribution. For this reason, in the present embodiment, the sensitivity of the difference between the element value of the NMF model [WH] _k , _j and the observation value x increases as the values of the observation values x _k , _j increase compared to the Poisson distribution. This is equivalent to providing sharper costs.

したがって、尤度関数がポアソン分布の場合のほうが、小さい観測データｘ_ｋ，_ｊにおける微量な誤差によって、推定結果全体が影響をうけやすく、一方、尤度関数がガウス分布の場合のほうが、観測データｘの外れ値や飛び値がある場合の影響を受けやすい。
特に音響信号のスペクトログラムのようにスパースなデータ（鋭いピーク成分がまばらにあり、それ以外は極端に小さい成分からなるデータ）に関しては、上記のいずれの影響も陽に受けることになるため、いずれか一方の尤度関数を仮定した場合、必然的に仮定した尤度関数の影響を受けることとなってしまう。 Therefore, when the likelihood function is a Poisson distribution, the entire estimation result is more easily affected by a small error in the small observation data x _k , _j , while the observation data is more likely when the likelihood function is a Gaussian distribution. It is easily affected by outliers and jump values of x.
Especially for sparse data such as spectrograms of acoustic signals (data consisting of sparsely sharp peak components and otherwise extremely small components), any of the above effects will be positively affected. When one likelihood function is assumed, it is inevitably influenced by the assumed likelihood function.

つまり、ＮＭＦの性能を向上させるためには、観測データが従う統計分布により即した目的関数をもとに、ＮＭＦを行うことが望ましい。そのためには、観測データに応じた統計分布を推定しながら、ＮＭＦパラメータである基底行列Ｗ，重み行列Ｈを推定する必要がある。
本発明は、ポアソン分布とガウス分布の混合分布により尤度関数を定義し、観測データに応じてＮＭＦパラメータだけでなく、観測データの統計分布の分布形も適応的に推定することができる。よって本発明は、これら統計分布と尤度関数とのバランスを観測データから適応的に調節し、適切な感度のコストを設けることが可能であり、いかなる応用場面でも最適化基準の選択に迫られることなく、統一的に利用可能なＮＭＦの枠組みを確率する可能性を発揮することができる。 That is, in order to improve the performance of NMF, it is desirable to perform NMF based on an objective function that conforms to the statistical distribution followed by the observation data. For this purpose, it is necessary to estimate the base matrix W and the weight matrix H, which are NMF parameters, while estimating the statistical distribution according to the observation data.
In the present invention, a likelihood function is defined by a mixed distribution of Poisson distribution and Gaussian distribution, and not only the NMF parameter but also the distribution form of the statistical distribution of the observation data can be adaptively estimated according to the observation data. Therefore, the present invention can adaptively adjust the balance between the statistical distribution and the likelihood function from the observation data, and can set an appropriate sensitivity cost, and is required to select an optimization criterion in any application situation. The possibility of probing an NMF framework that can be used in a unified manner can be exhibited.

［ハイブリッド尤度関数の定義］
次に、図５を参照して、ポアソン分布とガウス分布の混合分布について説明する。
図５に示すポアソン分布とガウス分布の混合分布は、

と定義される。ただし、ガウス分布は、

ポアソン分布は、

とする。
この式２７の混合分布の尤度関数は、式５に示すフロベニウスノルムおよび式６に示すＩダイバージェンスで定義される目的関数の両特性を含むハイブリッド尤度関数である。
ただし、ガウス分布とポアソン分布の分布混合比λは、０≦λ≦１である。すなわち、分布混合比λが大きいほど、混合分布はガウス分布に近似し、分布混合比λが小さいほど、混合分布はポアソン分布に近似する。 [Definition of hybrid likelihood function]
Next, a mixed distribution of Poisson distribution and Gaussian distribution will be described with reference to FIG.
The mixed distribution of Poisson distribution and Gaussian distribution shown in FIG.

Is defined. However, the Gaussian distribution is

Poisson distribution is

And
The likelihood function of the mixed distribution of Expression 27 is a hybrid likelihood function including both characteristics of the Frobenius norm shown in Expression 5 and the objective function defined by I divergence shown in Expression 6.
However, the distribution mixture ratio λ of the Gaussian distribution and the Poisson distribution is 0 ≦ λ ≦ 1. That is, the larger the distribution mixture ratio λ is, the closer the mixture distribution is to a Gaussian distribution, and the smaller the distribution mixture ratio λ is, the closer the mixture distribution is to a Poisson distribution.

したがって、対数尤度関数Ｌは、

となり、この対数尤度関数Ｌが最大となるような基底行列Ｗ、重み行列Ｈ、分布混合比λ、ガウス分布の分散σ^２を算出することが好ましい。なお、この分布混合比λとガウス分布の分散σ^２を自由パラメータとして扱うことで、観測データｘが従う確率分布をできるだけよくあらわすことができるようになる。つまり、統計分布と尤度関数とのバランスを観測データから適応的に調節し、適切な感度のコストを設けることが可能となる。
なお、本実施の形態において、パラメータ更新部３は、上述の対数尤度関数Ｌを目的関数として、この目的関数を最大化するように各パラメータを更新し、収束判定部４は、この目的関数の値の変化率が所定値以下となったか否かを判定している。 Therefore, the log likelihood function L is

Therefore, it is preferable to calculate the base matrix W, the weight matrix H, the distribution mixture ratio λ, and the variance σ ^{2 of the} Gaussian distribution that maximize the log likelihood function L. By treating the distribution mixture ratio λ and the Gaussian distribution variance σ ² as free parameters, the probability distribution followed by the observation data x can be represented as well as possible. That is, it is possible to adaptively adjust the balance between the statistical distribution and the likelihood function from the observation data, and to provide an appropriate sensitivity cost.
In the present embodiment, the parameter update unit 3 updates each parameter so as to maximize the objective function using the log likelihood function L described above as an objective function, and the convergence determination unit 4 uses the objective function. It is determined whether or not the rate of change of the value has become equal to or less than a predetermined value.

［ＮＭＦパラメータと分布パラメータの最尤推定］
次に、本実施の形態における乗法更新アルゴリズムを導くための、最適な補助関数の計算について説明する。
まず、対数関数の凸性に基づき、

という条件もとで、次のような不等式が成り立つ。

ただし、ガウス分布率行列Γは次式のとおりである。

また、

かつ

という条件のもとでは、次のような不等式が成り立つ。

さらに、

かつ

という条件のもとでは、次のような不等式が成り立つ。

[Maximum likelihood estimation of NMF parameters and distribution parameters]
Next, calculation of an optimal auxiliary function for deriving the multiplicative update algorithm in the present embodiment will be described.
First, based on the convexity of the logarithmic function,

Under the condition, the following inequality holds.

However, the Gaussian distribution rate matrix Γ is as follows.

Also,

And

Under the condition, the following inequality holds.

further,

And

Under the condition, the following inequality holds.

そこで、上述の不等式（式３２，３６，３９）において、いずれも

である場合に等号が成立するため、

が成り立つ。
だたし、３階のテンソルＭは

とする。 Therefore, in the above inequalities (Equations 32, 36, 39),

Since the equal sign holds when

Holds.
However, the tensor M on the 3rd floor is

And

以上より、次のような不等式が成り立つ。

従って、ガウス分布率γ_ｋ，_ｊが次のような場合、上述の不等式（式４３）が成り立つ。

また、ｍ_i，k，jが次のような場合、上述の不等式（式４４）が成り立つ。

From the above, the following inequality holds.

Therefore, when the Gaussian distribution rate γ _k , _j is as follows, the above inequality (formula 43) holds.

Further, when m _{i, k, j} is as follows, the above inequality (formula 44) holds.

このように、式４１，４４に基づき算出された関数Ｇ（Ｗ，Ｈ，λ，σ^２，Γ，Ｍ）は、対数尤度関数Ｌ（Ｗ，Ｈ，λ，σ^２）を効率的に局所最大化するための補助関数であり、式４５，４６による更新と、基底行列Ｗ、重み行列Ｈ、分布混合比λ、およびガウス分布の分散σ^２の更新を繰り返すことで、対数尤度関数Ｌ（Ｗ，Ｈ，λ，σ^２）を単調増加させることができる。
本実施の形態においては、上述の理論に基づき、ガウス分布率行列生成部２によるガウス分布率行列Γの更新、およびパラメータ更新部３による基底行列Ｗ、重み行列Ｈ、分布混合比λ、およびガウス分布の分散σ^２の更新を繰り返すことで、本発明にかかる非負値行列分解の数値計算装置は、対数尤度関数Ｌを単調増加させている。 As described above, the function G (W, H, λ, σ ² , Γ, M) calculated based on the equations 41 and 44 efficiently converts the log likelihood function L (W, H, λ, σ ² ). A logarithmic likelihood function that is an auxiliary function for local maximization, and is repeated by updating the formulas 45 and 46 and updating the base matrix W, the weight matrix H, the distribution mixture ratio λ, and the variance σ ² of the Gaussian distribution. L (W, H, λ, σ ² ) can be monotonously increased.
In the present embodiment, based on the above theory, the Gaussian distribution rate matrix generation unit 2 updates the Gaussian distribution rate matrix Γ, and the parameter update unit 3 uses the base matrix W, the weight matrix H, the distribution mixture ratio λ, and the Gaussian. By repeatedly updating the distribution variance σ ² , the numerical calculation device for non-negative matrix decomposition according to the present invention monotonically increases the log likelihood function L.

［乗法更新式の導入］
式４１により算出された補助関数Ｇに基づき算出される各パラメータの更新式について説明する。
まず、ガウス分布率行列Γ＝Γ´という条件のもとで、補助関数Ｇ（Ｗ，Ｈ，λ，σ^２，Γ，Ｍ）が最大となる分布混合比λを求める。

を解くと、

となる。
なお、ガウス分布率γ_ｋ，_ｉの式４５より

であるから

となる。 [Introduction of multiplication update formula]
An update formula for each parameter calculated based on the auxiliary function G calculated by Formula 41 will be described.
First, the distribution mixture ratio λ that maximizes the auxiliary function G (W, H, λ, σ ² , Γ, M) is obtained under the condition that the Gaussian distribution rate matrix Γ = Γ ′.

And solving

It becomes.
From Gaussian distribution rate γ _k , _i (Equation 45)

Because

It becomes.

次に、次式に基づき、基底行列Ｗ＝Ｗ´、重み行列Ｈ＝Ｈ´、ガウス分布率行列Γ＝Γ´、３階のテンソルＭ＝Ｍ´という条件のもとで、補助関数Ｇ（Ｗ，Ｈ，λ，σ^２，Γ，Ｍ）が最大となるガウス分布の分散σ^２を求める。

この式５１を、補助関数Ｇ（Ｗ，Ｈ，λ，σ^２，Γ，Ｍ）に代入して、

を解くと、

となる。 Next, on the basis of the following equation, the auxiliary function G () under the condition that the base matrix W = W ′, the weight matrix H = H ′, the Gaussian distribution rate matrix Γ = Γ ′, and the third-order tensor M = M ′ The variance σ ^{2 of the} Gaussian distribution that maximizes W, H, λ, σ ² , Γ, M) is obtained.

Substituting this equation 51 into the auxiliary function G (W, H, λ, σ ² , Γ, M),

And solving

It becomes.

次に、次式に基づき、重み行列Ｈ＝Ｈ´、ガウス分布率行列Γ＝Γ´、３階のテンソルＭ＝Ｍ´、ガウス分布の分散σ^２＝σ^´２という条件のもとで、補助関数Ｇ（Ｗ，Ｈ，λ，σ^２，Γ，Ｍ）が最大となる基底行列Ｗを求める。

この式５４に

を代入すると、

となる。ただし、行列の除法は、要素同士で商をとる演算とし、○中に中黒「・」で示す記号（数２３参照）は、Ｈａｄａｍａｒｄ積（行列の要素同士の積をとる演算）を、√内に「・」で示す表記は、「・」が行列のとき、行列の要素ごとの平方根をとる演算とする。なお、式５６において、「１（太字）」は要素がすべて１のＫ×Ｊの行列である。
この式５６の両辺に、

を乗じると、次に示す２次方程式となる。

ここで、

とすると、式５８の方程式の解は、

となる。
ただし、

のとき、すなわち、式４３の√の中が行列のとき、行列の要素ごとの平方根をとる演算子を表す。
以上により、基底行列Ｗの更新式は、

となる。式６２の解は、非負値行列であるため、式６４の基底行列Ｗ´が非負値行列であれば基底行列Ｗは必ず非負値行列に更新される。 Next, based on the following equation, under the condition of a weight matrix H = H ′, a Gaussian distribution rate matrix Γ = Γ ′, a third-order tensor M = M ′, and a Gaussian distribution σ ² = ^σ′2 . A basis matrix W that maximizes the auxiliary function G (W, H, λ, σ ² , Γ, M) is obtained.

In this formula 54

Substituting

It becomes. However, the matrix division is an operation that takes a quotient between elements, and a symbol indicated by a middle black “•” in ○ (see Equation 23) represents a Hadamard product (an operation that takes the product of elements of a matrix), √ The notation indicated by “•” is an operation that takes the square root of each element of the matrix when “•” is a matrix. In Expression 56, “1 (bold)” is a K × J matrix in which all elements are 1.
On both sides of Equation 56,

Is multiplied by the following quadratic equation.

here,

Then, the solution of equation 58 is

It becomes.
However,

In other words, that is, when √ in Equation 43 is a matrix, it represents an operator that takes the square root of each element of the matrix.
From the above, the update formula of the basis matrix W is

It becomes. Since the solution of Expression 62 is a non-negative matrix, if the base matrix W ′ of Expression 64 is a non-negative matrix, the base matrix W is always updated to a non-negative matrix.

次に、次式に基づき、基底行列Ｗ＝Ｗ´、ガウス分布率行列Γ＝Γ´、３階のテンソルＭ＝Ｍ´、ガウス分布の分散σ^２＝σ^´２という条件のもとで、補助関数Ｇ（Ｗ，Ｈ，λ，σ^２，Γ，Ｍ）が最大となる重み行列Ｈを求める。

この式６５に、式５５を代入すると、

となる。
この式６６の両辺に、

を乗じると、次に示す２次方程式となる。

Next, based on the following equation, under the condition that the base matrix W = W ′, the Gaussian distribution rate matrix Γ = Γ ′, the third-order tensor M = M ′, and the variance of the Gaussian distribution σ ² = σ ^′ ² , A weight matrix H that maximizes the auxiliary function G (W, H, λ, σ ² , Γ, M) is obtained.

Substituting equation 55 into equation 65,

It becomes.
On both sides of Equation 66,

Is multiplied by the following quadratic equation.

ここで、

とすると、式６８の方程式の解は、

となる。
以上により、重み行列Ｈの更新式は、

となる。基底行列Ｗと同様に、式７２の解は非負値行列であるため、式７３の重み行列Ｈ´が非負値行列であれば重み行列Ｈは必ず非負値行列に更新される。
本実施の形態において、パラメータ更新部３は、上述のとおり算出された式４８，５３，６４，７３に基づき、基底行列Ｗ、重み行列Ｈ、分布混合比λおよびガウス分布の分散σ^２を更新し、ガウス分布率行列生成部２により生成されたガウス分布率行列Γのもとで補助関数Ｇを最大化させている。つまり、この補助関数Ｇを最大化させるように更新を繰り返すことによって、収束判定部４は、目的関数Ｌを効率的に局所最大化させるための基底行列Ｗ、重み行列Ｈ、分布混合比λおよびガウス分布の分散σ^２を得て、パラメータの更新判定を実行することができる。 here,

Then, the solution of equation 68 is

It becomes.
From the above, the update formula of the weight matrix H is

It becomes. Similar to the base matrix W, since the solution of Expression 72 is a non-negative matrix, if the weight matrix H ′ of Expression 73 is a non-negative matrix, the weight matrix H is always updated to a non-negative matrix.
In the present embodiment, the parameter updating unit 3 updates the base matrix W, the weight matrix H, the distribution mixture ratio λ, and the Gaussian distribution σ ² based on the equations 48, 53, 64, 73 calculated as described above. The auxiliary function G is maximized under the Gaussian distribution rate matrix Γ generated by the Gaussian distribution rate matrix generation unit 2. That is, by repeating the update so as to maximize the auxiliary function G, the convergence determination unit 4 causes the base matrix W, the weight matrix H, the distribution mixture ratio λ, and the objective function L to be locally maximized efficiently. It is possible to obtain the variance σ ^{2 of the} Gaussian distribution and execute the parameter update determination.

［更新式の性質］
ここで、下に示す式７４および式７５のときの、

前記基底行列Ｗおよび重み行列Ｈの更新式（式６４，７３参照）の性質について説明する。
式７４のとき、ガウス分布率γ_ｋ，_ｊの式４５より

すなわち、ガウス分布率行列Γは、

である。
同様に、式７５のとき、ガウス分布率行列Γは、

である。 [Characteristics of update expression]
Here, in the case of Expression 74 and Expression 75 shown below,

The property of the update formula (see formulas 64 and 73) of the base matrix W and the weight matrix H will be described.
In the case of Expression 74, from Expression 45 of the Gaussian distribution rate γ _k , _j

That is, the Gaussian distribution rate matrix Γ is

It is.
Similarly, in Equation 75, the Gaussian distribution rate matrix Γ is

It is.

まず、ガウス分布率行列Γ´＝「１（太字）」すなわち、ガウス分布率行列Γ´がすべて１のＫ×Ｊの行列であるとき、

となり、これらを基底行列Ｗおよび重み行列Ｈの更新式（式６４，７３参照）に代入すると、基底行列Ｗおよび重み行列Ｈは、それぞれ

となり、この結果は、フロベニウスノルムで定義されたＦ（Ｘ，ＷＨ）を小さくする基底行列Ｗおよび重み行列Ｈの更新式（式７，８参照）とそれぞれ一致する。 First, when the Gaussian distribution rate matrix Γ ′ = “1 (bold)”, that is, when the Gaussian distribution rate matrix Γ ′ is a 1 × K × J matrix,

Substituting these into the update formulas of the base matrix W and the weight matrix H (see equations 64 and 73), the base matrix W and the weight matrix H are respectively

This result is in agreement with the update formulas of the base matrix W and the weight matrix H (see formulas 7 and 8) that reduce F (X, WH) defined by the Frobenius norm, respectively.

一方、

という条件のもとで、次式が成り立つ。

ただし、｜・｜は、「・」が行列のとき、行列の要素ごとの絶対値をとる演算子を表す。 on the other hand,

Under the condition:

However, | · | represents an operator that takes an absolute value for each element of the matrix when “·” is a matrix.

同様にして

では、

となり、この結果は、Ｉダイバージェンスで定義されたＩ（Ｘ，ＷＨ）を小さくする基底行列Ｗおよび重み行列Ｈの更新式（式９，１０参照）とそれぞれ一致する。 In the same way

Then

This result agrees with the update formulas of the base matrix W and the weight matrix H (see formulas 9 and 10) that reduce I (X, WH) defined by I divergence, respectively.

したがって、基底行列Ｗの更新式（式６４）および重み行列Ｈの更新式（式７３）は、２種類の乗法更新式を内包するハイブリッドな乗法更新式である。
本発明は、パラメータ更新部３により、このハイブリッドな乗法更新式を用いて基底行列Ｗおよび重み行列Ｈを更新している。
このため、本発明は、ポアソン分布とガウス分布の混合分布により尤度関数を定義し、観測データに応じてＮＭＦパラメータだけでなく、観測データの統計分布の分布形も適用的に推定することができる。 Therefore, the update formula (formula 64) of the base matrix W and the update formula (formula 73) of the weight matrix H are hybrid multiplicative update formulas including two types of multiplicative update formulas.
In the present invention, the parameter updating unit 3 updates the base matrix W and the weight matrix H using this hybrid multiplication update formula.
For this reason, the present invention defines a likelihood function by a mixed distribution of Poisson distribution and Gaussian distribution, and can appropriately estimate not only the NMF parameter but also the distribution form of the statistical distribution of the observation data according to the observation data. it can.

なお、上述の図３に示したポアソン分布の尤度関数の概形は、次式に示すポアソン分布の尤度関数の概形である。

また、上述の図４に示したガウス分布の尤度関数の概形は、次式に示すガウス分布の尤度関数の概形である。

Note that the approximate shape of the likelihood function of the Poisson distribution shown in FIG. 3 is the approximate shape of the likelihood function of the Poisson distribution shown in the following equation.

Further, the general form of the likelihood function of the Gaussian distribution shown in FIG. 4 is an outline of the likelihood function of the Gaussian distribution shown in the following equation.

ここで、以下上述の数式についての記号の定義について説明する。
非負値行列とは、全成分が非負（０以上）の値で構成される行列である。
次式に示す記号は、「全てのｋ，ｉに対して」という意味を表している。

次式に示す記号は、「行列Ｈ´の転置行列」という意味を表している。

次式に示す記号は、「左辺の式を右辺の式で定義する」という意味を表している。

次式に示す記号は、「Ｋ次元非負値ベクトルの空間」という意味を表している。

次式に示す記号は、「ベクトルｗの全ての成分が非負（０以上）である」という意味を表している。

次式に示す記号は、「左辺の式を右変の式で近似する」という意味を表している。

式９６に示す記号は、式９７に示す「フロベニウスノルム」を表している。

次式に示す記号は、「右辺を左辺に代入する（左辺の式を右辺の式で置き換える）」という意味を表している。

次式に示す記号は、「整数全体からなる集合」という意味を表している。

次式に示す記号は、「０から１の間（区間）に含まれる（０＜ガウス分布率γ_ｋ，_ｊ＜１）」という意味を表している。

Here, the definition of symbols for the above-described mathematical expressions will be described below.
A non-negative matrix is a matrix in which all components are composed of non-negative (0 or more) values.
The symbol shown in the following expression represents the meaning “for all k and i”.

The symbol shown in the following expression represents the meaning of “transposed matrix of matrix H ′”.

The symbol shown in the following expression represents the meaning that “the expression on the left side is defined by the expression on the right side”.

The symbol shown in the following expression represents the meaning of “K-dimensional non-negative vector space”.

The symbol shown in the following expression represents the meaning that “all components of the vector w are non-negative (0 or more)”.

The symbol shown in the following expression represents the meaning of “approximate the expression on the left side with a right-changing expression”.

The symbol shown in Expression 96 represents “Frobenius norm” shown in Expression 97.

The symbol shown in the following expression represents the meaning of “substitute the right side for the left side (replace the left side expression with the right side expression)”.

The symbol shown in the following formula represents the meaning of “set consisting of whole integers”.

The symbol shown in the following expression represents the meaning of “0 to 1 (interval) included (0 <Gaussian distribution ratio γ _k , _j <1)”.

なお、上述の初期パラメータ生成部１、ガウス分布率行列生成部２、パラメータ更新部３、収束判定部４、分解行列出力部５における、図２を用いて説明したそれぞれの動作の過程は、コンピュータに実行させるためのプログラムや、このプログラムとしてコンピュータ読み取り可能な記録媒体として利用可能であり、コンピュータシステムが読み出して実行することによって、上記処理が行われる。なお、ここでいう「コンピュータシステム」とは、ＣＰＵ及び各種メモリやＯＳ、周辺機器等のハードウェアを含むものである。
また、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、フラッシュメモリ等の書き込み可能な不揮発性メモリ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。 In the above-described initial parameter generation unit 1, Gaussian distribution rate matrix generation unit 2, parameter update unit 3, convergence determination unit 4, and decomposition matrix output unit 5, the operation processes described with reference to FIG. This program can be used as a computer-readable recording medium or a computer-readable recording medium as the program, and the computer system reads and executes the above-described processing. The “computer system” here includes a CPU, various memories, an OS, and hardware such as peripheral devices.
Further, the “computer system” includes a homepage providing environment (or display environment) if a WWW system is used.
The “computer-readable recording medium” means a flexible disk, a magneto-optical disk, a ROM, a writable nonvolatile memory such as a flash memory, a portable medium such as a CD-ROM, a hard disk built in a computer system, etc. This is a storage device.

さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムが送信された場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリ（例えばＤＲＡＭ（Dynamic Random Access Memory））のように、一定時間プログラムを保持しているものも含むものとする。
また、上記プログラムは、このプログラムを記憶装置等に格納したコンピュータシステムから、伝送媒体を介して、あるいは、伝送媒体中の伝送波により他のコンピュータシステムに伝送されてもよい。ここで、プログラムを伝送する「伝送媒体」は、インターネット等のネットワーク（通信網）や電話回線等の通信回線（通信線）のように情報を伝送する機能を有する媒体のことをいう。
また、上記プログラムは、前述した機能の一部を実現するためのものであっても良い。さらに、前述した機能をコンピュータシステムに既に記録されているプログラムとの組合せで実現できるもの、いわゆる差分ファイル（差分プログラム）であっても良い。 Further, the “computer-readable recording medium” means a volatile memory (for example, DRAM (Dynamic DRAM) in a computer system that becomes a server or a client when a program is transmitted through a network such as the Internet or a communication line such as a telephone line. Random Access Memory)), etc., which hold programs for a certain period of time.
The program may be transmitted from a computer system storing the program in a storage device or the like to another computer system via a transmission medium or by a transmission wave in the transmission medium. Here, the “transmission medium” for transmitting the program refers to a medium having a function of transmitting information, such as a network (communication network) such as the Internet or a communication line (communication line) such as a telephone line.
The program may be for realizing a part of the functions described above. Furthermore, what can implement | achieve the function mentioned above in combination with the program already recorded on the computer system, and what is called a difference file (difference program) may be sufficient.

本実施形態による非負値行列分解の数値計算装置の一例を示す概略ブロック図である。It is a schematic block diagram which shows an example of the numerical calculation apparatus of the nonnegative matrix decomposition | disassembly by this embodiment. 本実施形態による非負値行列分解の数値計算方法の一例を示すフローチャートである。It is a flowchart which shows an example of the numerical calculation method of the nonnegative matrix decomposition | disassembly by this embodiment. ポアソン分布の尤度関数の概形である。It is a rough form of the likelihood function of Poisson distribution. ガウス分布の尤度関数の概形である。It is a rough form of the likelihood function of Gaussian distribution. ポアソン分布とガウス分布の混合分布の概形である。It is an outline of a mixed distribution of Poisson distribution and Gaussian distribution.

Explanation of symbols

１初期パラメータ生成部
２ガウス分布率行列生成部
３パラメータ更新部
４収束判定部
５分解行列出力部
３１分布混合比更新部
３２ガウス分布分散更新部
３３基底行列更新部
３４重み行列更新部 DESCRIPTION OF SYMBOLS 1 Initial parameter generation part 2 Gaussian distribution rate matrix generation part 3 Parameter update part 4 Convergence determination part 5 Decomposition matrix output part 31 Distribution mixture ratio update part 32 Gaussian distribution dispersion | distribution update part 33 Basis matrix update part 34 Weight matrix update part

Claims

In the non-negative matrix decomposition numerical calculation method for calculating the parameters of the non-negative matrix decomposition model according to the observation data based on the observation data in which non-negative vectors storing non-negative data are arranged,
The initial parameter generator
Generating respective initial values of the base matrix, weight matrix, distribution mixture ratio representing the mixture ratio of Gaussian distribution and Poisson distribution, and variance of Gaussian distribution, which are parameters of the non-negative matrix decomposition model;
The Gaussian distribution rate matrix generation unit,
Based on the initial values of the parameters input from the initial parameter generation unit or the values of the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution updated by the parameter update unit Generate a Gaussian distribution rate matrix that represents the probability that each element of the data was generated from a Gaussian or Poisson distribution,
The distribution mixture ratio update unit of the parameter update unit,
Updating the distribution mixture ratio based on the Gaussian distribution rate matrix;
The Gaussian distribution update unit of the parameter update unit,
Updating the variance of the Gaussian distribution based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, and the observation data;
The base matrix update unit of the parameter update unit,
Updating the basis matrix based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the variance of the Gaussian distribution;
The weight matrix update unit of the parameter update unit,
Updating the weight matrix based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the variance of the Gaussian distribution;
The convergence judgment unit
Determining whether the basis matrix, the weight matrix, the distribution mixture ratio and the variance of the Gaussian distribution updated by the parameter update unit each satisfy a predetermined criterion;
When the updated base matrix, weight matrix, distribution mixture ratio, and variance of the Gaussian distribution do not satisfy predetermined criteria, the Gaussian distribution rate matrix generation unit generates the Gaussian distribution rate matrix. Output a signal to instruct, determine again the update of the parameter updated by the parameter update unit based on the generated Gaussian distribution rate matrix,
The decomposition matrix output section
When it is determined that the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution updated by the convergence determination unit satisfy predetermined criteria, respectively, the updated basis matrix and the A numerical calculation method for non-negative matrix decomposition, wherein a weight matrix is output.

The convergence determination unit satisfies a predetermined criterion that the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution satisfy the objective function defined by the Poisson distribution and the Gaussian mixture distribution. 2. The numerical calculation method of non-negative matrix decomposition according to claim 1, wherein whether or not it is present is determined.

In a numerical calculation device for non-negative matrix decomposition that calculates parameters of a non-negative matrix decomposition model according to observation data based on observation data in which non-negative vectors storing non-negative data are arranged,
An initial parameter generation unit for generating respective initial values of a basis matrix, a weight matrix, a distribution mixture ratio representing a mixture ratio of a Gaussian distribution and a Poisson distribution, and a variance of a Gaussian distribution, which are parameters of the non-negative matrix decomposition model;
Based on the initial value of the parameter input from the initial parameter generation unit, or the basis matrix updated by the parameter update unit, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution, A Gaussian distribution rate matrix generation unit for generating a Gaussian distribution rate matrix representing the probability that each element is generated from a Gaussian distribution among Gaussian distribution and Poisson distribution;
Basis matrix update unit, weight matrix update unit, and distribution mixture ratio update unit that update the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution input from the Gaussian distribution rate matrix generation unit, respectively. And a parameter updating unit having a Gaussian distribution updating unit,
It is determined whether the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution updated by the parameter update unit each satisfy a predetermined criterion, and the updated basis matrix, When the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution do not satisfy predetermined criteria, a signal that instructs the Gaussian distribution rate matrix generation unit to generate the Gaussian distribution rate matrix is output. A convergence determination unit that again determines update of the parameter updated by the parameter update unit based on the generated Gaussian distribution rate matrix;
When the convergence determination unit determines that the updated base matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution satisfy predetermined criteria, respectively, the base matrix and the weight A decomposition matrix output unit for outputting a matrix,
The parameter update unit
A distribution mixture ratio update unit for updating the distribution mixture ratio based on the Gaussian distribution rate matrix;
A Gaussian distribution update unit for updating a variance of the Gaussian distribution based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, and the observation data;
A base matrix update unit that updates the base matrix based on the variance of the Gaussian distribution rate matrix, the base matrix, the weight matrix, the observation data, and the Gaussian distribution;
A numerical value of non-negative matrix decomposition, comprising: a weight matrix updating unit that updates the weight matrix based on the variance of the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the Gaussian distribution Computing device.

The convergence determination unit satisfies a predetermined criterion that the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution satisfy the objective function defined by the Poisson distribution and the Gaussian mixture distribution. 4. The numerical calculation device for non-negative matrix decomposition according to claim 3, wherein whether or not the non-negative value matrix is decomposed is determined.

Based on observation data in which non-negative vectors containing non-negative data are arranged, a computer that calculates the parameters of the non-negative matrix decomposition model corresponding to this observation data,
Initial parameter generating means for generating respective initial values of a basis matrix, a weight matrix, a distribution mixture ratio representing a mixture ratio of a Gaussian distribution and a Poisson distribution, and a variance of a Gaussian distribution, which are parameters of the non-negative matrix decomposition model;
Based on the initial value of the parameter input from the initial parameter generation means, or the basis matrix updated by the update means, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution, A Gaussian distribution rate matrix generating means for generating a Gaussian distribution rate matrix representing the probability that an element is generated from a Gaussian distribution among Gaussian distribution and Poisson distribution;
Updating means for updating the distribution mixture ratio based on the Gaussian distribution rate matrix;
Updating means for updating a variance of the Gaussian distribution based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, and the observation data;
Gaussian distribution variance updating means for updating the basis matrix based on the variance of the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the Gaussian distribution;
Basis matrix update means for updating the weight matrix based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the variance of the Gaussian distribution;
The update of the parameter is determined whether or not each of the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution updated by the updating unit satisfies a predetermined criterion, and the updated When the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution do not satisfy predetermined criteria, a signal for instructing the Gaussian distribution rate matrix generation unit to generate the Gaussian distribution rate matrix is provided. A convergence determination unit that outputs and determines again the update of the parameter updated by the update unit based on the generated Gaussian distribution rate matrix;
When the convergence determining means determines that the updated base matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution satisfy predetermined criteria, respectively, the base matrix and the weight Decomposition matrix output means for outputting a matrix;
A program for running

Based on observation data in which non-negative vectors containing non-negative data are arranged, a computer that calculates the parameters of the non-negative matrix decomposition model corresponding to this observation data,
Initial parameter generating means for generating respective initial values of a basis matrix, a weight matrix, a distribution mixture ratio representing a mixture ratio of a Gaussian distribution and a Poisson distribution, and a variance of a Gaussian distribution, which are parameters of the non-negative matrix decomposition model;
Based on the initial value of the parameter input from the initial parameter generation means, or the basis matrix updated by the update means, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution, A Gaussian distribution rate matrix generating means for generating a Gaussian distribution rate matrix representing the probability that an element is generated from a Gaussian distribution among Gaussian distribution and Poisson distribution;
Updating means for updating the distribution mixture ratio based on the Gaussian distribution rate matrix;
Updating means for updating a variance of the Gaussian distribution based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, and the observation data;
Updating means for updating the basis matrix based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the variance of the Gaussian distribution;
Updating means for updating the weight matrix based on the Gaussian distribution rate matrix, the basis matrix, the weight matrix, the observation data, and the variance of the Gaussian distribution;
The update of the parameter is determined whether or not each of the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution updated by the updating unit satisfies a predetermined criterion, and the updated When the basis matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution do not satisfy predetermined criteria, a signal for instructing the Gaussian distribution rate matrix generation unit to generate the Gaussian distribution rate matrix is provided. A convergence determination unit that outputs and determines again the update of the parameter updated by the update unit based on the generated Gaussian distribution rate matrix;
When the convergence determining means determines that the updated base matrix, the weight matrix, the distribution mixture ratio, and the variance of the Gaussian distribution satisfy predetermined criteria, respectively, the base matrix and the weight A storage medium storing a program for executing decomposition matrix output means for outputting a matrix.