1-4-8
I ダイバージェンスに基づく独立低ランク行列分析の実験的評価 ∗☆最上伸一 (東大), 三井祥幹 (東大), 高宗典玄 (東大), 北村大地 (東大), 猿渡洋 (東大),高橋祐 (ヤマハ), 近藤多伸 (ヤマハ), 中嶋広明 (ヤマハ), 亀岡弘和 (NTT)1

はじめにブラインド音源分離（blind source separation:BSS）とは，複数の音源が混合された観測信号のみを用いて，混合前の個々の音源を推定する技術である．音源信号，観測信号及び分離信号をそれぞれ
sij = (sij1 , . . . , sijN )⊤ ∈ CN
(1)

xij = (xij1 , . . . , xijM )⊤ ∈ CM⊤

yij = (yij1 , . . . , yijN ) ∈ C

N(2)(3)とする．ここで，i = 1, . . . , I 及び j = 1, . . . , J はそれぞれ周波数ビン及び時間フレームのインデック
ス，n = 1, . . . , N 及び m = 1, . . . , M はそれぞれ音源，チャネルのインデックス，∗⊤ はベクトル及び行列の転置を表す．今，観測信号が
xij = Ai sij(4)と表現できると仮定する．ここで，Ai ∈ CM ×N は周波数ごとの混合行列である．さらに N = M かつ
Ai を正則と仮定すると，分離信号は
yij = Wi xij(5)と推定できる．ここで，Wi = (wi1 · · · wiN )H =H
A−1
i ，win は n 番目の音源の分離フィルタ，∗ はエルミート転置である．近年，音源間の統計的独立性と各音源の低ランク性を仮定した分離行列 Wi の推定法として，独立低ランク行列分析 (independent lowrank matrix analysis: ILRMA) [1] が提案されており，高精度な BSS を実現している．
一方，非負値行列因子分解 (nonnegative matrix factorization: NMF) を多チャネル信号に拡張した多チャネル NMF (multichannel NMF: MNMF) [2] は，全音源の低ランク時間周波数構造と音源毎の空間相関行列（混合系に対応する空間パラメータ）を同時に推定する BSS である．最適化の観点では，分離系の推定を行う ILRMA の方が MNMF よりも高速かつ高精度であることが報告されている [1]．従来の MNMF及び ILRMA は，板倉斎藤ダイバージェンスに基づ
く NMF に対応する原点対称複素正規分布生成モデルを仮定していたが，近年では原点対称複素 Student’s
t 分布に基づく MNMF (t-MNMF) [3] 及び ILRMA(t-ILRMA) [4] に一般化されている．さらに著者ら
は，複素 Poisson 分布 [5] を生成モデルに用いることで，一般化 Kullback–Leibler (KL) ダイバージェンス（I ダイバージェンス）規範で各音源を低ランク近似する ILRMA (KL-ILRMA) を提案している [6]．本
稿では区別のため，従来の ILRMA を IS-ILRMA と表記する．既存手法と KL-ILRMA の関係を Fig. 1に示す．
KL-ILRMA の音源モデル推定は一般化 KL ダイバージェンスの最小化に基づくため，その凸性より初期値等への頑健性が期待されるが，文献 [6] では十分な検証が行われていなかった．そこで本稿では，
KL-ILRMA の性能を評価するためにより詳細な実験を行う．まず，トイモデルを用いて，複素 Poisson 分布に従う音源分離に KL-ILRMA が適していること∗Fig. 1

Related methods of MNMF and ILRMA.

を示す．次に，音楽信号に対する音源分離実験を通して複数の従来手法と比較し，KL-ILRMA の有効性及び初期値に対する頑健性を示す．

2I ダイバージェンスに基づく ILRMA

2.1 生成モデル及びコスト関数
KL-ILRMA では，複素 Poisson 分布 [5] を音源モデルに仮定することで，一般化 KL ダイバージェンスに基づいて音源を低ランク近似する．複素 Poisson
分布の確率密度関数は，D = {z ∈ C | |z| ∈ N} を台とする複素数 z に対して次式で与えられる．
p−2 |z|pp |z|λe−λ(6)p(z) =p2π (|z| )!
NMF に基づく単一チャネル音源分離では，振幅スペクトログラムを用いた一般化 KL ダイバージェンス
規範 NMF (KL-NMF) が音源分離に適している [7]ことから，本稿でも p = 1 のみを考える．
KL-ILRMA では，各音源のスペクトログラムの複素時間周波数成分 yijn が，時間周波数ごとに異なるパラメータ λijn を持つ独立な複素 Poisson 分布に従って生成されると仮定する．従って，KL-ILRMA の音源モデルは以下の式で与えられる．
∏∏ |yijn |−1 λijn |yijn |p(yijn ) =e−λijn(7)2π|y|!
ijn
i,j,ni,j,nまた，λijn は低ランクな非負スペクトログラムとして，以下のようにモデル化される．
∑λijn =tikn vkjn
(8)k

