(57)【要約】【課題】高速で安定的に収束する非負値行列分解技術を提供する。【解決手段】Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化部を含む非負値行列分解最適化装置であって、最適化部が最適化に用いるコスト関数Jは、式J(A, B)=G(Y||ABT)+HA(A)+HB(B)（ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）で定義されるものであり、最適化部は、ブレグマン単調作用素分解に基づいて行列{A, B}を最適化する。【選択図】図７
2020126441.tif 000002
【特許請求の範囲】
【請求項１】
Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、
観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化部を含む非負値行列分解最適化装置であって、
前記最適化部が最適化に用いるコスト関数Jは、式
【数４９】
2020126441.tif 000051
（ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）
で定義されるものであり、
前記最適化部は、
次式により、行列Zを更新する行列Z更新部と、
【数５０】
2020126441.tif 000052
（ただし、HZは行列Zに対する正規化項、BD^+は関数D+を用いて定義されるブレグマンダイバージェンス）
次式により、双対変数~Vを更新する第１双対変数更新部と、
【数５１】
2020126441.tif 000053
次式により、行列Aを更新する行列A更新部と、
【数５２】
2020126441.tif 000054
次式により、行列Bを更新する行列B更新部と、
【数５３】
2020126441.tif 000055
次式により、双対変数~Xを更新する第２双対変数更新部と、
【数５４】
2020126441.tif 000056
を含む非負値行列分解最適化装置。
【請求項２】
Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、
観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化部を含む非負値行列分解最適化装置であって、
前記最適化部が最適化に用いるコスト関数Jは、式
【数５５】
2020126441.tif 000057
（ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）
で定義されるものであり、
前記最適化部は、
次式により、行列Zを更新する行列Z更新部と、
【数５６】
2020126441.tif 000058
（ただし、HZは行列Zに対する正規化項、BD^+は関数D+を用いて定義されるブレグマンダイバージェンス）
次式により、双対変数~Vを更新する第１双対変数更新部と、
【数５７】
2020126441.tif 000059
次式により、行列Aを更新する行列A更新部と、
【数５８】
2020126441.tif 000060
次式により、行列Bを更新する行列B更新部と、
【数５９】
2020126441.tif 000061
次式により、双対変数~Xを更新する第２双対変数更新部と、
【数６０】
2020126441.tif 000062
（ただし、ξは0<ξ<1を満たす定数）
を含む非負値行列分解最適化装置。
【請求項３】
請求項１または２に記載の非負値行列分解最適化装置であって、
前記行列Z更新部、前記行列A更新部、前記行列B更新部で用いるブレグマンダイバージェンスBD^+の定義に用いる関数D+は、式
【数６１】
2020126441.tif 000063
（ただし、Cはcijを(i, j)要素とする行列、hij=∇ij2G(cijold)（cijoldは現在の行列Cの(i, j) 要素、ε(>0)は所定の定数）
により与えられる
ことを特徴とする非負値行列分解最適化装置。
【請求項４】
請求項３に記載の非負値行列分解最適化装置であって、
関数D+は、hijを(i, j)要素とする行列Hに対してH=hAhBTを満たすベクトルhA, hBを用いて計算される
ことを特徴とする非負値行列分解最適化装置。
【請求項５】
Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、
非負値行列分解最適化装置が、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化ステップを実行する非負値行列分解最適化方法であって、
前記最適化ステップで最適化に用いるコスト関数Jは、式
【数６２】
2020126441.tif 000064
（ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）
で定義されるものであり、
前記最適化ステップは、
次式により、行列Zを更新する行列Z更新ステップと、
【数６３】
2020126441.tif 000065
（ただし、HZは行列Zに対する正規化項、BD^+は関数D+を用いて定義されるブレグマンダイバージェンス）
次式により、双対変数~Vを更新する第１双対変数更新ステップと、
【数６４】
2020126441.tif 000066
次式により、行列Aを更新する行列A更新ステップと、
【数６５】
2020126441.tif 000067
次式により、行列Bを更新する行列B更新ステップと、
【数６６】
2020126441.tif 000068
次式により、双対変数~Xを更新する第２双対変数更新ステップと、
【数６７】
2020126441.tif 000069
を含む非負値行列分解最適化方法。
【請求項６】
Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、
非負値行列分解最適化装置が、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化ステップを実行する非負値行列分解最適化方法であって、
前記最適化ステップで最適化に用いるコスト関数Jは、式
【数６８】
2020126441.tif 000070
（ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）
で定義されるものであり、
前記最適化ステップは、
次式により、行列Zを更新する行列Z更新ステップと、
【数６９】
2020126441.tif 000071
（ただし、HZは行列Zに対する正規化項、BD^+は関数D+を用いて定義されるブレグマンダイバージェンス）
次式により、双対変数~Vを更新する第１双対変数更新ステップと、
【数７０】
2020126441.tif 000072
次式により、行列Aを更新する行列A更新ステップと、
【数７１】
2020126441.tif 000073
次式により、行列Bを更新する行列B更新ステップと、
【数７２】
2020126441.tif 000074
次式により、双対変数~Xを更新する第２双対変数更新ステップと、
【数７３】
2020126441.tif 000075
（ただし、ξは0<ξ<1を満たす定数）
を含む非負値行列分解最適化方法。
【請求項７】
請求項１ないし４のいずれか１項に記載の非負値行列分解最適化装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、非負値行列を２つの非負値行列の積として表す非負値行列分解技術に関する。
【背景技術】
【０００２】
観測信号や観測データを表す非負値行列（以下、観測行列という）から当該行列を潜在的に構成する非負値行列群を推定する問題（非負値行列分解問題）が汎用的に利用されている。例えば、音声スペクトログラムまたは画像を入力として混在するノイズを除去する問題や、アンケート等の統計集約データを入力として内在するユーザー特性を分解して抽出する問題（トピックモデル）で実用的に利用されつつある。そして、こうした応用事例において処理時間の軽減やリアルタイム性が求められる場面も増えてきている。
【０００３】
ここでは、まず、上記非負値行列分解問題の定式化を行う。なお、明細書の記載表記の制約上、次のような表記に関する約束を導入する。^（キャレット）は上付き添字を表すこととし、例えば、xy^zはyzがxに対する上付き添字であることを表す。（アンダースコア）は下付き添字を表すこととし、例えば、xy_zはyzがxに対する上付き添字であることを表す。また、_ある文字xに対する^xや~xのような上付き添え字の”^”や”~”は、本来”x”の真上に記載されるべきであるが、^xや~xと記載することとする。
【０００４】
（非負の）観測行列Y∈RI×Jは、非負値行列分解可能な行列Z∈RI×Jとノイズを表す行列E∈RI×Jの和で構成されるとする。ここで、非負値行列Zは、例えば、ノイズを除去した音声スペクトログラムまたはノイズを除去した画像を表す行列となる。また、非負値行列Zは次式のように2つの非負値行列{A, B}（ただし、A∈RI×K, B∈RJ×K）の積で表されることとする。
【０００５】
【数１】
2020126441.tif 000003
【０００６】
ここで、Kは十分小さな数であり、KはI, Jより小さい数であるのが好ましい。また、Tは転置を表す。なお、行列{A, B}を分解要素行列と呼ぶこともある。
【０００７】
また、コスト関数J(A, B):RI×K×RJ×K→R∪{-∞, +∞}は、損失項Gと行列{A, B}に対する正規化項HA, HBの和で構成される。
【０００８】
【数２】
2020126441.tif 000004
【０００９】
ここで、損失項Gは2つ行列（式(1-2)では、行列Yと行列ABT）の誤差を測定するものである。また、正規化項HA, HBは、行列{A, B}の非負性を保証するものである。
【００１０】
以上より、解くべき問題は、コスト関数Jを最小化するような行列{A, B}を推定することであり、次式により定式化できる。
【００１１】
【数３】
2020126441.tif 000005
【００１２】
ここで、式(1-3)はコスト関数Jの下界(infimum)となる非負値行列{A, B}を推定する操作を表す。
【００１３】
以下、損失項と正規化項について詳細に説明する。まず、損失項について説明する。損失項Gには、任意の閉真凸関数（以後、単に凸関数と表す）を利用できる。例えば、（狭義の凸関数である）ブレグマンダイバージェンスによって損失項Gを表す場合、以下のようになる。
【００１４】
【数４】
2020126441.tif 000006
【００１５】
例えば、式(1-5)で表される（連続で微分可能かつ強凸である）狭義凸関数Φ(β)を用いることにすると、損失項GΦ^(β)として式(1-6)のβダイバージェンスを得る。βダイバージェンスは非負行列分解でよく利用される損失項の一つである。
【００１６】
【数５】
2020126441.tif 000007
【００１７】
なお、このβダイバージェンスは、β=1のときフロベニウスノルム、β=0のとき一般化KLダイバージェンス、β=-1のとき板倉−斎藤距離と一致する。
【００１８】
次に、正規化項について説明する。正規化項に対しては、少なくとも変数として扱う行列{A, B}が非負であることを保証する必要がある。このために、例えば、行列の全要素が0以上である場合に0を出力、その他の場合には+∞を出力する非負バリア関数を用いて正規化項を定義することができる。
【００１９】
【数６】
2020126441.tif 000008
【００２０】
ただし、
【００２１】
【数７】
2020126441.tif 000009
【００２２】
(0≦i≦I, 0≦j≦J, 0≦k≦K)である。
【００２３】
さらに、行列{A, B}の確率的な性質として、行列{A, B}がスパース性を持つことを仮定とすると、L1ノルムを用いた次式の正規化項とすることもできる。
【００２４】
【数８】
2020126441.tif 000010
【００２５】
ただし、μA, μBは事前に与えられた（正の）定数である。
【００２６】
続いて、上記非負値行列分解問題を解くための、非負値行列分解の最適化アルゴリズムについて説明する。ここでは、乗法更新法(multiplicative update rule)（非特許文献１）に基づく最適化アルゴリズムと交互方向乗数法(Alternating direction method of multipliers, ADMM)（非特許文献２）に基づく最適化アルゴリズムについて説明する。
【００２７】
まず、乗法更新法に基づく最適化アルゴリズムについて説明する。乗法更新法は、前向き勾配法の一種として表現できるものであり、現在の要素からコスト関数Jの勾配にステップサイズηを掛け合わせた量を減算するという処理で構成される。
【００２８】
【数９】
2020126441.tif 000011
【００２９】
例えば、式(1-6)のβダイバージェンスを損失項として用いる場合、コスト関数Jの勾配は次式で与えられる。
【００３０】
【数１０】
2020126441.tif 000012
【００３１】
ここで、ステップサイズηik, ηjk,を次式で定義すると、
【００３２】
【数１１】
2020126441.tif 000013
【００３３】
式(1-10)及び式(1-11)は以下のようになる。
【００３４】
【数１２】
2020126441.tif 000014
【００３５】
ここで、[z]+=max{z, 0}であり、ε(>0)は定数である。
【００３６】
したがって、乗法更新法に基づく最適化アルゴリズムでは、観測行列Yを入力とし、式(1-14)及び式(1-15)により行列{A, B}を所定の回数更新することにより、行列{A, B}を求める。なお、式(1-1)を用いることにより、ノイズを表す行列Eも求めることができる。
【００３７】
次に、ADMMに基づく最適化アルゴリズムについて説明する。なお、以下説明する方法ではコスト関数の形式が非特許文献２のそれとは異なるため、アルゴリズムの形式に差異があるものとなるが、本質的には等価なものである。
【００３８】
このアルゴリズムにおいて重要なアイデアは、式(1-2)の損失項に含まれる行列積ABTを補助変数Zとして置き換え、制約付凸最小化問題として扱う点にある。
【００３９】
【数１３】
2020126441.tif 000015
【００４０】
式(1-16)に関連するラグランジュ関数Lを次式で定義する。
【００４１】
【数１４】
2020126441.tif 000016
【００４２】
ただし、U∈RI×Jは双対変数である。また、<U, Z>=tr(UTZ)である。
【００４３】
ここで、（主問題を解く代わりに）双対問題を解くことにすると、非負値行列分解問題は次式の最適化問題として定式化される。
【００４４】
【数１５】
2020126441.tif 000017
【００４５】
ここで、式(1-18)はラグランジュ関数Lの下界(infimum)となる行列Z, A, Bを推定した後、推定した行列Z, A, Bに対するラグランジュ関数Lの上界(supremum)となる双対変数Uを推定する操作を表す。
【００４６】
式(1-18)の問題を解く、ADMMに基づく最適化アルゴリズムは図１、図２のようになる。図１の最適化アルゴリズムと図２の最適化アルゴリズムは等価なものである。また、図１・図２の最適化アルゴリズムに現れる変数{~X, ~V}は、双対変数Uを非線形変換して得られる変数である。この変数{~X, ~V}を用いることにより、アルゴリズム中の数式が単純化される。κ(>0)はステップサイズを表す定数である。図１・図２の最適化アルゴリズムは、いずれも、観測行列Yを入力とし、行列{A, B}を所定の回数更新することにより、行列{A, B}を求める。
【先行技術文献】
【非特許文献】
【００４７】
【非特許文献１】D. D. Lee and H. S. Seung, “Learning the parts of objects by non-negative matrix factorization”, Nature, Vol.401, No.6755, pp.788-791, 1999.
【非特許文献２】D. L. Sun and C. Fevotte, “Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence”, in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) 2014. IEEE, pp.6201-6205, 2014.
【発明の概要】
【発明が解決しようとする課題】
【００４８】
乗法更新法に基づく最適化アルゴリズムは、更新則が非常に単純で、実装が容易であるというメリットがある。また、収束速度も更新の初期段階ではとても速くかつ１ステップあたりの計算量も低いというメリットもある。しかし、更新回数が増えて不動点（すなわち、コスト関数を最小化する解）に近づくに従い、誤差の収束率が悪化することが報告されている（非特許文献２参照）。つまり、乗法更新法に基づく最適化アルゴリズムには、推定精度の高い解に近づきにくく、収束が遅いという問題がある。
【００４９】
一方、ADMMに基づく最適化アルゴリズムは、更新回数が増加しても、誤差の収束率が悪化しにくい傾向が実験的に示されている（非特許文献２参照）。しかし、更新に必要となるステップサイズ（図１・図２のアルゴリズムにおけるκ）を事前に決める必要があり、適切にステップサイズを設定しないと、その効果が得られないことが報告されている。つまり、ADMMに基づく最適化アルゴリズムには、ステップサイズが適切に設定されない場合など推定精度の高い解に近づきにくい場合があり、収束が安定しないという問題がある。
【００５０】
そこで本発明では、高速で安定的に収束する非負値行列分解技術を提供することを目的とする。
【課題を解決するための手段】
【００５１】
本発明の一態様は、Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化部を含む非負値行列分解最適化装置であって、前記最適化部が最適化に用いるコスト関数Jは、式J(A, B)=G(Y||ABT)+HA(A)+HB(B)（ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項）で定義されるものであり、前記最適化部は、ブレグマン単調作用素分解に基づいて行列{A, B}を最適化する。
【発明の効果】
【００５２】
本発明によれば、非負値行列を２つの非負値行列の積として表す非負値行列分解を高速かつ安定的に得ることができる。
【図面の簡単な説明】
【００５３】
【図１】ADMMを用いた最適化アルゴリズムの一例を示す図である。
【図２】ADMMを用いた最適化アルゴリズムの一例を示す図である。
【図３】本発明の実施形態における最適化アルゴリズム（B-P-R型）の一例を示す図である。
【図４】本発明の実施形態における最適化アルゴリズム（B-D-R型）の一例を示す図である。
【図５】非負値行列分解最適化装置１００の構成を示すブロック図である。
【図６】非負値行列分解最適化装置１００の動作を示すフローチャートである。
【図７】最適化部１１０の構成を示すブロック図である。
【図８】最適化部１１０の動作を示すフローチャートである。
【発明を実施するための形態】
【００５４】
以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【００５５】
各実施形態の説明に先立って、この明細書における表記方法について説明する。
【００５６】
^（キャレット）は上付き添字を表す。例えば、xy^zはyzがxに対する上付き添字であり、xy^zはyzがxに対する下付き添字であることを表す。また、_（アンダースコア）は下付き添字を表す。例えば、xy_zはyzがxに対する上付き添字であり、xy_zはyzがxに対する下付き添字であることを表す。
【００５７】
ある文字xに対する^xや~xのような上付き添え字の”^”や”~”は、本来”x”の真上に記載されるべきであるが、明細書の記載表記の制約上、^xや~xと記載しているものである。
【００５８】
＜技術的背景＞
ADMMではDouglas-Rachford型単調作用素分解を用いている。また、先述したように、ADMMに基づく最適化アルゴリズムはステップサイズなど所定のパラメータを事前に適切に設定する必要がある。
【００５９】
上記２点に着目し、本発明の実施形態における最適化アルゴリズムを、ADMMに基づく最適化アルゴリズムをベースとした高速かつ安定的に収束する最適化アルゴリズムとして構成する。具体的には、本発明の実施形態における最適化アルゴリズムを、（１）ステップサイズなどのパラメータ調整を不要とし、（２）Peaceman-Rachford型単調作用素分解も利用可能とした最適化アルゴリズムとして構成する。
【００６０】
まず、本発明の実施形態における最適化アルゴリズムの背景となる理論について説明する。式(1-18)で表される双対問題は、以下のように式変形することにより、凸共役関数（閉真凸関数）の和の最小化問題として扱える。
【００６１】
【数１６】
2020126441.tif 000018
【００６２】
ここで、G*, H*は次式で表される凸共役関数（閉真凸関数）である。
【００６３】
【数１７】
2020126441.tif 000019
【００６４】
式(2-1)の不動点は、次式のように凸共役関数{G*, H*}の劣微分の和が零行列を含むときに得られる。
【００６５】
【数１８】
2020126441.tif 000020
【００６６】
ここで、Ti:RI×J→RI×J(i=1, 2)は凸共役関数{G*, H*}の劣微分作用素を表す。すなわち、
【００６７】
【数１９】
2020126441.tif 000021
【００６８】
となり、Tiは、I×J行列を入力とし、I×J行列を出力する単調作用素となる。
【００６９】
パラメータ調整を不要とするためのメインとなるアイデアは、凸関数に適応して変数空間の計量を自動的に修正することにある。ここでは、そのための作用素∇D-1を導入する。ただし、関数D:RI×J→RI×Jは（連続で微分可能かつ強凸である）狭義凸関数であり、その微分作用素∇D:RI×J→RI×Jとその逆作用素∇D-1:RI×J→RI×Jは次式を満たすものとする。
【００７０】
【数２０】
2020126441.tif 000022
【００７１】
凸関数に適応して変数空間の計量を自動的に修正する方法については後ほど詳述する。
【００７２】
関数Dが式(2-5)を満たすとき、不動点に関する条件である式(2-4)の両辺に∇D-1を適用することで以下を得る。
【００７３】
【数２１】
2020126441.tif 000023
【００７４】
ここで、○は作用素の合成を表す。したがって、式(2-6)では、∇D-1とTiを合成している。
【００７５】
参考非特許文献１に従って式(2-6)を変形すると、式(2-7)のBregman Peaceman-Rachford (B-P-R)型単調作用素分解と式(2-8)のBregman Douglas-Rachford (B-D-R)型単調作用素分解の2種類の単調作用素分解が得られる。
【００７６】
【数２２】
2020126441.tif 000024
【００７７】
ここで、X∈RI×Jは双対変数Uの補助変数、ξは平均化作用素の平均化係数である。なお、ξにはξ∈(0, 1)を満たす任意の数を用いることができるが、1/2を用いることが多い。
（参考非特許文献１：K. Niwa and W. B. Kleijn, “Bregman monotone operator splitting”, arXiv preprint arXiv:1807.04871, 2018.）
【００７８】
式(2-7)、式(2-8)をみるとわかるように、B-D-R型単調作用素分解はB-P-R型単調作用素分解に平均化作用素を適用したものであるため、B-P-R型単調作用素分解の方が高速に収束することが予測される。
【００７９】
式(2-7)及び式(2-8)から再帰的な変数更新則を導出すると、式(2-9)〜式(2-13)のような更新処理が導出される。
【００８０】
【数２３】
2020126441.tif 000025
【００８１】
式(2-9)〜式(2-12)が式(2-7)から導出されるB-P-R型変数更新則、式(2-9)〜式(2-11)と式(2-13)が式(2-8)から導出されるB-D-R型変数更新則である。
【００８２】
式(2-9)〜式(2-13)には凸共役関数の劣微分作用素Tiなどが含まれているため、更なる式変更を行う。まず、式(2-9)を主変数と双対変数を交互に更新するように式変形を行う。式(2-9)に含まれる凸共役関数G*（式(2-2)参照）に関連する問題を以下で定義する。
【００８３】
【数２４】
2020126441.tif 000026
【００８４】
上記問題のラグランジュ関数Qを次式で定義する。
【００８５】
【数２５】
2020126441.tif 000027
【００８６】
ラグランジュ関数Qの劣微分がゼロベクトルを含むように変数{Z, U}は更新される。
【００８７】
【数２６】
2020126441.tif 000028
【００８８】
であることから、以下を得る。
【００８９】
【数２７】
2020126441.tif 000029
【００９０】
この変数Zと変数Uの関係を利用すると、式(2-9)のレゾルヴェント作用素R1の入出力変数{U, X}に対して、以下が成り立つ。
【００９１】
【数２８】
2020126441.tif 000030
【００９２】
この式を変形することで、以下を得る。
【００９３】
【数２９】
2020126441.tif 000031
【００９４】
これは以下の式を解くことと等価である。
【００９５】
【数３０】
2020126441.tif 000032
【００９６】
更新された変数Zを用いて、変数Uも更新される。
【００９７】
【数３１】
2020126441.tif 000033
【００９８】
式(2-10)と上記式を組み合わせ、変数Uを除去することで、双対補助変数Vの更新則が得られる。
【００９９】
【数３２】
2020126441.tif 000034
【０１００】
式を簡略化するために、以下の双対補助変数{~V, ~X}を導入する。
【０１０１】
【数３３】
2020126441.tif 000035
【０１０２】
双対補助変数{~V, ~X}を用いると、式(2-9)及び式(2-10)は以下のように表される。
【０１０３】
【数３４】
2020126441.tif 000036
【０１０４】
ここで、ブレグマンダイバージェンスBD^+は変数空間の計量を修正する狭義凸関数であり、以下で与えられる。
【０１０５】
【数３５】
2020126441.tif 000037
【０１０６】
同様に、式(2-11)〜式(2-13)は以下のように表される。
【０１０７】
【数３６】
2020126441.tif 000038
【０１０８】
式(2-14)〜式(2-19)が本発明の実施形態の最適化アルゴリズムの変数更新則を構成する式である。図３、図４に、B-P-R型の最適化アルゴリズム、B-D-R型の最適化アルゴリズムを示す。
【０１０９】
最後に、パラメータ調整をしなくても高速に収束することを実現するための狭義凸関数D+を設計する方法について説明する。具体的には、凸関数Gに合わせて関数D+を適応的に修正する方法を説明する。ここでは、一つの実装方法として、次式の一般化二次形式を用いる。なお、計算量低減のために、式(1-4)のように行列要素間の相関を無視したコストを想定している。
【０１１０】
【数３７】
2020126441.tif 000039
【０１１１】
ここで、zi,jは行列Zの(i,j)要素を表す。ε(>0)は十分に小さい数であり、関数Dが狭義凸関数であることを保証するために用いる。
【０１１２】
また、∇D+(0)=0を満たすことから、hijは、次式のように、関数Gの二次勾配に従うように設計するのが効率的である。
【０１１３】
【数３８】
2020126441.tif 000040
【０１１４】
ただし、zijoldは現在の行列Zの(i,j)要素を表す。
【０１１５】
このとき、関数D+と微分作用素が逆となる関数Dは、次式により与えられる。
【０１１６】
【数３９】
2020126441.tif 000041
【０１１７】
さらに計算量を削減し、高速化するためには、hijが次式の行列Hの(i,j)要素となるようにするとよい。
【０１１８】
【数４０】
2020126441.tif 000042
【０１１９】
ここで、hA∈RI, hB∈RJである。
【０１２０】
式(2-23)は（I×J）個の行列要素が（I+J）個のベクトル要素で表されることを示しており、変数更新の高速化に加えて、メモリ量抑制にも寄与する。
【０１２１】
なお、ベクトル{hA, hB}は、例えば、所定の数の関数Gの二次勾配の平均値を利用して決定することができる。
【０１２２】
なお、ここでは具体的な実験データは示さないが、図３・図４の最適化アルゴリズムを用いることにより、従来の最適化アルゴリズム（例えば、図１・図２の最適化アルゴリズム）と比べて、行列{A, B}の収束速度、行列{A, B}の推定精度がともに向上することが確認された。特に、図３のB-P-R型の最適化アルゴリズムで効果が高いことが確認された。
【０１２３】
＜第１実施形態＞
Yを観測信号または観測データを表す非負の観測行列、Zを式(1-1)を満たす観測行列Yを構成する非負の行列とする。非負値行列分解最適化装置１００は、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である非負の行列{A, B}を最適化する。
【０１２４】
以下、図５〜図６を参照して非負行列分解最適化装置１００を説明する。図５は、非負行列分解最適化装置１００の構成を示すブロック図である。図６は、非負行列分解最適化装置１００の動作を示すフローチャートである。図５に示すように非負行列分解最適化装置１００は、最適化部１１０と、記録部１９０を含む。記録部１９０は、非負行列分解最適化装置１００の処理に必要な情報を適宜記録する構成部である。記録部１９０は、例えば、入力となる観測行列Y、最適化対象となる行列{A, B}を記録する。
【０１２５】
図６に従い非負行列分解最適化装置１００の動作について説明する。
【０１２６】
Ｓ１１０において、最適化部１１０は、観測行列Yを入力として、Z=ABTを満たす行列Zの分解である行列{A, B}を最適化する。具体的には、最適化部１１０は、観測行列Yを入力として、コスト関数Jの最小化問題を解くことにより行列{A, B}を最適化する。ここで、コスト関数J(A, B):RI×K×RJ×K→R∪{-∞, +∞}は、次式で定義される。
【０１２７】
【数４１】
2020126441.tif 000043
【０１２８】
（ただし、Gは損失項、HAは行列Aに対する正規化項、HBは行列Bに対する正規化項である。）
【０１２９】
以下、図７〜図８を参照して最適化部１１０について説明する。図７は、最適化部１１０の構成を示すブロック図である。図８は、最適化部１１０の動作を示すフローチャートである。図７に示すように最適化部１１０は、初期化部１１１、行列Z更新部１１２と、第１双対変数更新部１１３と、行列A更新部１１４と、行列B更新部１１５と、第２双対変数更新部１１６と、カウンタ更新部１１７と、終了条件判定部１１８を含む。
【０１３０】
図８に従い最適化部１１０の動作について説明する。
【０１３１】
Ｓ１１１において、初期化部１１１は、カウンタtを初期化する。具体的には、t=1とする。また、初期化部１１１は、行列A、行列B、行列Z、双対変数~Xを初期化する。
【０１３２】
Ｓ１１２において、行列Z更新部１１２は、次式により、行列Zを更新する。
【０１３３】
【数４２】
2020126441.tif 000044
【０１３４】
（ただし、HZは行列Zに対する正規化項、BD^+は関数D+を用いて定義されるブレグマンダイバージェンスである。）
【０１３５】
なお、ブレグマンダイバージェンスBD^+の定義に用いる関数D+は、式
【０１３６】
【数４３】
2020126441.tif 000045
【０１３７】
（ただし、Cはcijを(i, j)要素とする行列、hij=∇ij2G(cijold)（cijoldは現在の行列Cの(i, j) 要素、ε(>0)は所定の定数）により与えられる。また、関数D+は、hijを(i, j)要素とする行列Hに対してH=hAhBTを満たすベクトルhA, hBを用いて計算してもよい。関数D+を計算する際、Zや~Xが用いられる。
【０１３８】
この関数D+は、行列A更新部１１４、行列B更新部１１５でも用いる。この場合、関数D+を計算する際、ABTや-~Vが用いられる。
【０１３９】
Ｓ１１３において、第１双対変数更新部１１３は、次式により、双対変数~Vを更新する。
【０１４０】
【数４４】
2020126441.tif 000046
【０１４１】
Ｓ１１４において、行列A更新部１１４は、次式により、行列Aを更新する。
【０１４２】
【数４５】
2020126441.tif 000047
【０１４３】
Ｓ１１５において、行列B更新部１１５は、次式により、行列Bを更新する。
【０１４４】
【数４６】
2020126441.tif 000048
【０１４５】
Ｓ１１６において、第２双対変数更新部１１６は、次式により、双対変数~Xを更新する。
【０１４６】
【数４７】
2020126441.tif 000049
【０１４７】
Ｓ１１７において、カウンタ更新部１１７は、カウンタtを1だけインクリメントする。具体的には、t←t+1とする。
【０１４８】
Ｓ１１８において、終了条件判定部１１８は、カウンタtが所定の更新回数T（Tは1以上の整数であり、例えば10万）に達した場合（つまり、t>Tとなり、終了条件が満たされた場合）は、そのときの行列{A, B}を出力して、処理を終了する。それ以外の場合、Ｓ１１２の処理に戻る。つまり、最適化部１１０は、Ｓ１１２〜Ｓ１１８の処理を繰り返す。
【０１４９】
なお、Ｓ１１４とＳ１１５（つまり、行列Aの更新処理と行列Bの更新処理）については、Ｓ１１４、Ｓ１１５の順に実行してもよいし、Ｓ１１５、Ｓ１１４の順に実行してもよい。また、Ｓ１１４とＳ１１５を並列に実行するようにしてもよい。
【０１５０】
（変形例）
第２双対変数更新部１１６では、B-P-R型更新則により双対変数~Xを更新したが、B-D-R型更新則により更新するようにしてもよい。この場合、Ｓ１１６において、第２双対変数更新部１１６は、次式により、双対変数~Xを更新する。
【０１５１】
【数４８】
2020126441.tif 000050
【０１５２】
（ただし、ξは0<ξ<1を満たす定数）
【０１５３】
以上、まとめると、最適化部１１０は、ブレグマン単調作用素分解に基づいて行列{A, B}を最適化するものである。
【０１５４】
本実施形態の発明によれば、非負値行列を２つの非負値行列の積として表す非負値行列分解を高速かつ安定的に得ることができる。また、凸関数Gに合わせて関数D+を適応的に修正するため、変数空間の計量を適応的に修正することができ、従来のADMMに基づく最適化アルゴリズムで必要であった（例えば、ステップサイズなどの）パラメータの調整が不要となる。
【０１５５】
＜補記＞
本発明の装置は、例えば単一のハードウェアエンティティとして、キーボードなどが接続可能な入力部、液晶ディスプレイなどが接続可能な出力部、ハードウェアエンティティの外部に通信可能な通信装置（例えば通信ケーブル）が接続可能な通信部、ＣＰＵ（Central Processing Unit、キャッシュメモリやレジスタなどを備えていてもよい）、メモリであるＲＡＭやＲＯＭ、ハードディスクである外部記憶装置並びにこれらの入力部、出力部、通信部、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部記憶装置の間のデータのやり取りが可能なように接続するバスを有している。また必要に応じて、ハードウェアエンティティに、ＣＤ−ＲＯＭなどの記録媒体を読み書きできる装置（ドライブ）などを設けることとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。
【０１５６】
ハードウェアエンティティの外部記憶装置には、上述の機能を実現するために必要となるプログラムおよびこのプログラムの処理において必要となるデータなどが記憶されている（外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくこととしてもよい）。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶される。
【０１５７】
ハードウェアエンティティでは、外部記憶装置（あるいはＲＯＭなど）に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてメモリに読み込まれて、適宜にＣＰＵで解釈実行・処理される。その結果、ＣＰＵが所定の機能（上記、…部、…手段などと表した各構成要件）を実現する。
【０１５８】
本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、上記実施形態において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【０１５９】
既述のように、上記実施形態において説明したハードウェアエンティティ（本発明の装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。
【０１６０】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【０１６１】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【０１６２】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記憶装置に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１６３】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【０１６４】
上述の本発明の実施形態の記載は、例証と記載の目的で提示されたものである。網羅的であるという意思はなく、開示された厳密な形式に発明を限定する意思もない。変形やバリエーションは上述の教示から可能である。実施形態は、本発明の原理の最も良い例証を提供するために、そして、この分野の当業者が、熟考された実際の使用に適するように本発明を色々な実施形態で、また、色々な変形を付加して利用できるようにするために、選ばれて表現されたものである。すべてのそのような変形やバリエーションは、公正に合法的に公平に与えられる幅にしたがって解釈された添付の請求項によって定められた本発明のスコープ内である。
【図１】
2020126441.tif 000076
【図２】
2020126441.tif 000077
【図３】
2020126441.tif 000078
【図４】
2020126441.tif 000079
【図５】
2020126441.tif 000080
【図６】
2020126441.tif 000081
【図７】
2020126441.tif 000082
【図８】
2020126441.tif 000083