1-P-8
Frobenius ノルム規準の非負値行列因子分解における乗法更新式に関する一考察∗○亀岡弘和, ルルージョナトン
日本電信電話 (株) NTT コミュニケーション科学基礎研究所

1序論

目的関数 G(z) に対し，非負値行列因子分解 (以後，NMF) とは，与えられたG(z) = min G+ (z, z̄)

非負値行列 Y = (ym,n )M ×N ∈ R≥0,M ×N に対し，z̄(5)が成り立つとき，G+ (z, z̄) を G(z) の補助関数，z̄ を補Y ≃ HU(1)

助変数と定義する。このとき，次の定理が成り立つ。となるような非負値行列 H = (hm,k )M ×K ∈ R≥0,M ×K ,定理 1 (補助関数法). 補助関数 G+ (z, z̄) を，z̄ に関し

U = (uk,n )K×N ∈ R≥0,K×N を決定することをいい，非負信号のブラインド信号分離への応用技術として脳・画て最小化するステップと，z に関して最小化するステップを繰り返すと，目的関数値は単調収束する。

像・音響信号処理などの分野で注目されている。
式 (1) の解法として，Y と HU の近さを誤差行列 Y −HU の Frobenius ノルムで測った場合の最適化問題
minimizesubject tof (H, U) = ∥Y − HU∥2FH ∈ R≥0,M ×K , U ∈ R≥0,K×N= argminz̄ G+ (zt , z̄) と し ，zt+1 =argminz G (z, z̄t+1 ) とする。ただし，t は反復計算のステップ数を表す。z = zt , z̄ = z̄t から z =証明： z̄t+1+(2)zt+1 , z̄ = z̄t+1 に更新されたときに，G(z) が増加しないことを示す。明らかに G(zt ) = G+ (zt , z̄t+1 ) であを解く効率的な反復アルゴリズムが Lee らによって与えり，zt+1 = argminz G+ (z, z̄t+1 ) より，G+ (zt , z̄t+1 ) ≥られている [1]。このアルゴリズムは，[]
YUT m,k]hm,k ← hm,k [HUUT m,k[ T ]H Y k,n]uk,n ← uk,n [ TH HU k,n(3)

G+ (zt+1 , z̄t+1 ) で あ る ．更 に 補 助 関 数 の 定 義 よ り，
G+ (zt+1 , z̄t+1 ) ≥ G(zt+1 ) なのだから，結局，G(zt ) ≥G(zt+1 ) である。
以上を踏まえ，式 (4) を導出する。ここで，yn :=(4)(y1,n , · · · , yM,n )T と H が与えられたとき，g(un ) :=∥yn − Hun ∥22のような乗算型の更新式からなるため，乗法更新法と呼ばれることがある。[·]i , [·]i,j はそれぞれベクトルの i 番(6)を un := (u1,n , · · · , uK,n )T に関して最小化する問題を目の要素，行列の {i, j} 番目の要素を表すものとする。 考えよう。式 (6) を直接最小化する un は M ≥ K のと
本報告では，Lee らによって導かれたこの反復アルゴ き un = (HT H)−1 HT yn で与えられるが，この場合 un
リズム (以後，Lee-Seung(LS) アルゴリズム)[1] がなぜの各要素が非負性を満たすかどうかは定かではない。LS効果的であるかについて一考察を与え，これをもとに， アルゴリズムのポイントは，g(un ) に対し要素同士の交より効率的な反復アルゴリズムが得られる可能性につい差項が存在しない二次形式の補助関数を立てる点にあて議論する。また，この議論を通して，複素 NMF[2] にり，これにより各要素について非負値制約に注意しながおける最適化アルゴリズムの効率化の可能性についてもら独立に更新式を求めることができる。このような補助言及する。

関数は次のようにして立てられる。まず，

2LS アルゴリズム [1]以下，LS アルゴリズムの導出の概要を示す。H, U はともに変数であるため，式 (2) は非負制約つき非線形最適化問題であり，解析的に解くことはできない。これにg + (un , u′n ) :=g(u′n ) + (un − u′n )T ∇g(u′n )1+ (un − u′n )T W(un − u′n )2と置く。ここで，W が対し，LS アルゴリズムは，制約を満たしながら目的関数を反復的に降下する手法であり，次に示す補助関数法の原理に基づく。
∗(7)

W = diag(w1 , · · · , wK )[] /wk = 2 HT Hu′n k u′k,n(8)(9)

Study on Multiplicative Update Formula for Nonnegative Matrix Factorization with Frobenius Norm Criterion. by Hirokazu Kameoka, Jonathan Le Roux (NTT Communication Science Laboratories)

日本音響学会講演論文集- 709 -2009年9月

で与えられるとき，を用いて無数に作ることができる。上記不等式の証明に(k)g(un ) ≤ g+(un , u′n )ついては [2] を参照されたい。この不等式は，βm の値(10)(k)に依らず αm が([] )]1 [(k)(k)が言える (W−2HT H が半正定値であることによる [1])。
αm=hm,k uk + βmym − Hu mym
明らかに u′n = un のとき g(un ) = g + (un , u′n ) であるた
め，g + (un , u′n ) は g(un ) の補助関数としての定義を満 のとき等号成立 [2] するため，式 (16) の右辺)2∑ 1 (たす。従って，定理 1 より，g + (un , u′n ) を u′n と un に(k)G + (u, α) :=αy−humm,kkm(k)
関してそれぞれ最小化するステップ
m,k βm(17)(18)

u′n ← un(11)は g(u) の補助関数としての定義を満たし，αm が補助un ← u′n − W−1 ∇g(u′n )(12)変数となる。G + (u, α) はベクトル要素 uk 同士の交差項(k)のない二次関数の和になっており，このような性質の補を行えば g(un ) を降下させることができる。式 (11), (12)をまとめて，さらに要素ごとの表現に書き直すと，[ T ]H yn k]uk,n ← uk,n [ T(13)H Hun k(k)助関数が βm の任意性の分，無数に定義できたことになる。ここで注目すべきは，2 章で立てた補助関数が，hm,k uk(k)]βm=[Hu m(19)となり，式 (4) を得る。g + (un , u′n ) は un の要素 uk,n ごと置いた場合の G + (u, α) の特殊形に相当している点でとの二次関数に分離されており，式 (13) の更新値はそある。これを以下で確認しよう。の放物線の頂点の座標に該当する。従って，式 (13) がも式 (17) に式 (19) を代入するとし負であれば uk,n ≥ 0 の制約下で g+(un , u′n ) が最小となるのは uk,n = 0 のときである。このため，正確には{}[ T ]H yn k] ,0uk,n ← max uk,n [ T(14)H Hun kと書かれるべきであるが，H, yn , un の要素がすべて非負値であれば式 (13) の更新では un は必ず非負値に更新されるため，H と U の初期値が非負であることが前提となっている場合は更新式の表現は式 (13) のままで十分である。このように更新値の非負性が逐次的に保証される点は LS アルゴリズムの特徴の一つである。式 (3)も以上と同様の考え方により得られる。

3hm,k uk(k)]αm=[Hu m(20)(k)(k)が得られ，βm = αm となる。式 (19), (20) での uk は式 (18) の uk と区別する必要があるため u′k と置いた上で式 (19), (20) を式 (18) に代入すると)∑ hm,k u′ (] 2u [[] k ym − k′ Hu′ m
G + (u, α) =uk
Hu′ m(21)m,kを得る。これは，式 (7) を要素ごとの表記に書き直したものと等しい。なお，ここでは H を定数と見なしている
ので α を補助変数と見なすことは u′ を補助変数と見なすことと同じである。以上より，G + (u, α) は g + (u, u′ )を包含する一般形であることが分かった。

一般化 LS アルゴリズム

さて，定理 1 より，以下，同様に un に関する更新法についてのみ議論す(k)αm← argmin G + (u, α)る。また，表記の簡単化のため n は以後省略する。

(22)(k)αm以上で見てきたように，LS アルゴリズムのポイント

uk ← argmin G + (u, α)は，g(u) を直接最小化しようとするのではなく，u の(23)uk

要素同士の交差項が存在せず要素ごとの二次関数に分離を行えば g(u) を降下させることができる。ここで，式された形態をもつ補助関数を設計することにあった。実(22) は式 (17) で与えられ，式 (23) は∑(k)(k)
αmym hm,k /βmは，以上のような性質の補助関数は，∑(k)αm= 1,k∑(k)(k)
βm= 1, βm>0(15)uk =m∑k(k)(24)m(k)を満たす任意の αm , βm のもとで成り立つ不等式)2∑ 1 ((k)αy−hug(u) ≤mm,kkm(k)m,k βm(k)h2m,k /βm(16)で与えられる。hm,k の更新方法の一般形についても以上と同様にして得られる。hm,k の更新式を含めて式 (22),(23) により構成される反復アルゴリズムを一般化 LS アルゴリズムと呼ぶことにする。

日本音響学会講演論文集- 710 -2009年9月

LS アルゴリズムの正当性について

4を式 (27) に代入すると(k)3 章では LS アルゴリズムは一般化 LS アルゴリズムに(k)おいて βm を式 (19) と置いた特殊ケースに相当することを示した。ここでは，その合理性について考察する。
(k)
まず，βm の役割について考えよう。式 (18) を見る(k)と明らかなように，βm は，どの u においても g(u) を決して下回らない G + (u, α) の選び方の中で，uk ごとのαm(k)(k)βm=∑= αm(k)αmkとなる。これを式 (17) に代入すると，最終的にhm,k uk(k)(k)]αm= βm=[Hu m二次関数の頂点の座標を調節する変数と見ることができる。従って，式 (23) において G + (u, α) をどれだけ降下(29)(30)が導かれる。

(k)

できるかは，βm の選び方次第となる。以上より，uk の以上より，LS アルゴリズムは，一般化 LS アルゴリ二次関数の頂点の値ができるだけ小さくなるように βm(k)ズムの枠組から見たとき，各ステップにおいて目的関数を決めることができれば，反復アルゴリズムを効率化での降下量を最大にする βm が選択されているわけではきるはずである。

ないが，式 (26) のある種の近似解が選択されているこ(k)

以上の観点に立ち，各ステップで βm を最適に決定(k)とが明らかになった。このことは，LS アルゴリズムがする問題を定式化する。G (u, α) は，効率的に動作することに対する一つの正当性を示すも{(のである。一方で，LS アルゴリズムはどの場面でも効+G (u, α) =+∑)
uk(k)mk−2∑ hm,k 2βm)((k)
∑ αmym hm,kuk +(k)βmm

率的であるとは限らず，より安定した効率性を備えた反2(k) 2∑ αmym 2}

復アルゴリズムを検討する余地が残されていることが分かった。

(k)
βmmと書けるので，G (u, α) の u に関して最小となる値は+(k)βm の関数となり，{(k) 2∑ ∑ αmym 2
F(β) :=(k)βmmk()2/()}(k)∑ αm∑ hm,k 2ym hm,k−(k)(k)βmmm βm以上の考察をヒントにし，ここでは LS アルゴリズムの改良の可能性について議論する。
前章で示したとおり，各ステップで目的関数の降下量(k)を最大にする βm は解析的に求めることはできないが，(k)(25)

最適な βm を数値的に探索することはできる。すなわ
ち，一般化 LS アルゴリズムにおいて，式 (22) の後段で行われる式 (23) の更新を u, β に関する最適化問題と捉(k)で与えられる。βm には式 (15) の制約があるため，minimize F(β)∑(k)
subject toβm
= 1 for all m

改良 LS アルゴリズムの可能性について5えるわけである。ここで，G + (u, α) を u に関して最小化するステップと β に関して最小化するステップを交互(26)に繰り返すアプローチをとることにすると，反復アルゴリズムは表 1 のとおりとなる。なお，式 (31) は u を固kを解く問題に帰着する。この最適化問題は解析的に解く

定した状態で式 (15) の制約下において G + (u, α) を最小ことはできないが，実は式 (19) の結果は，F(β) の第一にする β の解である。以下に示したアルゴリズムはあく項を最小化する解となっている。これを以下で確認しよまで改良案であり，真に改良アルゴリズムとなっているう。式 (15) の制約に対応するラグランジュ未定乗数 λかどうかは数値実験で評価する必要がある。これは今後を導入し，ラグランジュ関数をの検討課題である。

L(β, λ) :=(k) 2∑ ∑ αmym 2(k)km+βm∑mλm(∑)(k)
βm−1k(k)と定義する。これを βm に関して偏微分し，0 と置くと(k)αm ym(k)βm= √λmが得られるので，式 (15) に代入し，()2(k)∑ αm∑ym2(k)√= 1 ⇒ λm = ymαmλmkk

日本音響学会講演論文集

複素 NMF アルゴリズム [2] の効率化

6

本稿での議論は，Frobenius ノルム規準の NMF に限らず，我々が以前提案した複素 NMF[2] の問題設定にも当てはまる。以下，複素 NMF の概要を示す。

(27)

複素 NMF とは，観測信号の時間周波数成分 Y :=(Yω,t )Ω×T ∈ CΩ×T に対し，以下に示すモデル
Fω,t =(28)- 711 -

∑

(k)

(k)

Hω(k) Ut ejϕω,t(32)k

2009年9月

表1をもとにした反復アルゴリズムを提案したが，β をどの改良 LS アルゴリズム (案)ように決定すべきかが重要課題として残されていた。ここでは本稿の考察に従ってこの課題を以下で議論する。

1. H, U を初期設定する。

複素 NMF における β の決定方法において，5 章まで2. β を以下のとおりに初期設定する。の議論との唯一の違いは複素数を扱っている点である。
この点に注意して，H, U , ϕ の更新による目的関数の降hm,k uk,n(k)]βm,n=[HU m,n下量を最大にする β の最適化問題を 3 章と同様に定式化する。例えば U の更新においては，β の目的関数は{¯ (k) ¯2 ¯¯
3. α を以下の式に従って更新する。
¯ ¯Yω,t ¯2∑ ∑ ¯αm([] )] F(β) :=1 [(k)(k)(k)←ym,n − HU m,nαm,nhm,k uk,n + βm,nβω,tωk,tym,n)/(()}]2∑ Hω(k) 2∑ Hω(k) [ (k)(k)∗ ∗ jϕω,t−Re αω,t Yω,t e
4. 以下の 2 ステップを収束するまで繰り返す。
(k)(k)βω,tωω βω,t(a) H, U の行列要素を以下の式に従って更新する。となる。F(β) を最小化する β は解析的に求めることは∑/(k)(k)できないが，4 章と同様 F(β) の第一項を最小化するには
αm,nym,n uk,n βm,nn¯ (k) ¯∑hm,k ←
¯αω,t ¯/ (k)2(k)
uk,n βm,nβω,t = ∑ ¯ (k) ¯(35)¯αω,t ¯n∑
/ (k)(k)kym,n hm,k βm,nαm,nuk,n ←n∑のように更新すれば良いことが示される。この β の更新/ (k)hm,k 2 βm,n則は，[2] で採用した更新則と異なるため，今後それぞnれの更新則を採用した場合のアルゴリズムの振る舞いを(b) β を以下の式に従って更新する。
¯ (k)¯¯αm,n ym,n − hm,k uk,n ¯(k)βm,n ← ∑ ¯¯(k)¯αm,nym,n − hm,k uk,n ¯

比較検討したい。

(31)

まとめ

7

本報告では，Frobenius ノルム規準の非負値行列因子k分解 (NMF) のための反復アルゴリズム (LS アルゴリズ

収束後，2. に戻る。

ム)[1] を包含するアルゴリズムの一般形 (一般化 LS アルゴリズム) を示し，それにより導かれる解釈をもとに，を，U ができるだけスパースになるようにフィッティンより効率的なアルゴリズム (改良 LS アルゴリズム) をグするものである。各パラメータの説明については [2]実現できる可能性について考察した。また，その考察をを参照されたい。複素 NMF では，具体的には，もとに，一般化 LS アルゴリズムと類似のアプローチに

f (H, U, ϕ) :=

∑¯
∑ (k)
¯
¯Yω,t − Fω,t ¯2 + 2λ|Ut |pω,tよりパラメータの最適化を行う複素 NMF [2] における(33)

参考文献を目的関数としてminimize
subject to

計算の効率化に向けて，その展望について議論した。

k,tf (H, U, ϕ)∑Hω(k) = 1 (k = 1, · · · , K)

(34)

x

[1] D.D. Lee and H.S. Seung, “Algorithms for nonnegative matrix factorization,” in Advances inNeural Information Processing Systems 13 (Proc.
NIPS*2000), pp. 556–562, 2000.のような最適化問題を解くことになる。[2] ではこの最適化問題の解法として，補助関数[2] H. Kameoka, N. Ono, K. Kashino, S. Sagayama,f + (H, U, ϕ, α, V )¯∑ ∑ 1 ¯¯ (k)(k)(k) (k) jϕω,t ¯ 2:=Y−HUe¯α¯ω,tω,ttω(k)ω,t k βω,t}∑{(k)(k) 2(k)+λp|Vt |p−2 Ut+ (2 − p)|Vt |p”Complex NMF: A New Sparse Representation forAcoustic Signals,” In Proc. 2009 IEEE International Conference on Acoustics, Speech and Signal
Processing (ICASSP2009), pp. 3437-3440, 2009.

k,t

日本音響学会講演論文集- 712 -2009年9月