特許7167746 | 知財ポータル「IP Force」

知財求人 - 知財ポータルサイト「IP Force」

青山学院大学 (神奈川県相模原市中央区淵野辺)

▶ 日本電信電話株式会社の特許一覧

特許7167746非負値行列分解最適化装置、非負値行列分解最適化方法、プログラム

書誌要約請求の範囲詳細な説明課題実施例実施するための形態図面の説明

目に優しい文字サイズ小中大 PDF Top

< >

1
2
3
4
5
6
7
8

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2022-10-31

(45)【発行日】2022-11-09

(54)【発明の名称】非負値行列分解最適化装置、非負値行列分解最適化方法、プログラム

(51)【国際特許分類】

G06N 99/00 20190101AFI20221101BHJP

G06F 17/16 20060101ALI20221101BHJP

【ＦＩ】

G06N99/00 180

G06F17/16 Z

【請求項の数】 7

(21)【出願番号】P 2019018425

(22)【出願日】2019-02-05

(65)【公開番号】P2020126441

(43)【公開日】2020-08-20

【審査請求日】2021-05-28

(73)【特許権者】

【識別番号】000004226

【氏名又は名称】日本電信電話株式会社

(74)【代理人】

【識別番号】100121706

【弁理士】

【氏名又は名称】中尾直樹

(74)【代理人】

【識別番号】100128705

【弁理士】

【氏名又は名称】中村幸雄

(74)【代理人】

【識別番号】100147773

【弁理士】

【氏名又は名称】義村宗洋

(72)【発明者】

【氏名】丹羽健太

(72)【発明者】

【氏名】原田登

【審査官】杉浦孝光

(56)【参考文献】

【文献】特開２０１４－２０３４５９（ＪＰ，Ａ）

【文献】CHRETIEN, Stephane , et al.，ESTIMATING FEATURES WITH MISSING VALUES AND OUTLIERS: A BREGMAN-PROXIMAL POINT ALGORITHM FOR ROBUST，arXiv.org，2015年02月09日，[online], [検索日 2020.02.26]，インターネット：<URL: https://arxiv.org/abs/1502.02523>

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０６Ｎ９９／００

Ｇ０６Ｆ１７／１６

(57)【特許請求の範囲】

【請求項1】

Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、
観測行列Yを入力として、Z=AB^Tを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化部を含む非負値行列分解最適化装置であって、
前記最適化部が最適化に用いるコスト関数Jは、式

【数49】

（ただし、Gは損失項、H_Aは行列Aに対する正規化項、H_Bは行列Bに対する正規化項）
で定義されるものであり、
前記最適化部は、
次式により、行列Zを更新する行列Z更新部と、

【数50】

（ただし、H_Zは行列Zに対する正規化項、B_D^+は関数D⁺を用いて定義されるブレグマンダイバージェンス）
次式により、双対変数~Vを更新する第１双対変数更新部と、

【数51】

次式により、行列Aを更新する行列A更新部と、

【数52】

次式により、行列Bを更新する行列B更新部と、

【数53】

次式により、双対変数~Xを更新する第２双対変数更新部と、

【数54】

を含む非負値行列分解最適化装置。

【請求項2】

【数55】

【数56】

【数57】

次式により、行列Aを更新する行列A更新部と、

【数58】

次式により、行列Bを更新する行列B更新部と、

【数59】

次式により、双対変数~Xを更新する第２双対変数更新部と、

【数60】

（ただし、ξは0<ξ<1を満たす定数）
を含む非負値行列分解最適化装置。

【請求項3】

請求項１または２に記載の非負値行列分解最適化装置であって、
前記行列Z更新部、前記行列A更新部、前記行列B更新部で用いるブレグマンダイバージェンスB_D^+の定義に用いる関数D⁺は、式

【数61】

（ただし、Cはc_ijを(i, j)要素とする行列、h_ij=∇_ij ²G(c_ij ^old)（c_ij ^oldは現在の行列Cの(i, j) 要素、ε(>0)は所定の定数）
により与えられる
ことを特徴とする非負値行列分解最適化装置。

【請求項4】

