1-P-5
I ダイバージェンス規準複素 NMF ∗☆鏡 英章 1 , 亀岡弘和 2 , △湯川正裕 1( 慶應大・総合デザイン工, 2 NTT CS 研)
1

1

はじめに

近年，非負値行列因子分解 (Non-negative MatrixFactorization; NMF[1]) に基づく音源分離手法が注目されている。観測信号の振幅（あるいはパワー）スペクトログラムを非負値行列と見なし，これを二つの非負値行列の積で近似することは，各時刻における観測スペクトルを複数のスペクトルテンプレートの重みつき和でフィッティングしていることに相当する。
従って，NMF により求まる二つの非負値行列を用いることで観測スペクトログラムを各スペクトルテンプレートに対応する音源成分に分解することができる。このアプローチではスペクトルの加法性を仮定しているが，振幅スペクトルやパワースペクトルは信号を非線形変換したものであるため，スペクトルは実際には加法的ではない。このため，NMF により観測スペクトログラムを加法的な成分に分解したとしても必ずしも各音源に対応したスペクトログラムが得られるとは限らない。
この課題を解決するため，NMF 法における各スペクトルテンプレートに位相スペクトルを時変パラメータとして付加することで，NMF 法と類似した考え方により複素スペクトル領域での信号分解を可能にする
複素 NMF[2] と呼ぶ枠組が提案されている。従来，複素 NMF におけるパラメータ推定問題（NMF 法における行列分解に相当）は，観測複素スペクトログラムと複素スペクトログラムモデルとの Euclid 距離（差の絶対値の二乗）を規準とした最適化問題として定式化されており，効率的なパラメータ推定アルゴリズムが提案されている。NMF では Euclid 距離，I ダイバージェンス（一般化 Kullback-Leibler ダイバージェンス），板倉齋藤擬距離，β ダイバージェンスを規準とした場合のパラメータ推定アルゴリズムがそれぞれ提案されているが，経験的に I ダイバージェンスを規準とした場合に高い音源分離性能が得られることが知られている [3, 4]。このことから，複素 NMF においても Euclid 距離以外の規準を用いた場合のパラメータ推定アルゴリズムを導出することができればより高い音源分離性能が得られる可能性がある。
本稿では NMF において I ダイバージェンス規準で音源分離性能が高い傾向にあるのと同様，提案するI ダイバージェンス規準複素 NMF により従来の複素NMF と比べて高い音源分離性能を得られることを確認した。

2既存手法

2.1 NMF 法
周波数と時刻に対応するインデックスを k =1, . . . , K, m = 1, . . . , M とし，観測信号の複素スペクトログラムを Yk,m ∈ C，Yk,m の絶対値またはその絶対値二乗を Zk,m ∈ R+ とする。NMF 法は観測振幅（またはパワー）スペクトログラム Z = (Zk,m )K×M ∈K×M
K×L
R+を非負行列 H ∈ R+, U ∈ RL×Mの積で+近似する手法で，Z の各要素をXZk,m ≃Hk,l Ul,m(1)l∗と近似することが目的となる。H の各列ベクトルをスペクトルテンプレートと解釈すると，式 (1) は，各時刻の観測スペクトルをスペクトルテンプレートの重みつき和で近似していることに相当する。ただし，lはスペクトルテンプレートに対応するインデックスである。式 (1) の分解表現が得られれば，例えば Wiener
フィルタ
Hk,l Ul,mX̂l,k,m = Yk,m P(2)
l Hk,l Ul,mにより各構成音の複素スペクトログラム X̂l,k,m を得ることができる。
式 (1) の分解表現を得る問題は Z と HU の誤差の規準を用いた最適化問題として定式化することができる。誤差規準としては Euclid 距離 DEU ，I ダイ
バージェンス DKLX2
DEU (Z | HU ) :=(Zk,m − (HU )k,m ) (3)k,mDKL (Z | HU )
X
Zk,m−Zk,m + (HU )k,m (4):=Zk,m log(HU )k,mk,mなどが用いられる。以上の目的関数を停留点へ収束させる H, U は次節で述べる補助関数法に基づく反復アルゴリズムにより探索することができる。ただし，(·)k,m は行列の k 行 m 列目の要素を表す。

2.2 補助関数法
目的関数を F とするときF (θ) = min F + (θ, θ̄)

(5)

θ̄