ここで，tikn 及び vkjn は，それぞれ基底行列 Tn ∈RI×K
及びアクティベーション行列 Vn ∈ RK×Jの++成分である．R+ は非負値実数の集合を表す．k =1, . . . , K は基底インデックスであり，基底の本数 Kは通常 I, J と比較して非常に小さい値とする．式 (7)H
xij からの負対数尤度は，yijn = win[∑HH
L=log(|winxij |!) + log|win
xij |i,j,nExperimental evaluation of independent low-rank matrix analysis based on I divergence by Shinichi
Mogami, Yoshiki Mitsui, Norihiro Takamune, Daichi Kitamura, Hiroshi Saruwatari (The University ofTokyo), Yu Takahashi, Kazunobu Kondo, Hiroaki Nakajima (Yamaha), Hirokazu Kameoka (NTT).

日本音響学会講演論文集- 423 -2018年3月

H− |winxij | log

− 2J

∑

∑

tikn vkjn +∑k]tikn vkjnklog |det Wi | + const.

(9)iHHとなる．log(|winxij |!) は |winxij | に関する離散的な関数であるが，ここに Stirling の近似log(x!) ≈ x log x − x(10)を導入することにより，負対数尤度の式 (9) を初等関H数で近似しつつ win
xij の変域を複素数平面全域に拡張した以下のコスト関数を得る．
[∑HHJ =|winxij | log|winxij |i,j,nH
H− |winxij | + log|winxij |H− |winxij | log

− 2J

∑∑

tikn vkjn +k∑]tikn vkjn

klog |det Wi | + const.

(11)i

これは，一般化 KL ダイバージェンス DKL (y | x) =y log y − y log x − y + x を用いて∑[∑
HJ =DKL (|winxij | | k tikn vkjn )i,j,n]
H+ log|winxij |∑− 2J
log |det Wi | + const.

(12)iとも書けることから，音源間の独立性最大化と一般化 KL ダイバージェンスに基づく低ランク音源モデル推定を同時に達成するコストであると解釈できる．