請求項３に記載の非負値行列分解最適化装置であって、
関数D⁺は、h_ijを(i, j)要素とする行列Hに対してH=h_Ah_B ^Tを満たすベクトルh_A, h_Bを用いて計算される
ことを特徴とする非負値行列分解最適化装置。

【請求項5】

Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、
行列Z更新部と、第１双対変数更新部と、行列A更新部と、行列B更新部と、第２双対変数更新部とを含む非負値行列分解最適化装置が、観測行列Yを入力として、Z=AB^Tを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化ステップを実行する非負値行列分解最適化方法であって、
前記最適化ステップで最適化に用いるコスト関数Jは、式

【数62】

（ただし、Gは損失項、H_Aは行列Aに対する正規化項、H_Bは行列Bに対する正規化項）
で定義されるものであり、
前記最適化ステップは、
前記行列Z更新部が、次式により、行列Zを更新する行列Z更新ステップと、

【数63】

（ただし、H_Zは行列Zに対する正規化項、B_D^+は関数D⁺を用いて定義されるブレグマンダイバージェンス）
前記第１双対変数更新部が、次式により、双対変数~Vを更新する第１双対変数更新ステップと、

【数64】

前記行列A更新部が、次式により、行列Aを更新する行列A更新ステップと、

【数65】

前記行列B更新部が、次式により、行列Bを更新する行列B更新ステップと、

【数66】

前記第２双対変数更新部が、次式により、双対変数~Xを更新する第２双対変数更新ステップと、

【数67】

を含む非負値行列分解最適化方法。

【請求項6】

【数68】

【数69】

【数70】

前記行列A更新部が、次式により、行列Aを更新する行列A更新ステップと、

【数71】

前記行列B更新部が、次式により、行列Bを更新する行列B更新ステップと、

【数72】

前記第２双対変数更新部が、次式により、双対変数~Xを更新する第２双対変数更新ステップと、

【数73】

（ただし、ξは0<ξ<1を満たす定数）
を含む非負値行列分解最適化方法。

【請求項7】

請求項１ないし４のいずれか１項に記載の非負値行列分解最適化装置としてコンピュータを機能させるためのプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、非負値行列を２つの非負値行列の積として表す非負値行列分解技術に関する。

【背景技術】

【0002】

観測信号や観測データを表す非負値行列（以下、観測行列という）から当該行列を潜在的に構成する非負値行列群を推定する問題（非負値行列分解問題）が汎用的に利用されている。例えば、音声スペクトログラムまたは画像を入力として混在するノイズを除去する問題や、アンケート等の統計集約データを入力として内在するユーザー特性を分解して抽出する問題（トピックモデル）で実用的に利用されつつある。そして、こうした応用事例において処理時間の軽減やリアルタイム性が求められる場面も増えてきている。

【0003】

ここでは、まず、上記非負値行列分解問題の定式化を行う。なお、明細書の記載表記の制約上、次のような表記に関する約束を導入する。^（キャレット）は上付き添字を表すこととし、例えば、x^{y^z}はy^zがxに対する上付き添字であることを表す。（アンダースコア）は下付き添字を表すこととし、例えば、x^y_zはy_zがxに対する上付き添字であることを表す。また、_ある文字xに対する^xや~xのような上付き添え字の”^”や”~”は、本来”x”の真上に記載されるべきであるが、^xや~xと記載することとする。

【0004】

（非負の）観測行列Y∈R^I×Jは、非負値行列分解可能な行列Z∈R^I×Jとノイズを表す行列E∈R^I×Jの和で構成されるとする。ここで、非負値行列Zは、例えば、ノイズを除去した音声スペクトログラムまたはノイズを除去した画像を表す行列となる。また、非負値行列Zは次式のように2つの非負値行列{A, B}（ただし、A∈R^I×K, B∈R^J×K）の積で表されることとする。

【0005】

【数1】

【0006】

ここで、Kは十分小さな数であり、KはI, Jより小さい数であるのが好ましい。また、Tは転置を表す。なお、行列{A, B}を分解要素行列と呼ぶこともある。

【0007】

また、コスト関数J(A, B):R^I×K×R^J×K→R∪{-∞, +∞}は、損失項Gと行列{A, B}に対する正規化項H_A, H_Bの和で構成される。