を満たす F + を F の補助関数，θ̄ を補助変数と呼ぶ。
定理 1 補助関数をθ̄ ← argmin F + (θ, θ̄), θ ← argmin F + (θ, θ̄) (6)θ̄θと更新するとき，目的関数 F は非増加である。
証明
θ̄ l+1 = argminθ̄ F + (θ l , θ̄), θ l+1 =+argminθ F (θ, θ̄ l+1 )と す る 。F (θ l )=minθ̄ F + (θ l , θ̄) = F + (θ l , θ̄ l+1 ) ≥ F + (θ l+1 , θ̄ l+1 ) ≥F (θ l+1 ). 
つまり式 (5) を満たすような補助関数を見つけることができれば，式 (6)(6) の更新を繰り返すことで目的関数を停留点へ収束させることができる。従来のNMF の反復アルゴリズムは上述の原理に基づいて導かれる。

2.3 複素 NMF
式 (1) では，振幅スペクトルまたはパワースペクトルが加法的であるという仮定の下で混合音スペクトルがモデル化されているが，振幅スペクトルもパワースペクトルも実際には加法的ではない。そこで，加法性が成り立つ複素スペクトル領域で NMF 法と同様の考え方による信号分解を可能にする枠組として提案されたのが複素 NMF[2] である。
前述の NMF のモデルにおける各スペクトルテン
プレートに時変の位相スペクトル ejφl,k,m ∈ C をパラ

Complex nonnegative matrix factorization with I divergence. by Hideaki Kagami (Keio University),
Hirokazu Kameoka (NTT Communication Science Laboratories), Masahiro Yukawa (Keio University).

日本音響学会講演論文集- 433 -2016年9月

メータとして付加することでX
Fk,m =Hk,l Ul,m ejφl,k,m(7)lのような形の複素スペクトログラムのモデルを立てることができる。観測信号の複素スペクトログラムを Y = (Yk,m )K×M ∈ CK×M とすると，複素 NMF
はこのモデルと観測複素スペクトログラムとの間のEuclid 距離XX
2IEU (H,U ,φ)=|Yk,m −Fk,m | +2λ|Ul,m |p (8)k,mStep 0) 初期値設定
Step 1) 補助変数の更新
Xl,k,m=Hk,l Ul,m ejφl,k,m+βl,k,m (Yk,m − Fk,m )(14)Vl,m = Ul,m(15)Step 2)

2|Yk,m − Fk,m |X 1
2Xl,k,m − Hk,l Ul,m ejφl,k,m(9)≤βl,k,mlPが成り立つ。ただし，l Xl,k,m = Yk,m , 0 < βl,k,m <P
1, l βl,k,m = 1 である。次に第二項に関し，0 < p ≤2 より，
22|Ul,m |p ≤ p|Vl,m |p−2 Ul,m+ (2 − p)|Vl,m |p (10)が成り立つ。ただし，Vl,m ≥ 0 である。以上より，
Xl,k,m ∈ C，Vl,m ∈ R を補助変数にもつ+IEU(H, U , φ, X) :=X
12Xl,k,m − Hk,l Ul,m ejφl,k,mβl,k,ml,k,mX
2+λ(p|Vl,m |p−2 Ul,m+ (2 − p)|Vl,m |p )Ul,m =(11)のような補助関数を設計することができる。また，Pβl,k,m は βl,k,m > 0, l βl,k,m = 1 を満たす任意定数である。以上の補助関数を最小にする補助変数 Xl,k,m ,Vl,m は，Hk,ljφl,k,m∗
k βl,k,m ℜ Xl,k,m e2P Hk,lp−2k βl,k,m +λp|Vl,m |[],Xl,k,mejφl,k,m = |Xl,k,m|

3(17)
(18)

問題の定式化

従来の複素 NMF では各構成音の位相スペクトログ
ラム ejφl,k,m はモデルパラメータ，Xl,k,m は各構成音の複素スペクトログラムの成分を意味する補助変数であった。ここで，ejφl,k,m と Xl,k,m に対する扱いを逆転させ，ejφl,k,m を補助変数，Xl,k,m をモデルパラメータと扱った場合の，上述の複素 NMF アルゴリズ
ムと等価な手続きを行う反復アルゴリズムについて考える。理由は後述するが，実はこの反復アルゴリズ
ムは次の最適化問題 (P0 )(P0 ) minimize JEU (H, U , X)
X
12(|Xl,k,m | − Hk,l Ul,m ):=
βl,k,ml,k,mX
subject to Yk,m =Xk,m(19)(20)