2.2 音源モデルの更新式
式 (12) で与えられるコスト関数 J を Tn , Vn につH
いてみれば，これらが関与するのは DKL (|winxij | |∑tv)の項のみである．従って，
T及び
Vn
ikn
kjn
nk
は KL-NMF [8] と同様に以下の乗算更新式を反復することで最適化できる．
∑ |wH xij |
∑ in
vkjn
′′k′ tik n vk jnj∑(13)
tikn ← tikn ·j vkjn
∑ |wH xij |∑ intikn′′k′ tik n vk jni∑vkjn ← vkjn ·(14)i tikn
上記の一般化 KL ダイバージェンス最小化問題は各変数に関して凸であることが知られており，それが従来の MNMF や ILRMA とは異なる利点であると期待される．
2.3 分離行列の更新式
従来手法の IS-ILRMA は，分離行列の更新に iterative projection (IP) [9] と呼ばれる手法を用いた．
しかしながら，IP はコスト関数が − log |det Wi | とwin の二次形式の和で構成される場合にのみ適用でき，KL-ILRMA のコスト関数 (11) の形には適用できない．そこで，補助関数法を用いて分離行列 Wi のH反復更新式を導出する．まず，式 (11) の log|winxij |日本音響学会講演論文集に関して次の接線不等式
1HHlog|winxij | ≤(|winxij | − αijn ) + log αijnαijn(15)を適用すると，次の補助関数 J1 が得られる．
]∑[ 1HHJ ≤ J1 =|win
xij |2 + dijn |winxij |αijni,j,n∑− 2Jlog |det Wi | + const.
(16)
idijn =1αijn+ log αijn − log∑tikn vkjn − 2k(17)ここで αijn > 0 は補助変数であり，式 (15) は αijn =H|winxij | のときに等号成立する．
続いて，win についての微分が可能となるようにH
dijn |winxij | の上界関数を設計することを考える．
Fig. 2 が示す通り，d の符号で場合分けを行うこと
によって，任意の y0 に対して d |y| に点 y0（Fig. 2 では黒の点線）で一致するような上界関数を，y について 2 次（放物面）または 1 次（平面）の形で構成できHる．よって，dijn |winxij | に対する上界関数の式は

dijn1H|winxij |2 + dijn βijn


2β2
ijn


(dijn ≥ 0)H
dijn |win xij | ≤()1
HHx
dijn ωijn win
xij + ωijn win
ij
2

(dijn < 0)(18)と表される．ここで，∗ は ∗ の複素共役を表す．βijn >H0 は実数の補助変数であり，等号は βijn = |winxij |にて成立する．また，ωijn は |ωijn | = 1 を満たす複素H
H
数の補助変数であり，等号は ωijn = win
xij /|winxij |のときに成立する．式 (18) を式 (16) に適用すると，次の補助関数 J2 が得られる．
∑[]HHHJ1 ≤ J2 = Jwin
Uin win + win
rin + rin
win
i,n∑

− 2J1∑Uin =J jrin =log |det Wi | + const.

i(1max(0, dijn )+αijn
2βijn

1 ∑ ωijn
Jj2(19))xij xHij(20)

min(dijn , 0)xij(21)式 (19) は win の一次項を含むため，依然としてIP を適用することはできない．このような一次項を含む関数の最小化には，著者らの一部が開発した
vectorwise coordinate descent (VCD) [10] が適用できる．VCD では，まず変数行列 Wi の特定の行ベク
トル win に着目し，余因子展開によって式 (19) を以下の形式に変形する．
∑[]
HHH
J2 = Jwin
Uin win + win
rin + rinwini,n−J∑

2log|bHin win | + const.

(22)i

ここで，bin は Wi の余因子行列 Bi = (bi1 , . . . , biN )の列ベクトルであり，en を n 番目の要素のみが 1 の- 424 -2018年3月

Table 1 Generating conditions of toy models
Size of spectrogram
I = 257, J = 514
Number of bases
K = 10Scale parameter c
1, 2, 3, 4, 5, 10, 15, 20
Kurtosis of gamma
distribution(b)(a)

Source 2

Fig. 2 Shapes of auxiliary function (blue) of original function d |y| (green) with contact point y0(black dotted line). When d > 0, d |y| can be majorized by paraboloid of revolution (left). Whend < 0, d |y| can be majorized by plane (right).

を解くと，停留点は次式で与えられる [10]．
[]
1

−1
Ub−r
(if ubr = 0)√ininin)][ ubb (√4ubbubr
win =−1bin − rinU1− 1+2
 in 2ubb|ubr |(otherwise)(24)−1−1H
ここで，ubb = bH
in Uin bin , ubr = bin Uin rin である．

式 (24) に対して各補助変数の等号成立条件を代入することで，win の更新式が次のように得られる．
1H
dijn ←+ log|winxij |Hx ||winij∑− log
tikn vkjn − 2(25)1∑Uin ←J j(k1max(0, dijn )+Hx |Hx ||win
2|winijij)
xij xHij(26)

r̃in ←Hx
1 −1 ∑ 1 winij
Uinmin(dijn , 0)xij (27)
Hx |J2|w
in ijjw̃in ← (Wi Uin )−1 enHuww ← w̃inUin w̃inHuwr ← w̃in Uin r̃in(28)
(29)

(30)
w̃in


− r̃in(if uwr det Wi = 0)√ uww()√uwr4uww
win ←1− 1+w̃in − r̃in
2|uwr |
 2uww(otherwise)(31)

KL-ILRMA では，音源モデル Tn 及び Vn の更新
式 (13), (14) と，分離行列 Wi の更新式 (25)–(31) を交互に繰り返すことによって，式 (11) で表されるコスト関数 J の最小化を達成する．これらの更新式は全て補助関数法で導出されていることから，更新における J の単調非増加性が保証されている．

日本音響学会講演論文集

2 m Source 2Source 14040

E2A impulse
response
T60 = 300 ms

5.66 cm

2mSource 14020E2A impulse
response
T60 = 300 ms

2.83 cm

Fig. 3 Recording conditions of impulse responses
E2A obtained from RWCP database [12]: (a) IR1and (b) IR2.

SDR improvement [dB]

単位ベクトルとして bin = (det Wi )Wi−1 en と書ける．bin , Uin 及び rin は win と独立であるため，式(22) は他の行ベクトルを固定することで win のみの関数となり，win をブロック単位とした座標降下法が適用できる．式 (22) の win に関する停留点条件
()
∂J2
bin
= − H+ Uin win + rin = 0 (23)H∂winwin binT̃n : 12, Ṽn : 2302826242220181614IS-ILRMAKL-ILRMA

051015
Scale parameter c20Fig. 4 Average SDR improvements of IS-ILRMAand KL-ILRMA for toy model separation.

3

性能比較実験

3.1 トイモデルを用いた性能評価
本稿では，初めにトイモデルを用いて，KL-ILRMAが複素 Poisson 分布に従う音源の分離に適していることを確かめる．まず，独立なガンマ分布に従って行I×K
列 T̃n ∈ R+及び Ṽn ∈ RK×J(n = 1, 2) の各成分+をランダムに生成し，T̃n Ṽn の最大値を 1 に正規化することでモデル行列 Rn を得る．次に，Rn をスケール倍した行列 cRn の各 i, j 成分をパラメータ λ とする複素 Poisson 分布から，分離信号のスペクトログラ
ム Yn = (yijn )i,j を生成する．以上の手順により，複素 Poisson 分布に従う音源のトイモデルを得た．ト
イモデルの生成条件の詳細を Table 1 に示す．複素Poisson 分布 (6) の形状はパラメータ λ のスケールによって変化し，λ が大きいほど振幅 |z| が平均値 λ 付近の値を取りやすくなる．したがって，スケールパラメータ c の値が小さいとき Yn の振幅値はスパースになり，逆に大きいときは Yn の振幅値は cRn に近づき低ランク性の強い行列となる．
以上より得られる各音源のトイモデルをそれぞれ時間信号に変換し，Fig. 3 の (a) に示すインパルス応答 (IR1) を畳み込むことで，本実験で用いる混合音を作成した．これを従来手法の IS-ILRMA と提案
手法の KL-ILRMA で分離し，その性能を signal-todistortion ratio (SDR) [11] 改善量によって評価した．
なお，各 ILRMA の音源モデルの基底数 K はトイモデル生成時のパラメータ行列 T̃n Ṽn のランクと同一の 10 本とし，反復回数は 200 回とした．また，分離行列の初期値には理想解に対して 5 % の複素正規乱数を加えたものを用い，音源モデル Tn 及び Vn の初期値には一様乱数を用いた．

- 425 -

2018年3月

Table 3 Experimental conditions
Sampling frequency

16 000 Hz

STFT

512-ms-long Hamming
window with 128-ms shift
200
1020

Number of iterations
Number of bases
Number of trials

3.2 音楽信号の分離実験による既存手法との比較
次に音楽信号分離について，KL-ILRMA と従来手法である MNMF, t-MNMF, IS-ILRMA 及
び t-ILRMA との性能比較を行う．音楽信号として ，SiSEC2011 [13] で 公 開 さ れ て い る 音 楽 信 号
berlin-roads のうち A. guitar/vocal 及び A. guitar/piano の 2 種類の組合せを用いた．さらに各々に対して Fig. 3 の (a), (b) に示す 2 種類のインパルス応答 （IR1 及び IR2） を畳み込むことで，計 4 種類の混合音 No. 1–No. 4 を作成した．各混合音における実験条件の組合せを Table 2 に示す．IS-ILRMA,t-ILRMA 及び KL-ILRMA では，分離行列の初期値を単位行列，音源モデル Tn 及び Vn の初期値を一様乱数とした．MNMF 及び t-MNMF では，文献 [2] と同様の初期値を用いた．t-MNMF の自由度パラメータは ν = 1，t-ILRMA の自由度パラメータ及びドメインパラメータはそれぞれ ν = 3, p = 2 とし，t-ILRMA
では文献 [4] に基づいて ν に関する tempering 処理を行った．その他の実験条件を Table 3 に示す．
各音源における手法ごとの SDR 改善量の平均をFig. 5 に示す．エラーバーは標準偏差を表す．No. 2–No. 4 の 3 種類の混合音に対して，KL-ILRMA は他手法と比較して平均的に優れた性能を示している．また，いずれの混合音においても KL-ILRMA の標準偏差は他手法と比較して小さい．これは，凸なダイバージェンスを用いて音源モデルを推定する KL-ILRMAが音源初期値対して頑健であることを示している．

おわりに

本稿では，凸な一般化 KL ダイバージェンスを用いて音源モデルを低ランク近似する KL-ILRMA の性能を詳細な実験によって比較した．トイモデル実

日本音響学会講演論文集

MNMF
t-MNMF
IS-ILRMA
t-ILRMA
KL-ILRMA

1086420−2−4 ∼ ∼∼ ∼−12

スケールパラメータ c に対する各手法の性能をFig. 4 に示す．ただし，Tn 及び Vn の初期値の乱数
シードを変えて 20 回試行した際の平均 SDR 改善量を示し，エラーバーは標準偏差を表す．KL-ILRMA の性能は IS-ILRMA の性能を任意の c で上回っており，KL-ILRMA が複素 Poisson 分布に従う信号の分離に有効であることが示されている．また，KL-ILRMAの標準偏差は最大で 0.03 dB 程度と IS-ILRMA に比べて小さく，KL-ILRMA はより頑健な音源モデルの推定が可能であることが示唆されている．なお，提案 KL-ILRMA においては，式 (10) の近似を導入しているためスケールの大きな所では生成分布との不整合が生じているが，20 dB 程度の SDR は確保されていることより実用上は妥当な近似であることも示されている．

412SDR improvement [dB]

Table 2 Music sources obtained from SiSEC2011
Index
Source (1/2)
Impulse responce
No. 1 A. guitar/vocal
IR1No. 2 A. guitar/vocal
IR2
No. 3 A. guitar/piano
IR1
No. 4 A. guitar/piano
IR2−14No. 1∼∼No. 2No. 3Source

No. 4Fig. 5 Average SDR improvements for MNMF, tMNMF, IS-ILRMA, t-ILRMA, and KL-ILRMA.
験では，複素 Poisson 分布に従う音源に KL-ILRMAが適していることが示された．また，音楽信号の分離実験では，他の分離手法と比較して KL-ILRMA が主に頑健性の面で優れていることを示した．
謝辞 本研究は，セコム科学技術振興財団，JSPS 科研費
17H06572 の助成を受けた．

参考文献[1] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, andH. Saruwatari, “Determined blind source separation unifying independent vector analysis and nonnegative matrix factorization,” IEEE/ACM Trans. ASLP, vol. 24, no. 9, pp.
1626–1641, 2016.
[2] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “Multichannel extensions of non-negative matrix factorization withcomplex-valued data,” IEEE Trans. ASLP, vol. 21, no. 5,pp. 971–982, 2013.
[3] K. Kitamura, Y. Bando, K. Itoyama, and K. Yoshii, “Student’s t multichannel nonnegative matrix factorization forblind source separation,” in Proc. IWAENC, 2016.
[4] S. Mogami, D. Kitamura, Y. Mitsui, N. Takamune,
H. Saruwatari, and N. Ono, “Independent low-rank matrix
analysis based on complex Student’s t-distribution for blindaudio source separation,” in Proc. MLSP, 2017.
[5] H. Kameoka, “Towards a statistical audio signal processing
framework based on the I-divergence,” in Proc. 2011 Spring
Meeting of Acoustical Society of Japan, 2011, pp. 813–814.
[6] 最上伸一, 三井祥幹, 高宗典玄, 北村大地, 猿渡洋, 高橋祐, and 近藤多伸, “I ダイバージェンスを用いた独立低ランク行列分析,” in 第 32
回信号処理シンポジウム, 2017, pp. 354–359.
[7] D. Kitamura, H. Saruwatari, K. Yagi, K. Shikano, Y. Takahashi, and K. Kondo, “Music signal separation based onsupervised nonnegative matrix factorization with orthogonality and maximum-divergence penalties,” IEICE Trans.
Fundamentals, vol. E97-A, no. 5, pp. 1113–1118, 2014.
[8] D. D. Lee and H. S. Seung, “Algorithms for non-negative
matrix factorization,” in Proc. NIPS, 2000, pp. 556–562.
[9] N. Ono, “Stable and fast update rules for independent vector
analysis based on auxiliary function technique,” in Proc.
WASPAA, 2011, pp. 189–192.
[10] 三井祥幹, 高宗典玄, 北村大地, 猿渡洋, 高橋祐, and 近藤多伸, “空間事前情報を用いた独立低ランク行列分析,” in 第 32 回信号処理シンポジウム, 2017, pp. 360–365.
[11] E. Vincent, R. Gribonval, and C. Févotte, “Performancemeasurement in blind audio source separation,” IEEETrans. ASLP, vol. 14, no. 4, pp. 1462–1469, 2006.
[12] S. Nakamura, K. Hiyane, F. Asano, T. Nishiura, and T. Yamada, “Acoustical sound database in real environments forsound scene understanding and hands-free speech recognition,” in Proc. LREC, 2000, pp. 965–968.
[13] S. Araki, F. Nesta, E. Vincent, Z. Koldovský, G. Nolte,A. Ziehe, and A. Benichoux, “The 2011 signal separation evaluation campaign (SiSEC2011):-audio source separation,” in Proc. LVA/ICA, 2012, pp. 414–422.

- 426 -2018年3月