【0008】

【数2】

【0009】

ここで、損失項Gは2つ行列（式(1-2)では、行列Yと行列AB^T）の誤差を測定するものである。また、正規化項H_A, H_Bは、行列{A, B}の非負性を保証するものである。

【0010】

以上より、解くべき問題は、コスト関数Jを最小化するような行列{A, B}を推定することであり、次式により定式化できる。

【0011】

【数3】

【0012】

ここで、式(1-3)はコスト関数Jの下界(infimum)となる非負値行列{A, B}を推定する操作を表す。

【0013】

以下、損失項と正規化項について詳細に説明する。まず、損失項について説明する。損失項Gには、任意の閉真凸関数（以後、単に凸関数と表す）を利用できる。例えば、（狭義の凸関数である）ブレグマンダイバージェンスによって損失項Gを表す場合、以下のようになる。

【0014】

【数4】

【0015】

例えば、式(1-5)で表される（連続で微分可能かつ強凸である）狭義凸関数Φ^(β)を用いることにすると、損失項G_Φ^(β)として式(1-6)のβダイバージェンスを得る。βダイバージェンスは非負行列分解でよく利用される損失項の一つである。

【0016】

【数5】

【0017】

なお、このβダイバージェンスは、β=1のときフロベニウスノルム、β=0のとき一般化KLダイバージェンス、β=-1のとき板倉－斎藤距離と一致する。

【0018】

次に、正規化項について説明する。正規化項に対しては、少なくとも変数として扱う行列{A, B}が非負であることを保証する必要がある。このために、例えば、行列の全要素が0以上である場合に0を出力、その他の場合には+∞を出力する非負バリア関数を用いて正規化項を定義することができる。

【0019】

【数6】

【0020】

ただし、

【0021】

【数7】

【0022】

(0≦i≦I, 0≦j≦J, 0≦k≦K)である。

【0023】

さらに、行列{A, B}の確率的な性質として、行列{A, B}がスパース性を持つことを仮定とすると、L1ノルムを用いた次式の正規化項とすることもできる。

【0024】

【数8】

【0025】

ただし、μ_A, μ_Bは事前に与えられた（正の）定数である。

【0026】

続いて、上記非負値行列分解問題を解くための、非負値行列分解の最適化アルゴリズムについて説明する。ここでは、乗法更新法(multiplicative update rule)（非特許文献１）に基づく最適化アルゴリズムと交互方向乗数法(Alternating direction method of multipliers, ADMM)（非特許文献２）に基づく最適化アルゴリズムについて説明する。

【0027】

まず、乗法更新法に基づく最適化アルゴリズムについて説明する。乗法更新法は、前向き勾配法の一種として表現できるものであり、現在の要素からコスト関数Jの勾配にステップサイズηを掛け合わせた量を減算するという処理で構成される。

【0028】

【数9】

【0029】

例えば、式(1-6)のβダイバージェンスを損失項として用いる場合、コスト関数Jの勾配は次式で与えられる。

【0030】

【数10】

【0031】

ここで、ステップサイズη_ik, η_jk,を次式で定義すると、

【0032】

【数11】

【0033】

式(1-10)及び式(1-11)は以下のようになる。

【0034】

【数12】

【0035】

ここで、[z]₊=max{z, 0}であり、ε(>0)は定数である。

【0036】

したがって、乗法更新法に基づく最適化アルゴリズムでは、観測行列Yを入力とし、式(1-14)及び式(1-15)により行列{A, B}を所定の回数更新することにより、行列{A, B}を求める。なお、式(1-1)を用いることにより、ノイズを表す行列Eも求めることができる。

【0037】

次に、ADMMに基づく最適化アルゴリズムについて説明する。なお、以下説明する方法ではコスト関数の形式が非特許文献２のそれとは異なるため、アルゴリズムの形式に差異があるものとなるが、本質的には等価なものである。

【0038】

このアルゴリズムにおいて重要なアイデアは、式(1-2)の損失項に含まれる行列積AB^Tを補助変数Zとして置き換え、制約付凸最小化問題として扱う点にある。

【0039】

【数13】

【0040】

式(1-16)に関連するラグランジュ関数Lを次式で定義する。