を停留点へ収束させるためのアルゴリズムとなっている。このことを以下で確認する。表記の簡単化のため cl,k,m = ejφl,k,m とすると，複素 NMF の補助関数+IEUの第一項は，X1(|Xl,k,m |2 − 2Hk,l Ul,m Re[c∗l,k,m Xl,k,m ]
βl,k,ml,k,m2
2+Hk,lUl,m)+βl,k,m (Yk,m −Fk,m ) (12)Vl,m = Ul,mP(16)

ll,mjφl,k,m

モデルパラメータの更新
PUl,mjφl,k,m∗]m βl,k,m ℜ[Xl,k,m e,Hk,l =U2
Pl,mm βl,k,ml,mを最小化する問題として定式化される。ただし，上式の第二項は U の要素がスパースになるように誘導する目的の正則化項であり，0 < p ≤ 2, λ > 0 とする。
この目的関数を最小にする H, U , φ を解析的に得ることはできないが，補助関数法により目的関数を停留点へ収束させる反復アルゴリズムを導くことができる。まず第一項に関し，二次関数は凸関数であるため，Jensen の不等式

Xl,k,m=Hk,l Ul,m e

Table 1 複素 NMF(式 (8)) 最小化アルゴリズム(13)となる。式 (12) は，各ステップで推定された l 番目の構成音に対応する複素スペクトログラムモデル
Hk,l Ul,m ejφl,k,m に全体のモデル化誤差 Yk,m −Fk,m のβl,k,m 倍を足したものを Xl,k,m と置いていることを意P味する。このことから，Xl,k,m は l Xl,k,m = Yk,mとなるように分配された l 番目の構成音に対応する複素スペクトログラムの成分と解釈することができる。同様に，以上の補助関数を最小にするパラメー
タ H, U , φ も解析的に求めることができる。なお，
Pβl,k,m > 0, l βl,k,m = 1 を満たす範囲であればβl,k,m を各ステップで更新したとしてもアルゴリズ
ムの収束性は保証されており，文献 [2] では βl,k,m =P Hk,l Ul,mHk,n Un,m のような更新式が採用されている。以上(21)と書ける。ここで，g(z) = −Re(c z) (|c| = 1) という関数と f (z) = −|z| という関数を Re(z) と Im(z) の 2変数の関数として見たとき，g(z) が f (z) に c の方向で接する接平面になっていることに注意すると，∗−Re(c∗ z) ≥ −|z|が成り立つ。従って，
X1(|Xl,k,m |2βl,k,m(22)l,k,m22− 2Hk,l Ul,m Re[c∗l,k,m Xl,k,m ] + Hk,lUl,m)X
1(|Xl,k,m |2≥
βl,k,m(23)l,k,mnより，複素 NMF アルゴリズムは Table 1 のようにまとめることができる。

22− 2Hk,l Ul,m |Xl,k,m | + Hk,lUl,m)X
1(|Xl,k,m | − Hk,l Ul,m )2=
βl,k,m(24)(25)

l,k,mが言え，cl,k,m を補助変数，Xl,k,m をモデルパラメー+タと見なせば，IEU
は最適化問題 (P0 ) の目的関数の日本音響学会講演論文集- 434 -

2016年9月

補助関数にもなっていることが分かる。従って，(P0 )+の停留点を探索するアルゴリズムも IEUを用いて導くことができ，前章の複素 NMF アルゴリズムにおいて式 (14) と式 (18) を入れ替えたものとなる。以上の関係より，(P0 ) を複素 NMF の双対形式と呼ぶことにする。
以上の複素 NMF の双対形式では，|Xl,k,m | とHk,l Ul,m の Euclid 距離が目的関数となっているが，|Xl,k,m | と Hk,l Ul,m はいずれも非負値のため，Euclid 距離の代わりに I ダイバージェンス規準を用いて両者の誤差を測ることもできる。そこで，

minimize JKL (X, H, U )X
|Xl,k,m |:=|Xl,k,m | log
Hk,l Ul,m(26)

l,k,m



− |Xl,k,m | + Hk,l Ul,m + R(U )X
subject to Yk,m =Xl,k,m(27)(28)

lのような最適化問題を考える。この問題を I ダイバージェンス規準の複素 NMF と定義する。ただし，R(U )は U に対するスパース正則化項であり，
XR(U ) = λ|Ul,m |p(29)l,mることで，2が成立する。明らかに等号は ξl,k,m = |Xl,k,m | のとき成立する。式 (32) の第二項には依然として微分不可能な点が含まれる。ここで，|Xl,k,m | の係数を Dl,k,m =ξl,k,mlog Hk,lUl,m −2 と置くと，Dl,k,m は ξl,k,m , Hk,l , Ul,mの値に応じて正にも負にもなりうる。そこで，l,k,mという不等式の形にまとめることができる。ただし，
Dl,k,m1
+(Dl,k,m ≥ 0)
 2bξl,k,ml,k,m(38)Al,k,m := 1(Dl,k,m < 0)
ξl,k,m
0(Dl,k,m ≥ 0)
(39)
Bl,k,m :=− Dl,k,m cl,k,m /2 (Dl,k,m < 0)である。
また，正則化項 R(U ) については，

R(U ) ≤ R+ (U , V ) X
λ(p|Vl,m |p−1 (Ul,m − Vl,m ) + |Vl,m |p )



l,m



(0 < p ≤ 1):=(40)Xλ
p−2 2p(p|V|U+(2−p)|V|)l,ml,ml,m2
l,m


(1 < p ≤ 2)

パラメータ推定アルゴリズムの導出

4.1 補助関数の設計
JKL (X, H, U ) は Xl,k,m に関して微分不可能な点を含んでおり，JKL (X, H, U ) を最小化する Xl,k,mは解析的に求まらない。しかし，下記に示す補助関数を用いることで JKL (X, H, U ) を停留点へ収束させる Hk,l , Ul,m , Xl,k,m を得るためのアルゴリズムを導くことができる。まず，対数関数は凹関数より，対数関数の接線が対数関数を下回らないことを利用す
ると，|Xl,k,m | log |Xl,k,m |− |Xl,k,m | log(Hk,l Ul,m )−|Xl,k,m |
(30)

1
(|Xl,k,m | − ξl,k,m )≤ |Xl,k,m | log ξl,k,m +ξl,k,m− |Xl,k,m | log(Hk,l Ul,m )−|Xl,k,m | (31)2
ξl,k,m|Xl,k,m |+ |Xl,k,m | log−2(32)=
ξl,k,mHk,l Ul,mが成り立つ。
以上より，+JKL(X, H, U , ξ, b, c, V )2∗:= Al,k,m |Xl,k,m | −2ℜ[Bl,k,mXl,k,m ]1+Hk,l Ul,m + R+ (U , V )+ξl,k,m|Xl,k,m |bl,k,m+2bl,k,m2

4.2 更新式
PXl,k,m には l Xl,k,m = Yk,m という制約があるため，ラグランジアン
L(X)
XX
2
∗
:=Al,k,m |Xl,k,m | − 2ℜ[Bl,k,mXl,k,m ]l2λk,m(33)2λ∗k,m(34)X++k,mX

Xl,k,m − Yk,m!

X∗∗
Xl,k,m− Yk,m!

lk,mが成り立つ（等号は bl,k,m = |Xl,k,m |, cl,k,m =Xl,k,m|Xl,k,m | ∈ C のとき成立する）ことを利用し，Dl,k,mの符号に合わせて上記のいずれかの不等式を適用す

日本音響学会講演論文集

k,mX

および

− |Xl,k,m | ≤ −ℜ[c∗l,k,m Xl,k,m ](41)

を JKL の補助関数として用いることができる。

2|Xl,k,m | ≤

(35)

(36)
12∗


|Xl,k,m | + Dl,k,m ℜ[cl,k,m Xl,k,m ] ξl,k,m
(Dl,k,m < 0)
12∗
= Al,k,m |Xl,k,m | − 2ℜ[Bl,k,mXl,k,m ] +(37)ξl,k,m≤とする。

4|Xl,k,m |+ Dl,k,m |Xl,k,m |ξl,k,m
Dl,k,m1
12

+|Xl,k,m | +

2bξ
ξ

l,k,ml,k,ml,k,m(D≥ 0)

l(42)+を用いて補助関数 JKLを最小にする Xl,k,m の更新式∗を求めることがきでる。L(X) の Xl,k,mについての偏微分を 0 と置くことで，∂L= 2Al,k,m Xl,k,m − 2Bl,k,m + 2λ∗k,m = 0∗
∂Xl,k,m(43)- 435 -2016年9月

∴ Xl,k,m =Bl,k,m − λ∗k,mAl,k,m(44)Pを得る。これを l Xl,k,m = Yk,m に代入することで，未定乗数 λ∗k,m をX

Xl,k,m =l

X Bl,k,m − λ∗k,ml∴ λ∗k,m =Al,k,mP Bl,k,m= Yk,ml Al,k,m − Yk,mP
1l Al,k,m(45)(46)

より得ることができる。式 (46) を式 (44) に代入する
ことで，Xl,k,m の更新式
P Bl,k,m Yk,m − l Al,k,m1 PBl,k,m +Xl,k,m =1Al,k,ml Al,k,m(47)+が得られる。次に，補助関数 JKLを最小にする Hk,l ,+Ul,m の更新式は，JKLをそれぞれの変数に関して偏微分したものを 0 と置くことで+X  |Xl,k,m |
∂JKL−
(48)
=+ Ul,m = 0∂Hk,lHk,lmP
m |Xl,k,m |∴ Hk,l = P(49)m Ul,m+∂JKL∂Ul,m
 X
|Xl,k,m |+ λp|Vl,m |p−1+H−k,l

Ul,m
k

= 0 (0 < p ≤ 1)

(50)
=X

|Xl,k,m |p−2
+ Hk,l + λp|Vl,m | Ul,m−Ul,mk= 0 (1 < p ≤ 2)

∴ Ul,mP
k |Xl,k,m |P(0 < p ≤ 1)
p−1
Hk,l + λp|Vl,m |kp P
P
( k Hk,l )2 +4λp|Vl,m |p−2 k |Xl,k,m |=(51)
2λp|Vl,m |p−2


P



k Hk,l(1 < p ≤ 2)
−2λp|Vl,m |p−2のように得られる。
以上をまとめると次のようになる。
Step 0) 初期値設定
Step 1) 補助変数の更新
ξl,k,m== |Xl,k,m |, bl,k,m=|Xl,k,m |,Xcl,k,m= |Xl,k,m, Vl,k,m = Ul,k,m ,l,k,m |Step 2) モデルパラメータの更新
式 (47),(49),(51).