【0041】

【数14】

【0042】

ただし、U∈R^I×Jは双対変数である。また、<U, Z>=tr(U^TZ)である。

【0043】

ここで、（主問題を解く代わりに）双対問題を解くことにすると、非負値行列分解問題は次式の最適化問題として定式化される。

【0044】

【数15】

【0045】

ここで、式(1-18)はラグランジュ関数Lの下界(infimum)となる行列Z, A, Bを推定した後、推定した行列Z, A, Bに対するラグランジュ関数Lの上界(supremum)となる双対変数Uを推定する操作を表す。

【0046】

式(1-18)の問題を解く、ADMMに基づく最適化アルゴリズムは図１、図２のようになる。図１の最適化アルゴリズムと図２の最適化アルゴリズムは等価なものである。また、図１・図２の最適化アルゴリズムに現れる変数{~X, ~V}は、双対変数Uを非線形変換して得られる変数である。この変数{~X, ~V}を用いることにより、アルゴリズム中の数式が単純化される。κ(>0)はステップサイズを表す定数である。図１・図２の最適化アルゴリズムは、いずれも、観測行列Yを入力とし、行列{A, B}を所定の回数更新することにより、行列{A, B}を求める。

【先行技術文献】

【非特許文献】

【0047】

【文献】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.

【発明の概要】

【発明が解決しようとする課題】

【0048】

乗法更新法に基づく最適化アルゴリズムは、更新則が非常に単純で、実装が容易であるというメリットがある。また、収束速度も更新の初期段階ではとても速くかつ１ステップあたりの計算量も低いというメリットもある。しかし、更新回数が増えて不動点（すなわち、コスト関数を最小化する解）に近づくに従い、誤差の収束率が悪化することが報告されている（非特許文献２参照）。つまり、乗法更新法に基づく最適化アルゴリズムには、推定精度の高い解に近づきにくく、収束が遅いという問題がある。

【0049】

一方、ADMMに基づく最適化アルゴリズムは、更新回数が増加しても、誤差の収束率が悪化しにくい傾向が実験的に示されている（非特許文献２参照）。しかし、更新に必要となるステップサイズ（図１・図２のアルゴリズムにおけるκ）を事前に決める必要があり、適切にステップサイズを設定しないと、その効果が得られないことが報告されている。つまり、ADMMに基づく最適化アルゴリズムには、ステップサイズが適切に設定されない場合など推定精度の高い解に近づきにくい場合があり、収束が安定しないという問題がある。

【0050】

そこで本発明では、高速で安定的に収束する非負値行列分解技術を提供することを目的とする。

【課題を解決するための手段】

【0051】

本発明の一態様は、Yを観測信号または観測データを表す非負の観測行列、Zを観測行列Yを構成する非負の行列とし、観測行列Yを入力として、Z=AB^Tを満たす行列Zの分解である非負の行列{A, B}を最適化する最適化部を含む非負値行列分解最適化装置であって、前記最適化部が最適化に用いるコスト関数Jは、式J(A, B)=G(Y||AB^T)+H_A(A)+H_B(B)（ただし、Gは損失項、H_Aは行列Aに対する正規化項、H_Bは行列Bに対する正規化項）で定義されるものであり、前記最適化部は、ブレグマン単調作用素分解に基づいて行列{A, B}を最適化する。

【発明の効果】

【0052】

本発明によれば、非負値行列を２つの非負値行列の積として表す非負値行列分解を高速かつ安定的に得ることができる。

【図面の簡単な説明】

【0053】

【図1】ADMMを用いた最適化アルゴリズムの一例を示す図である。

【図2】ADMMを用いた最適化アルゴリズムの一例を示す図である。

【図3】本発明の実施形態における最適化アルゴリズム（B-P-R型）の一例を示す図である。

【図4】本発明の実施形態における最適化アルゴリズム（B-D-R型）の一例を示す図である。

【図5】非負値行列分解最適化装置１００の構成を示すブロック図である。

【図6】非負値行列分解最適化装置１００の動作を示すフローチャートである。

【図7】最適化部１１０の構成を示すブロック図である。

【図8】最適化部１１０の動作を示すフローチャートである。

【発明を実施するための形態】

【0054】

以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。

【0055】

各実施形態の説明に先立って、この明細書における表記方法について説明する。

【0056】

^（キャレット）は上付き添字を表す。例えば、x^{y^z}はy^zがxに対する上付き添字であり、x_y^zはy^zがxに対する下付き添字であることを表す。また、_（アンダースコア）は下付き添字を表す。例えば、x^y_zはy_zがxに対する上付き添字であり、x_{y_z}はy_zがxに対する下付き添字であることを表す。

【0057】

ある文字xに対する^xや~xのような上付き添え字の”^”や”~”は、本来”x”の真上に記載されるべきであるが、明細書の記載表記の制約上、^xや~xと記載しているものである。

【0058】

＜技術的背景＞
ADMMではDouglas-Rachford型単調作用素分解を用いている。また、先述したように、ADMMに基づく最適化アルゴリズムはステップサイズなど所定のパラメータを事前に適切に設定する必要がある。

【0059】

上記２点に着目し、本発明の実施形態における最適化アルゴリズムを、ADMMに基づく最適化アルゴリズムをベースとした高速かつ安定的に収束する最適化アルゴリズムとして構成する。具体的には、本発明の実施形態における最適化アルゴリズムを、（１）ステップサイズなどのパラメータ調整を不要とし、（２）Peaceman-Rachford型単調作用素分解も利用可能とした最適化アルゴリズムとして構成する。

【0060】

まず、本発明の実施形態における最適化アルゴリズムの背景となる理論について説明する。式(1-18)で表される双対問題は、以下のように式変形することにより、凸共役関数（閉真凸関数）の和の最小化問題として扱える。

【0061】

【数16】

【0062】

ここで、G^*, H^*は次式で表される凸共役関数（閉真凸関数）である。

【0063】

【数17】

【0064】

式(2-1)の不動点は、次式のように凸共役関数{G^*, H^*}の劣微分の和が零行列を含むときに得られる。

【0065】

【数18】

【0066】

ここで、T_i:R^I×J→R^I×J(i=1, 2)は凸共役関数{G^*, H^*}の劣微分作用素を表す。すなわち、

【0067】

【数19】

【0068】

となり、T_iは、I×J行列を入力とし、I×J行列を出力する単調作用素となる。

【0069】

パラメータ調整を不要とするためのメインとなるアイデアは、凸関数に適応して変数空間の計量を自動的に修正することにある。ここでは、そのための作用素∇D^-1を導入する。ただし、関数D:R^I×J→R^I×Jは（連続で微分可能かつ強凸である）狭義凸関数であり、その微分作用素∇D:R^I×J→R^I×Jとその逆作用素∇D^-1:R^I×J→R^I×Jは次式を満たすものとする。

【0070】

【数20】

【0071】

凸関数に適応して変数空間の計量を自動的に修正する方法については後ほど詳述する。

【0072】

関数Dが式(2-5)を満たすとき、不動点に関する条件である式(2-4)の両辺に∇D^-1を適用することで以下を得る。

【0073】

【数21】

【0074】

ここで、○は作用素の合成を表す。したがって、式(2-6)では、∇D^-1とT_iを合成している。

【0075】

参考非特許文献１に従って式(2-6)を変形すると、式(2-7)のBregman Peaceman-Rachford (B-P-R)型単調作用素分解と式(2-8)のBregman Douglas-Rachford (B-D-R)型単調作用素分解の2種類の単調作用素分解が得られる。

【0076】

【数22】

【0077】

ここで、X∈R^I×Jは双対変数Uの補助変数、ξは平均化作用素の平均化係数である。なお、ξにはξ∈(0, 1)を満たす任意の数を用いることができるが、1/2を用いることが多い。
（参考非特許文献１：K. Niwa and W. B. Kleijn, “Bregman monotone operator splitting”, arXiv preprint arXiv:1807.04871, 2018.）

【0078】

式(2-7)、式(2-8)をみるとわかるように、B-D-R型単調作用素分解はB-P-R型単調作用素分解に平均化作用素を適用したものであるため、B-P-R型単調作用素分解の方が高速に収束することが予測される。

【0079】