5NMF(KL-NMF)，Euclid 距離規準 NMF(EU-NMF)，
板倉齋藤擬距離規準 NMF(IS-NMF) による各音源の平均 Signal to Noise Ratio (SNR) 改善値を比較した。
音響信号は全て 44.1kHz のモノラル信号で，フレーム長 32ms，フレームシフト 16ms で短時間 Fourier変換を行いスペクトログラムを算出した。実験は辞書 H を通常の NMF（反復回数:100）で作成し，5-分割交差検証（反復回数:30）を行った。また正則化パラメータ λ は EU-CNMF,KL-CNMF いずれも λ = 1とした。音源ごとの平均 SNR 改善値を Fig.1 に示す。
全ての音源で KL-CNMF，EU-CNMF が最も良い結果になっている訳ではないが，提案手法が複素 NMFに比べて高い性能になることが確認できる。

6

おわりにNMF 法に基づく音源分離手法では，振幅スペクトログラムの加法性という必ずしも成り立たない仮定のもので分離を行っているという問題があった。これを解決するため複素スペクトル領域での信号分解を可能にする複素 NMF という枠組みが提案されている。
NMF 法では経験的に I ダイバージェンスを規準とした場合に高い音源分離性能が得られることが知られているため，本稿では複素 NMF の双対形式に基づきI ダイバージェンス規準の複素 NMF アルゴリズムを導き，提案する I ダイバージェンス規準複素 NMF により従来の複素 NMF と比べて高い音源分離性能を得られることを確認した。
謝辞 本研究は JSPS 科研費 15K06081, 15H02757,26730100 の助成を受けて行われた。

参考文献

実験結果

提案手法の有効性を確認するため，SiSEC 2013
Database の”professionally produced music recordings”（5 音源からなるポピュラー音楽の音響信号）を用いて半教師あり音源分離実験を行い，I ダイバージェンス規準複素 NMF(KL-CNMF)，Euclid 距離基
準複素 NMF(EU-CNMF)，I ダイバージェンス規準

日本音響学会講演論文集Fig. 1 5 音源 (Source1-5) からなる音響信号分離後の平均 SNR 改善値。

[1] D.D. Lee and H.S. Seung, “Algorithms fornonnegative matrix factorization”, Proc. NIPS,
2000.
[2] H. Kameoka et al., “Complex NMF: A newsparse representation for acoustic signals”, Proc.
ICASSP 2009.
[3] M. Nakano et al.,“Convergence-guaranteed multiplicative algorithms for nonnegative matrix
factorization with β-divergence.” choice, 2010.
[4] D. Fitzgerald et al., “On the use of the Beta Divergence for Musical Source Separation”, Proc.
Irish Signals Syst. Conf., 2009.

- 436 -2016年9月