式(2-7)及び式(2-8)から再帰的な変数更新則を導出すると、式(2-9)～式(2-13)のような更新処理が導出される。

【0080】

【数23】

【0081】

式(2-9)～式(2-12)が式(2-7)から導出されるB-P-R型変数更新則、式(2-9)～式(2-11)と式(2-13)が式(2-8)から導出されるB-D-R型変数更新則である。

【0082】

式(2-9)～式(2-13)には凸共役関数の劣微分作用素T_iなどが含まれているため、更なる式変更を行う。まず、式(2-9)を主変数と双対変数を交互に更新するように式変形を行う。式(2-9)に含まれる凸共役関数G^*（式(2-2)参照）に関連する問題を以下で定義する。

【0083】

【数24】

【0084】

上記問題のラグランジュ関数Qを次式で定義する。

【0085】

【数25】

【0086】

ラグランジュ関数Qの劣微分がゼロベクトルを含むように変数{Z, U}は更新される。

【0087】

【数26】

【0088】

であることから、以下を得る。

【0089】

【数27】

【0090】

この変数Zと変数Uの関係を利用すると、式(2-9)のレゾルヴェント作用素R₁の入出力変数{U, X}に対して、以下が成り立つ。

【0091】

【数28】

【0092】

この式を変形することで、以下を得る。

【0093】

【数29】

【0094】

これは以下の式を解くことと等価である。

【0095】

【数30】

【0096】

更新された変数Zを用いて、変数Uも更新される。

【0097】

【数31】

【0098】

式(2-10)と上記式を組み合わせ、変数Uを除去することで、双対補助変数Vの更新則が得られる。

【0099】

【数32】

【0100】

式を簡略化するために、以下の双対補助変数{~V, ~X}を導入する。

【0101】

【数33】

【0102】

双対補助変数{~V, ~X}を用いると、式(2-9)及び式(2-10)は以下のように表される。

【0103】

【数34】

【0104】

ここで、ブレグマンダイバージェンスB_D^+は変数空間の計量を修正する狭義凸関数であり、以下で与えられる。

【0105】

【数35】

【0106】

同様に、式(2-11)～式(2-13)は以下のように表される。

【0107】

【数36】

【0108】

式(2-14)～式(2-19)が本発明の実施形態の最適化アルゴリズムの変数更新則を構成する式である。図３、図４に、B-P-R型の最適化アルゴリズム、B-D-R型の最適化アルゴリズムを示す。

【0109】

最後に、パラメータ調整をしなくても高速に収束することを実現するための狭義凸関数D⁺を設計する方法について説明する。具体的には、凸関数Gに合わせて関数D⁺を適応的に修正する方法を説明する。ここでは、一つの実装方法として、次式の一般化二次形式を用いる。なお、計算量低減のために、式(1-4)のように行列要素間の相関を無視したコストを想定している。

【0110】

【数37】

【0111】

ここで、z_i,jは行列Zの(i,j)要素を表す。ε(>0)は十分に小さい数であり、関数Dが狭義凸関数であることを保証するために用いる。

【0112】

また、∇D⁺(0)=0を満たすことから、h_ijは、次式のように、関数Gの二次勾配に従うように設計するのが効率的である。

【0113】

【数38】

【0114】

ただし、z_ij ^oldは現在の行列Zの(i,j)要素を表す。

【0115】

このとき、関数D⁺と微分作用素が逆となる関数Dは、次式により与えられる。

【0116】

【数39】

【0117】

さらに計算量を削減し、高速化するためには、h_ijが次式の行列Hの(i,j)要素となるようにするとよい。

【0118】

【数40】

【0119】

ここで、h_A∈R^I, h_B∈R^Jである。

【0120】

式(2-23)は（I×J）個の行列要素が（I+J）個のベクトル要素で表されることを示しており、変数更新の高速化に加えて、メモリ量抑制にも寄与する。

【0121】

なお、ベクトル{h_A, h_B}は、例えば、所定の数の関数Gの二次勾配の平均値を利用して決定することができる。

【0122】

なお、ここでは具体的な実験データは示さないが、図３・図４の最適化アルゴリズムを用いることにより、従来の最適化アルゴリズム（例えば、図１・図２の最適化アルゴリズム）と比べて、行列{A, B}の収束速度、行列{A, B}の推定精度がともに向上することが確認された。特に、図３のB-P-R型の最適化アルゴリズムで効果が高いことが確認された。

【0123】

＜第１実施形態＞
Yを観測信号または観測データを表す非負の観測行列、Zを式(1-1)を満たす観測行列Yを構成する非負の行列とする。非負値行列分解最適化装置１００は、観測行列Yを入力として、Z=AB^Tを満たす行列Zの分解である非負の行列{A, B}を最適化する。

【0124】

以下、図５～図６を参照して非負行列分解最適化装置１００を説明する。図５は、非負行列分解最適化装置１００の構成を示すブロック図である。図６は、非負行列分解最適化装置１００の動作を示すフローチャートである。図５に示すように非負行列分解最適化装置１００は、最適化部１１０と、記録部１９０を含む。記録部１９０は、非負行列分解最適化装置１００の処理に必要な情報を適宜記録する構成部である。記録部１９０は、例えば、入力となる観測行列Y、最適化対象となる行列{A, B}を記録する。

【0125】

図６に従い非負行列分解最適化装置１００の動作について説明する。

【0126】

Ｓ１１０において、最適化部１１０は、観測行列Yを入力として、Z=AB^Tを満たす行列Zの分解である行列{A, B}を最適化する。具体的には、最適化部１１０は、観測行列Yを入力として、コスト関数Jの最小化問題を解くことにより行列{A, B}を最適化する。ここで、コスト関数J(A, B):R^I×K×R^J×K→R∪{-∞, +∞}は、次式で定義される。

【0127】

【数41】

【0128】

（ただし、Gは損失項、H_Aは行列Aに対する正規化項、H_Bは行列Bに対する正規化項である。）

【0129】

以下、図７～図８を参照して最適化部１１０について説明する。図７は、最適化部１１０の構成を示すブロック図である。図８は、最適化部１１０の動作を示すフローチャートである。図７に示すように最適化部１１０は、初期化部１１１、行列Z更新部１１２と、第１双対変数更新部１１３と、行列A更新部１１４と、行列B更新部１１５と、第２双対変数更新部１１６と、カウンタ更新部１１７と、終了条件判定部１１８を含む。

【0130】

図８に従い最適化部１１０の動作について説明する。

【0131】

Ｓ１１１において、初期化部１１１は、カウンタtを初期化する。具体的には、t=1とする。また、初期化部１１１は、行列A、行列B、行列Z、双対変数~Xを初期化する。

【0132】

Ｓ１１２において、行列Z更新部１１２は、次式により、行列Zを更新する。

【0133】

【数42】

【0134】

（ただし、H_Zは行列Zに対する正規化項、B_D^+は関数D⁺を用いて定義されるブレグマンダイバージェンスである。）

【0135】

なお、ブレグマンダイバージェンスB_D^+の定義に用いる関数D⁺は、式

【0136】

【数43】

【0137】

（ただし、Cはc_ijを(i, j)要素とする行列、h_ij=∇_ij ²G(c_ij ^old)（c_ij ^oldは現在の行列Cの(i, j) 要素、ε(>0)は所定の定数）により与えられる。また、関数D⁺は、h_ijを(i, j)要素とする行列Hに対してH=h_Ah_B ^Tを満たすベクトルh_A, h_Bを用いて計算してもよい。関数D⁺を計算する際、Zや~Xが用いられる。

【0138】

この関数D⁺は、行列A更新部１１４、行列B更新部１１５でも用いる。この場合、関数D⁺を計算する際、AB^Tや-~Vが用いられる。

【0139】

Ｓ１１３において、第１双対変数更新部１１３は、次式により、双対変数~Vを更新する。

【0140】

【数44】

【0141】

Ｓ１１４において、行列A更新部１１４は、次式により、行列Aを更新する。

【0142】

【数45】

【0143】

Ｓ１１５において、行列B更新部１１５は、次式により、行列Bを更新する。

【0144】

【数46】

【0145】

Ｓ１１６において、第２双対変数更新部１１６は、次式により、双対変数~Xを更新する。

【0146】

【数47】

【0147】

Ｓ１１７において、カウンタ更新部１１７は、カウンタtを1だけインクリメントする。具体的には、t←t+1とする。

【0148】

Ｓ１１８において、終了条件判定部１１８は、カウンタtが所定の更新回数T（Tは1以上の整数であり、例えば10万）に達した場合（つまり、t>Tとなり、終了条件が満たされた場合）は、そのときの行列{A, B}を出力して、処理を終了する。それ以外の場合、Ｓ１１２の処理に戻る。つまり、最適化部１１０は、Ｓ１１２～Ｓ１１８の処理を繰り返す。

【0149】

なお、Ｓ１１４とＳ１１５（つまり、行列Aの更新処理と行列Bの更新処理）については、Ｓ１１４、Ｓ１１５の順に実行してもよいし、Ｓ１１５、Ｓ１１４の順に実行してもよい。また、Ｓ１１４とＳ１１５を並列に実行するようにしてもよい。

【0150】

（変形例）
第２双対変数更新部１１６では、B-P-R型更新則により双対変数~Xを更新したが、B-D-R型更新則により更新するようにしてもよい。この場合、Ｓ１１６において、第２双対変数更新部１１６は、次式により、双対変数~Xを更新する。

【0151】

【数48】

【0152】

（ただし、ξは0<ξ<1を満たす定数）

【0153】

以上、まとめると、最適化部１１０は、ブレグマン単調作用素分解に基づいて行列{A, B}を最適化するものである。

【0154】

本実施形態の発明によれば、非負値行列を２つの非負値行列の積として表す非負値行列分解を高速かつ安定的に得ることができる。また、凸関数Gに合わせて関数D⁺を適応的に修正するため、変数空間の計量を適応的に修正することができ、従来のADMMに基づく最適化アルゴリズムで必要であった（例えば、ステップサイズなどの）パラメータの調整が不要となる。

【0155】

＜補記＞
本発明の装置は、例えば単一のハードウェアエンティティとして、キーボードなどが接続可能な入力部、液晶ディスプレイなどが接続可能な出力部、ハードウェアエンティティの外部に通信可能な通信装置（例えば通信ケーブル）が接続可能な通信部、ＣＰＵ（Central Processing Unit、キャッシュメモリやレジスタなどを備えていてもよい）、メモリであるＲＡＭやＲＯＭ、ハードディスクである外部記憶装置並びにこれらの入力部、出力部、通信部、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部記憶装置の間のデータのやり取りが可能なように接続するバスを有している。また必要に応じて、ハードウェアエンティティに、ＣＤ－ＲＯＭなどの記録媒体を読み書きできる装置（ドライブ）などを設けることとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。

【0156】

ハードウェアエンティティの外部記憶装置には、上述の機能を実現するために必要となるプログラムおよびこのプログラムの処理において必要となるデータなどが記憶されている（外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくこととしてもよい）。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶される。

【0157】

ハードウェアエンティティでは、外部記憶装置（あるいはＲＯＭなど）に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてメモリに読み込まれて、適宜にＣＰＵで解釈実行・処理される。その結果、ＣＰＵが所定の機能（上記、…部、…手段などと表した各構成要件）を実現する。

【0158】

本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、上記実施形態において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。

【0159】

既述のように、上記実施形態において説明したハードウェアエンティティ（本発明の装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。

【0160】

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ－ＲＡＭ（Random Access Memory）、ＣＤ－ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ－Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ－ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。

【0161】

また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ－ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

【0162】

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記憶装置に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

【0163】

また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

【0164】

上述の本発明の実施形態の記載は、例証と記載の目的で提示されたものである。網羅的であるという意思はなく、開示された厳密な形式に発明を限定する意思もない。変形やバリエーションは上述の教示から可能である。実施形態は、本発明の原理の最も良い例証を提供するために、そして、この分野の当業者が、熟考された実際の使用に適するように本発明を色々な実施形態で、また、色々な変形を付加して利用できるようにするために、選ばれて表現されたものである。すべてのそのような変形やバリエーションは、公正に合法的に公平に与えられる幅にしたがって解釈された添付の請求項によって定められた本発明のスコープ内である。

【図1】