Vol.2013-MUS-100 No.8
2013/8/31

情報処理学会研究報告
IPSJ SIG Technical Report

モノラル音響信号に対する音源分離のための
無限半正定値テンソル分解
吉井 和佳1,a)

富岡 亮太2,b)

持橋 大地3,c)

後藤 真孝1,d)

概要：本稿では，モノラルの音響信号に対して音源分離を行うための新しい因子分解法について述べる．
この種の音源分離タスクにおいては，非負値行列分解 (NMF) を利用するのが現在の主流である．NMF では，与えられた混合音の振幅あるいはパワースペクトログラムを非負値行列とみなし，それを少数の基底ス
ペクトルとそれらの時間方向のアクティベーションとの積に分解することができる．しかし，NMF の性質
上，入力は非負値行列でなくてはならず，混合音の複素スペクトログラムに本来備わっている位相情報は考慮されていなかった．そのため，音源信号を再合成するには，音源信号のスペクトログラムは混合音のスペクトログラムと同一の位相をもつという強い仮定をおく必要があった．本研究では，周波数領域での位相復元を回避するため，混合音を時間領域で直接分解することができる半正定値テンソル分解 (PSDTF)を提案する．PSDTF は NMF の本質的な拡張となっているため，NMF と同様に効率的な乗法更新アルゴ
リズムに基づく最尤推定が可能であり，ガンマ過程に基づくノンパラメトリックベイズ拡張（基底数の無限化）も可能である．実験の結果，PSDTF は NMF より高品質な音源分離ができることを確かめた．

1. はじめに音楽音響信号に対する音源分離は，音楽情報検索 (MIR)

トログラムの位相は同一であると仮定することが一般的であるため，非適切な位相をもつスペクトログラムから高品
質な音源信号を再合成することはできなかった．

を支える基礎技術のひとつである．高品質な音源分離が実

これまで，実際の時間領域信号に対応する「無矛盾な」複

現できれば，歌唱者の声質や性別，楽器構成などの楽曲の素スペクトログラムを推定するため，多数の研究が行われ内容に基づいて，より緻密にユーザの好みに合う楽曲を検てきた．Griﬃn ら [6] は，与えられた振幅スペクトログラ

索できるようになる [1]．また，能動的音楽鑑賞 [2] の一

ムの位相を復元するため，その振幅スペクトログラムにで形態として，混合音中に含まれる既存の楽器パートを自分きるだけ近い振幅スペクトログラムをもつ時間領域信号を好みに編集しながら音楽鑑賞を楽しむことができるシステ

推定できる反復 STFT 法を提案している．Le Roux ら [7]ム [1–4] を実現するうえでも不可欠な技術である．

は，与えられた複素スペクトログラムの矛盾性を評価する

モノラル音響信号に対する音源分離タスクでは，非負ためのコスト関数を導出し，それを最小化するための効率

値行列分解 (Nonnegative Matrix Factorization: NMF) [5]

的なアルゴリズムを提案している [8]．一方，亀岡ら [9] は，を利用することが主流になっている．NMF は，入力とな混合音の複素スペクトログラムを複素成分の加法性に基づる非負値行列（振幅あるいはパワースペクトログラム）をいて直接分解することができる複素 NMF を提案している．

二つの非負値行列（基底スペクトルの集合と対応するアク

しかし，得られる音源信号の複素スペクトログラムは必ず

ティベーションの集合）の積で近似することができる．そしも無矛盾であるとは限らないため，Le Roux らのコストの後，ウィナーフィルタを用いて，混合音の複素スペクト関数を組み込む方法が提案されている [10]．ここで，複素ログラムを音源信号の複素スペクトログラムの和に分解す

スペクトログラムが無矛盾性を満たすことは、得られる音る処理が行われる．このとき，混合音と音源信号のスペク

源信号が高品質であることを必ずしも意味しないことに注意する必要がある．このことは，周波数領域における位相

1
2
3
a)
b)
c)
d)

産業技術総合研究所
1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan
東京大学
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
統計数理研究所
10-3 Midori-cho, Tachikawa, Tokyo 190-8562, Japank.yoshii(at)aist.go.jp
tomioka(at)mist.i.u-tokyo.ac.jp
daichi(at)ism.ac.jp
m.goto(at)aist.go.jp

c 2013 Information Processing Society of Japan

復元は容易ではないことを示唆している．
本研究では，位相復元を根本的に不要とするため，半正定値テンソル分解 (Positive Semideﬁnite Tensor Factoriza-

tion: PSDTF) [11] と呼ぶ新しい因子分解法を提案する．
PSDTF では，モノラルの混合音を時間領域において音源
信号に直接分解することができる．一般的に，時間領域に

1

Vol.2013-MUS-100 No.8
2013/8/31

情報処理学会研究報告
IPSJ SIG Technical Report
PSDTF

行列

基底

分解

合に

結
の凸

2.1 非負値行列分解: 離散フーリエ変換行列: フレーム における局所信号: フレーム における複素スペクトルNMF の目標は，非負値行列 X = [x1 , · · · , xN ] ∈ RM ×N
を低ランク近似する，すなわち二つの非負値行列 W =NMFは半正定値行列の対角要素のみに着目(本来存在するはずの周波数成分間の相関は無視)[w1 , · · · , wK ] ∈ RM ×K お よ び H = [h1 , · · · , hK ]T ∈解

分

RK×N の積で近似することである．

に

合

結

凸

の

ル

[dB]

ト

ク

ベ

NMF

X ≈ WH = Y
def

(1)

底

テンソルデータ: 局所的な共分散行列の集合 行列データ: パワースペクトルの集合 基(半正定値行列の集合)(非負値ベクトルの集合)

図1

周波数領域における音源分離：PSDTF は NMF の自然な拡張

ここで，wk および hk はそれぞれ基底ベクトルおよび対応するアクティベーションベクトルである．ただし，

K  min(M, N ) は基底数，Y = [y 1 , · · · , y N ] ∈ RM ×N
おいて，音響信号の位相や波形を何らかの関数形を用いては再構成行列を表す．行列積表現である式 (1) は，各 n に明示的に表現することは容易ではないため，我々は音響信関するベクトル和として書きくだすことができる．

号の統計的な性質に着目する．まず，少数の定常な基底信号の存在を考え，それぞれが異なるカーネル（基底カーネ
ルと呼ぶ）をもつ定常なガウス過程に従うと仮定する．与

xn ≈

K


hkn wk = y n
def

k=1

(2)


えられた混合音は，局所的にみると複数の基底信号の線形

ここで，y kn = hkn wk と定義すると，y n =和で構成されていると仮定する．ただし，結合係数は各時立する．観測ベクトル xn と再構成ベクトル y n との間の刻ごとに異なっている．このとき，混合音は局所的にガウ

誤差 C(xn |y n ) を評価する尺度として，Bregman ダイバース過程に従い，そのカーネルは基底カーネルの凸結合で与ジェンス [14] が広く利用されている．

えられる．我々の目標は，観測データである混合音の局所的な共分散から，それらを構成する少数の基底カーネルを求めることである．具体的には，乗法更新アルゴリズムに基づく最尤推定あるいは変分ベイズ法に基づくベイズ推定k y kn が成Cφ (xn |y n ) = φ(xn ) − φ(y n ) − φ (y n )T (xn − y n ) (3)
ここで，φ は厳密に凸な関数である．Bregman ダイバージェ
ンスは常に非負であり，xn = y n であるときに限りゼロをm (xm log xm − xm )が可能である．さらに，基底数が不明な場合でも，ガンマとる．その特別な形として，φ(x) =過程を用いて理論上無限個の基底の存在を仮定するノンパの 場 合 の Kullback-Leibler (KL) ダ イ バ ー ジ ェ ン ス や ，φ(x) = − m log xm の場合の Itakura-Saito (IS) ダイバー

ラメトリックベイズモデルを構成できる．
たが，実は周波数領域において等価な PSDTF を構成するジェンスがよく知られている．コスト関数 Cφ (X|Y ) =n Cφ (xn |y n ) を最小化する W および H を求めるため，ことができ，PSDTF は NMF の自然な拡張となっているこ

乗法更新アルゴリズム [15] が提案されている．これは，あこれまで時間領域における PSDTF について議論してきとが明らかとなる．図 1 に示すように，PSDTF では，各る確率モデルの最尤推定を行うことに相当する．一方，W

時刻における複素スペクトルの直積，すなわち半正定値行および H に適切な事前分布（通常はガンマ分布）を導入

列を少数の半正定値基底行列の和に分解する．一方，NMF

してベイズ推定を行う方法も提案されている [16, 17]．

では，上記直積行列の対角成分（パワースペクトル），すな
わち非負値ベクトルを少数の非負値基底ベクトルの和に分

2.2 音源分離への応用

解する．ここでの重要な違いは，PSDTF では異なる周波

我々の目的は，与えられた混合音を K 個の音源信号の和数ビン間の相関が考慮されていることである．このことのに分解することである．NMF を用いる場合，音源分離は

妥当性は，スペクトログラムを得るための標準的な方法で周波数領域で行うことになる．いま，混合音の「非負値」ス

ある短時間フーリエ変換やウェーブレット変換では，周波ペクトログラムを X ∈ RM ×N とし，M および N は周波数ビン間を完全に無相関化することはできないという事実

数ビン数およびフレーム数を表すものとする．このとき，に基づいている．したがって，調波構造などの周波数ビン

NMF を用いて X ≈ W H という分解を行うと，w k ∈ RM

間の強い相関を考慮しながら音源分離を行うことで，高品および hk ∈ RN はそれぞれ基底スペクトルおよび時間方

質な音源信号を復元することができる．

向の音量変化（アクティベーション）を表す．

2. 周波数領域における音源分離

潜在変数である音源信号の推定を行うには，NMF の確率モデルを考える必要がある．いま，混合音の複素スペクトロ

本章では，周波数領域における音源分離を行ううえで，非グラムを S = [s1 , · · · , sN ] ∈ CM ×N とし，k 番目の音源信

負値行列分解 (NMF) の理論的な裏付けについて説明する．

号の複素スペクトログラムを S k = [sk1 , · · · , skn ] ∈ CM ×N

特に，NMF のよく知られた 2 つの変種である KL-NMF [12]とする．混合音が K 個の音源信号の瞬時混合であると仮および IS-NMF [13] について議論する．

定すると，以下が成り立つ．

c 2013 Information Processing Society of Japan

2

Vol.2013-MUS-100 No.8
2013/8/31

情報処理学会研究報告
IPSJ SIG Technical Report

Sn =

K


Sk

i.e.,

k=1

sn =

K


skn

2.4 IS-NMF に基づく音源分離
(4)

k=1

観測変数として S が与えられると，潜在変数である S k の
期待値は次式で計算することができる．

E[sknm |snm ] =yknm
wkm hknsnm = snm
ynm
k wkm hkn

IS-NMF は，混合音の「パワー」スペクトログラムを分解するために用いることが一般的である [13]．すなわち，

xnm = |snm |2 とする．コスト関数は IS ダイバージェンスに基づくものであり，次式で与えられる．

(5)

この処理はウィナーフィルタリングと呼ばれ，S k の位相
は S の位相と同一であるという仮定がおかれている．最後に，重畳加算合成法 [18] を用いれば，E[S k |S] から k 番目の音源信号を復元することができる．しかし，得られた音
源信号から複素スペクトログラムを再度計算してみても，

E[S k |S] とは等しくはならない．なぜなら，E[S k |S] は実際の時間領域信号に対応する「無矛盾な」スペクトログラ
ムではないからである．これは NMF に基づく音源分離の原理的な欠陥であり，音源信号の合成品質に限界がある根
本的な原因のひとつとなっている．

CIS (xn |y n ) =

M 

xnm
m=1

ynm

− log

xnm
−1
ynm



式 (10) はスケール不変であるので，理論的には IS-NMFの方が音源分離に適していることが知られている．
確率モデル化のためには，潜在変数 sknm が yknm を分散パラメータとする複素ガウス分布に従うことを仮定する．

sknm |yknm ∼ Nc (0, yknm )ここで，snm =k sknm かつ ynm =(11)k yknm であることを思い出すと，複素ガウス分布の再生性から次式を得る．

snm |ynm ∼ Nc (0, ynm )2.3 KL-NMF に基づく音源分離
KL-NMF は，混合音の「振幅」スペクトログラムを分解するために用いることが一般的である [12]．すなわち，

xnm = |snm | とする．コスト関数は KL ダイバージェンスに基づくものであり，次式で与えられる．

M 

xnm
x
CKL (xn |y n ) =− xnm + ynm (6)nm logynm
m=1

ただし，式 (6) はスケール不変ではないことに注意する．すなわち，任意の実数 α に対して CKL (X|Y ) = CKL (αX|αY )とはならない．このことは，与えられた混合音全体の音量を変化させると音源分離結果も変化することを意味しており，理論的には妥当であるとはいえない．
確率モデル化のためには，潜在変数 |sknm | が yknm を平均パラメータとするポアソン分布に従うことを仮定する．
|sknm |yknm ∼ Poisson(yknm )(7)
ここで，振幅の加法性が成立することを仮定すると，式 (4)
は |snm | = k |sknm | と書きなおせる．すなわち，混合音のスペクトログラムの位相は音源信号のスペクトログラムの位相と同一であると仮定することを意味する．これは通常成立しないが，実用上は都合がよい．いま，xnm = |snm |かつ ynm = k yknm であることを思い出すと，ポアソン
分布の再生性から次式を得る．
xnm ynm ∼ Poisson(ynm )(8)|snm | = k |sknm | を満たすポアソン変数 {|sknm |}K
k=1

(10)

(12)

したがって，xnm = |snm |2 は指数分布に従うことが分かる．

xnm |ynm ∼ Exponential(ynm )
snm =



(13)

K
k sknm を満たすガウス変数 {sknm }k=1 に基づく確率モデルを用いると，snm が与えられたときの sknmの期待値は式 (5) で，分散は次式で求めることができる．
−1
V[sknm |snm ] = yknm − yknm ynm
yknm

(14)

3. 時間領域における音源分離
本章では初心に戻り，音源分離を時間領域における分解問題として考えなおす．具体的には，音源信号の重畳に対する確率モデルを提案し，潜在変数である音源信号を確率的な枠組みで推定する方法について議論する．

3.1 問題設定
我々の目標は，与えられた混合音を K 個の音源信号の和に分解することである．いま，観測データとして N 個の実数値ベクトルの集合 O = [o1 , · · · , oN ] ∈ RM ×N を考える．ここで，on ∈ RM は，長さ M の窓を用いて混合音から切り出された n 番目のフレームにおける局所的な信号である．このとき，on の分解は次式で与えられる．

on =

K


okn

(15)

k=1

に基づく確率モデルを用いると，xnm = |snm | が与えられ

ここで，okn は k 番目の音源信号から切り出された n 番目たときの |sknm | の期待値は次式で与えられる．
−1
E[|sknm ||sknm |] = yknm ynm|sknm |のフレームにおける局所信号である．離散フーリエ変換行

最終的に，位相の保存性から式 (5) を得る.

c 2013 Information Processing Society of Japan(9)

列を F ∈ CM ×M とすると，sn = F on かつ skn = F okn
であるから，式 (15) で与えられる時間領域表現は式 (4) で与えられる周波数領域表現と等価である．

3

Vol.2013-MUS-100 No.8
2013/8/31

情報処理学会研究報告
IPSJ SIG Technical Report

短時間窓を用いて切り出した局所的な信号
時間2
okn |π, V ∼ N (0, πkn
V k)

(17)

さらに，式 (15) に着目すると，ガウス分布の再生性から時間変化する振幅をもつ正弦波
図 2から okn もやはりガウス分布に従うことが分かる．

局所信号 okn および okn は異なる位相や振幅をもっているが，同一の周期に従っていることに着目する.

モノラルの混合音を分離することは不良設定問題である
ため，何らかの制約が必要になる．本研究では，k 番目のon もガウス分布に従うことが分かる．

 K

2
on |π, V ∼ N 0,
πkn V k(18)

k=1

式 (18) には，基底信号の具体的な波形を表す {φkn }K
k=1 が音源信号における局所的な波形 {xkn }N
n=1 は一見異なって含まれておらず，基底信号を周辺化することであらゆるいるものの，何らかの統計的な性質（周期性や白色性など）{φkn }K
k=1 の位相を明示的に考える必要がなく，より頑健を共有していると仮定する．例えば，図 2 に示すように，
音源信号が時間変化する振幅をもつ正弦波であれば，局所
信号 okn および okn (n = n ) は異なる位相や振幅をもっているものの，定常な周期で振動している．そこで，okn波形の可能性を考慮した表現となっている．したがって，
2≥ 0，
な音源分離ができると期待できる．いま，hkn = πkn
T*1
X n = on on  0，Y n =k hkn V k  0 とすると ，

式 (18) から X n の対数尤度は次式で与えられる*2 ．

1
1clog p(X n |Y n ) = − log |Y n | − tr(X n Y −1n )22に対して，非定常な因子 πkn と定常な因子 φkn とに基づく
分解 okn = πkn φkn を考える．ここで， φkn は k 番目の定常な「基底信号」から切り出された n 番目のフレームにおける局所信号であり，πkn はその係数である．これらは非負値に限定されていないことに注意する．
結局，観測データ O が与えられたときに，各 k について
局所信号 {okn }N
n=1 を求めることが目標となる．k 番目の音源信号は，重畳加算法合成法 [18] を用いて得ることができる．これらの処理は全て時間領域で実行されるため，短(19)

いま，観測データとしてテンソル X = [X 1 , · · · , X n ] ∈
n log p(X n |Y n )RM ×M ×N が与えられたとき，対数尤度が最大となるような H = [h1 , · · · , hK ]T ∈ RK×N およびV = [V 1 , · · · , V K ] ∈ RM ×M ×K を求めたい．この問題は，あとの 4 章に示すように，入力となる X n がランク 1 の半正定値行列 (X n = on oTn ) に限定された PSDTF の特別な場合である. したがって，4.3 節で述べる乗法更新アルゴリ
ズムを利用して最尤推定することができる．

時間フーリエ変換やその逆変換などは必要ない．

3.3 確率モデルに基づく音源分離

3.2 確率モデルの定式化
式 (15) で与えられる分解に対する確率モデルを定式化したい．定常な基底信号に対しては，局所信号 {φkn }N
n=1が同じ共分散をもつことが期待できるので，同一の多次元
ガウス分布を事前分布として考えることができる．

φkn |V ∼ N (0, V k ),H および V が求まれば，局所信号 okn = πkn φkn を確
率的な枠組みのもとで推定することができる．ここで，πkn
や φkn を先に求めておく必要はなく，直接 okn の事後分布を計算することができる．具体的には，式 (17) および式 (18) から，okn の事後分布もガウス分布となることが分(16)

かり，その平均と分散は次式で与えられる．

ここで，共分散行列 V k ∈ RM ×M は対角行列に限定されておらず，任意の半正定値行列でよい．音響信号は通常ゼロを中心とする実数として記録されるため，平均ベクトルは0 とした．式 (16) は，k 番目の基底信号がカーネル V k をもつ定常ガウス過程に従うことを意味する．本来，音響信
号は離散的な時刻上に存在する不連続な系列ではなく，連続時間上に存在するなめらかな関数とみなすべきであるので，連続関数に対する確率分布を考えることは本質的である．本研究では，任意の M 点における周辺分布が式 (16)で与えられるガウス分布となるため，基底信号が従う確率分布はガウス過程に他ならない．例えば，もし V k が周期
カーネルであれば，局所信号 {φkn }N
n=1 は異なる位相をもつ一方，同一の周期に従うことが期待できる．

(20)

V[okn |on ]=Y kn − Y kn Y −1
n Y kn(21)

ここで，Y kn = hkn V k  0 とした (Y n =k Y kn  0)．

式 (20) は，観測信号 on のうち，位相や波形を明示的に考えることなく，カーネル V k で表現できる成分のみを通過させるウィナーフィルタである．時間領域における式 (20)および式 (21) は，周波数領域における式 (5) および式 (14)と形式上は非常に類似している．ただし，我々の提案法ではベクトル okn ごとに独立な計算になっているのに対し，
従来はベクトル skn の要素 sknm ごとに独立な計算になっている点が異なる．したがって，提案法では，周波数ビン間の相関構造を考慮しながら音源分離を行っていることが示唆される（詳細は次節参照のこと）．

次に，観測信号 on の尤度について考える．okn と φkn の*1間には線形性 okn = πkn φkn が成立しているため，式 (16)*2c 2013 Information Processing Society of JapanE[okn |on ]=Y kn Y −1n on

Ψ が半正定値行列であるとき Ψ  0 と書く．
c
= は定数項を除いて等号が成立することを示す．

4

Vol.2013-MUS-100 No.8
2013/8/31

情報処理学会研究報告
IPSJ SIG Technical Report

3.4 フーリエトリックに基づく近似的高速計算
従来と同様に周波数領域での定式化について議論する（図 1）．離散フーリエ変換行列を F ∈ CM ×M とすると，
式 (18) から，観測信号 on の線形写像である複素スペクト
ル F on は複素ガウス分布に従うことが分かる．
 K


H
F on |H, V ∼ Nc 0,
hkn F V k F

(22)

k=1

ここで，複素スペクトルの要素間の完全な共分散構造が考慮されていることに注意する．
式 (22) は式 (11) で与えられる IS-NMF の確率モデルの

X≈

K


hk ⊗ V k = Y
def

(23)

k=1

ここで，⊗ はクロネッカー積を表す．式 (23) は，各 n に関する行列和として書きくだすことができる．

Xn ≈

K


hkn V k = Y n
def

(24)

k=1

ここで，hkn ≥ 0 は n 番目の要素 X n における k 番目の基
底行列 V k の重みである．観測行列 X n と再構成行列 Y nとの間の誤差 Cφ (X n |Y n ) を評価する尺度として，本研究では Bregman 行列ダイバージェンス [14] を利用する．

自然な拡張となっていることを示している．数学的には，

V k が循環行列であれば，F V k F H は対角行列となること

Cφ (X n |Y n )が知られている（V k は対角化されるという）．このとき，= φ(X n ) − φ(Y n ) − tr ∇φ(Y n )T (X n − Y n )式 (22) は，複素スペクトルの各要素が独立であると仮定している式 (11) と等価になる．自明な例は，V k が単位行列である場合，すなわち，φkn が定常な白色ガウス性雑音である場合である．一方，V k が周期カーネルであり，窓幅 M がその周期より十分に大きければ，V k は循環行列と類似した斜行型の要素配置をもつ．しかし，厳密には V k
は循環行列とは異なるため，複素スペクトルの要素間に相関が発生してしまう (図 4 および図 5 の実験結果を参照)．
ただし，音響信号は白色成分と周期成分とで構成されているとみなせるため，我々の提案法の近似的高速算法として

IS-NMF は有用かつ妥当であるといえる．

4. 半正定値テンソル分解
本章では，半正定値テンソル分解 (Positive Semideﬁnite(25)

ここで，φ は微分可能で厳密に凸な関数である．Bregman
行列ダイバージェンスは常に非負であり，X n = Y n であるときに限りゼロをとる．本稿では，その特別な形として，

φ(Z) = − log |Z| である場合の Log-Determinant (LD) ダイバージェンス [19] を用いる場合に着目する．

CLD (X n |Y n ) + tr X n Y −1 − M= − log X n Y −1
n
n

(26)

非 負 の 実 数 x, y 上 に 定 義 さ れ る IS ダ イ バ ー ジ ェ ン ス

CIS (x|y) = − log(x/y) + x/y − 1 は，M = 1 とした LD
ダイバージェンスの特別な場合である．
我々の目標は，コスト関数 CLD (X|Y ) =n CLD (X n |Y n )を最小化するような H = [h1 , · · · , hK ] ∈ RN ×K および

Tensor Factorization: PSDTF) とよぶ新しい因子分解法

V = [V 1 , · · · , V K ] ∈ RM ×M ×K を求めることである．たについて説明する．NMF は N 個の非負値ベクトル（行列だし，H に関しては非負値制約が，V に関しては半正定データ）を K 個の非負値ベクトルの凸結合で近似するのと値制約が課されていることに注意する．このような因子分

同様に，PSDTF は N 個の半正定値行列（テンソルデータ）

解法を LD-PSDTF と呼ぶことにする．

を K 個の半正定値行列の凸結合で近似する．したがって，

PSDTF は NMF のエレガントな拡張になっており，乗法
更新アルゴリズムやノンパラメトリックベイズモデルに基4.2 補助関数法
本研究では，解析的な計算を可能にするため，補助関

づく基底数の無限化などの様々な技術が応用可能である．

数法 [15] を用いて，コスト関数 CLD (X|Y ) を Y （H およ

4.1 問題設定

F (θ) を θ に関して最小化すべき関数であるとすると，はじめに，我々が取り組む問題について定義する．観
測データとして，3 階のテンソル X = [X 1 , · · · , X n ] ∈び V ）に関して間接的に最小化することを考える．いま，F (θ) ≤ F + (θ, φ)(27)

RM ×M ×N が与えられるものとする．ここで，テンソルを満たす F + (θ, φ) を F (θ) の補助関数と呼ぶ．ここで，φの各要素 X n ∈ RM ×M は実対称半正定値行列であるとす

は補助パラメータである．このとき，以下の反復更新則

る．X n ∈ CM ×M が複素エルミート半正定値行列である
場合でも同様に扱うことができるが，本稿では簡単のため

X n ∈ RM ×M の場合について議論する．

φnew ←argminφ F + (θ old , φ)(28)←argminθ F + (θ, φ

(29)

θ

new

new

)

PSDTF の目標は，それぞれの半正定値行列 X n を K 個を用いると，F (θ) は単調非増加であることが証明できる．

の半正定値行列 {V k }K
k=1（基底行列とよぶ）の凸結合で近このアルゴリズムの収束性は保証されており，IS-NMF の似することである．

ベイズ学習でも同様の手法が利用されている [17].

c 2013 Information Processing Society of Japan

5

Vol.2013-MUS-100 No.8
2013/8/31

情報処理学会研究報告
IPSJ SIG Technical Report

CLD (X|Y ) に対する補助関数 ULD (X|Y ) を得るには，半正定値行列を変数とする関数の凸性および凹性に基づく不等式が必要になる．まず，f (Z) = log |Z| が凹関数である
ことに着目すると，f (Z) に対して 1 次のテイラー展開を次に，式 (32) を V k について微分してゼロとおき，式 (33)を代入すると，次式を得る．
old
V k P k V k = V old
k Qk V k(35)

ここで，V oldk は V k の現在の値である．P k および Qk は行うことで，次式を得る．

log |Z| ≤ log |Ω| + tr(Ω−1 Z) − M(30)

ここで，Ω は任意の半正定値行列（展開点）であり，M は Zのサイズである．等号成立条件は，Ω = Z で与えられる．
次に，任意の半正定値行列 A に対して g(Z) = tr(Z −1 A)

半正定値行列であり，次式で与えられる．

Pk =

N


hkn Y −1
n

Qk =

n=1

N


−1
hkn Y −1
(36)n X nY nn=1

ここで，Qk は半正定値行列であるので，ある下三角行列は凸関数であることに着目すると，澤田らの提案する不等Lk についてコレスキー分解 Qk = Lk LTk が可能である．

式 [20] を適用可能である．

これを用いると，式 (35) を解析的に解くことができ，最終


tr

K

k=1 Z k

−1

 
K
T
A ≤
tr Z −1k Φk AΦk

的に以下の更新則を得る．

(31)

k=1

ここで，{Z k }K
{Φk }K
k=1 は任意の半正定値行列であり，k=1 はV k ← V k Lk (LTk V k P k V k Lk )− 2 LTk V k
1(37)

行列の半正定値性の定義に従えば，ある行列 A が実対称
足すと単位行列になるような補助行列である ( k Φk = I)．

等号成立条件は，Φk = Z k ( k Z k )−1 で与えられる．

半正定値行列であることは，A = ZZ T を満たす実行列 Z

式 (30) および式 (31) を用いると，式 (26) に対する補助て V k の半正定値性は自然に保たれていることが分かる．

関数 ULD (X n |Y n ) を導くことができる．

が存在することと同値である．したがって，式 (37) におい

4.4 IS-NMF との関連性c

CLD (X n |Y n ) = log |Y n | + tr X n Y −1n−M
≤ log |Ωn | + tr Y n Ω−1nT+ k tr Y −1
(32)
kn Φkn X n Φkn≤ log |Ωn | + k tr hkn V k Ω−1
−M
n

def
−1
T
+ k tr h−1
kn V k Φkn X n Φkn = ULD (X n |Y n )LD-PSDTF は IS-NMF の自然な拡張であるが，X n お
よび V k が対角行列であれば，LD-PSDTF は IS-NMF と等しくなる．このとき，任意の半正定値行列の対角成分は非負値ベクトルであるので，式 (26) のコスト関数はCLD (X n |Y n ) = CIS (diag(X n )|diag(Y n ))(38)

ここで，Ωn は半正定値行列であり，{Φkn }K
k=1 は足すと単とでき，式 (34) および式 (37) で与えられる乗法更新則

位行列となるような補助行列である．等号が成立する，すは IS-NMF のための収束保証付きの乗法更新則 [15] と一

なわち，ULD (X n |Y n ) を最小化するときの条件（補助パラ

致する．IS-NMF に対しては，Expectation-Maximization

メータの反復更新則）は次式で与えられる．

(EM) アルゴリズムや収束保証がない乗法更新アルゴリズ

Ωn = Y nΦkn = Y kn Y −1n(33)

ムも利用可能であることが知られている [15]．LD-PSDTF
でも同様であるが，経験的な収束性の良さから，本稿では

4.3 乗法更新アルゴリズム

補助関数 ULD (X|Y ) = n ULD (X n |Y n ) を単調減少させることができる乗法更新アルゴリズムを導く．ここでは，スケールの任意性を除くため，tr(V k ) = 1 という制約をおく（あとで説明する LD-PSDTF のベイズ学習においてはこのような制約は必要ない）．tr(V k ) = s である場合
は，CLD (X n |Y n ) および ULD (X n |Y n ) の値を変化させずに，V k ← 1s V k かつ hkn ← shkn と更新できる．
まず，式 (32) を hkn について微分してゼロとおき，式 (33)

収束保証付きの乗法更新アルゴリズムを提案した．

LD-PSDTF を音源分離に利用するには，X n = on oTn とすればよい．式 (19) で与えられる対数尤度の符号を反転させると，式 (26) で与えられるコスト関数と定数項を除いて等しい．したがって，LD-PSDTF を用いて H および V の最尤推定を行うことができる．ただし，本来の LD-PSDTF
はランク 1 に限らず，任意のランクの行列 X n に対して適用可能であることを明記しておく．

5. 無限半正定値テンソル分解

を代入すると，以下の更新則を得る．

hkn ← hkn

−1
tr Y −1n V kY n Xn

tr Y −1n Vk

本章では，理論上は可算無限個の基底行列を内包する無(34)

限半正定値テンソル分解 (iPSDTF) について述べる．これまで，適切な基底数 K は事前に指定するものであるとしてこれは，hkn に非負係数を乗ずる乗法更新則となっており，議論してきたが，実際には容易ではない場合が多い．しかhkn の非負性は自然に保たれている．

し，異なる K で PSDTF を何度も行い，何らかの基準のも

c 2013 Information Processing Society of Japan

6

Vol.2013-MUS-100 No.8
2013/8/31

情報処理学会研究報告
IPSJ SIG Technical Reportとで最適な K をあとから選ぶ古典的なモデル選択では，計

5.1.2 無限モデル

算コストが非常に高くなる問題があった．一方，ノンパラ

次に，基底数が K → ∞ である場合のベイズモデルを設

メトリックベイズモデルに基づく iPSDTF では，観測デー計する．重要なことは，無限個の基底行列が存在するとしタに合わせて，無限個存在する基底行列のうちのごく一部ても，観測データを表現するのに実質的には少数の基底行のみが実質的にアクティベートされる機構を備えている．

列しか利用されないという点である．これを実現するため，

IS-NMF に対するノンパラメトリックベイズモデル [17] を
5.1 確率モデルの定式化参考にする．いま，K → ∞ とすると，基底行列の重みを

LD-PSDTF のベイズモデルについて説明する．まず，
式 (24) を少し変形した以下の分解を考える．

Xn ≈

K


表す θ = [θ1 , · · · , θ∞ ]T は無限次元の非負値ベクトルとなり，そのうちの一部のみが有意に大きな値をもち，それ以外はほとんどゼロとなるようにしたい．したがって，θ にθk hkn V k = Y n
def

(39)

k=1

対する事前分布としてガンマ過程を用いるのが自然である．本研究ではまず，以下のガンマ分布を考える．

ここで，θk ≥ 0 は k 番目の基底行列の重みであり，θk = 1 とすれば通常の LD-PSDTF と一致する．θ = [θ1 , · · · , θK ]Tの導入は基底数を無限化する際に重要な役割を果たす．

5.1.1 有限モデル
まず，基底数 K が有限である場合のベイズモデルを設計する．具体的には，H および V を確率変数とみなし，なんらかの事前分布を導入する．θ については θk = 1 とし

θk ∼ G(αc/K, α)(43)

ここで，α ≥ 0 および c ≥ 0 は正の超パラメータであり，
Eprior [θk ] = c/K および Eprior [ k θk ] = c となっている．

K → ∞ とすれば，以下のガンマ過程が得られる．
G ∼ GaP(α, G0 )(44)ておく．非負値 hkn ≥ 0 および半正定値行列 V k  0 に対

ここで，G0 はある空間 Θ 上に定義された基底測度であり，する事前分布としては，ガンマ分布およびウィシャート分

G(Θ) = c を満たす（確率測度とは限らない）．サンプル

布を用いることが一般的である．

される G は Θ 上の離散測度となることが知られており，

hkn ∼ G(a0 , b0 )V k ∼ W(ν0 , V 0 )(40)

E[G] = G0 かつ V[G] = G0 /α となっている．空間 Θ の微小区間への分割を {Θ1 , · · · , Θ∞ } とすると，G(Θk ) = θk

ここで，a0 および b0 は，ガンマ分布の形状およびレートとなっている．α は集中度と呼ばれ，α が小さくなるほどパラメータであり，ν0 および V 0 はウィシャート分布の自θ はよりスパースになる (偏り V[G] が大きくなる)．

由度および尺度行列である．
与えられた半正定値行列 {X n }N
n=1 は，それぞれ独立に
ウィシャート分布に従うと仮定する．

 K



νX n θ, H, V ∼ W ν,θk hkn V k

は扱えないが，K を α に比べて十分大きな値に設定すれば，式 (43) はガンマ過程の良い近似となる．このとき，基

(41)

k=1

ここで，ν はウィシャート分布の自由度である．X n に νを乗じる理由は，観測行列 X n の期待値が再構成行列 Y nと等しくなるようにするためである (E[X n ] = Y n )．このとき，ν−1
M であれば，M[X n ] = ν−M
Y n ≈ Y n とな
νる．一方，ν < M であれば，M[X n ] は定義されず，X n
はランク落ちになる．また，M = ν = 1 であれば，式 (41)は，式 (13) で与えられる指数分布と等価になる（IS-NMFと等価）．X n の対数尤度は次式で与えられる．

ν −M −1
log p(X n |Y n )=C(ν) +log |X n |2νν− log |Y n | − tr X n Y −1n2
2

最終的に，GaP-LD-PSDTF のベイズモデルは，式 (40)，
式 (41) および式 (43) で与えられる．計算機上では K → ∞(42)

底測度 G0 は一様な測度であることを仮定している．

5.2 変分ベイズ法
我々の目標は，観測データ X が与えられたとき，ベイズの公式を用いて確率変数の事後分布 p(θ, H, V |X) =p(X, θ, H, V )/p(X) を計算することである．しかし，周p(X, θ, H, V )dθdHdV の計算は解析
辺尤度 p(X) =的に行えないため，変分ベイズ法を用いて p(θ, H, V |X)
を近似的に求めることにする．まず，次式のような因子分解が可能な関数形をもつ分布 q(θ, H, V ) を考える．

N


K


q(θk )q(θ, H, V ) =q(hkn ) q(V k )
(45)
k=1

n=1

ここで，C(ν) は ν のみを含む定数項であり，観測行列 X n

真の事後分布 p(θ, H, V |X) ではこのような変数間の独立のみに依存する第二項も定数項となる．式 (42) と式 (26)を比較すると，対数尤度 log p(X|Y ) = n log p(X n |Y n )性は成立しないが，q(θ, H, V ) と p(θ, H, V |X) の間の KL

を Y に関して最大化する問題は，コスト関数 CLD (X|Y ) =n CLD (X n |Y n ) を最小化する問題と等価である.
c 2013 Information Processing Society of Japan

ダイバージェンスを最小化するような q(θ, H, V ) を求めたい．これは，対数周辺尤度 log p(X) の変分下限 L を最大化することと同値であることが知られている．

7

Vol.2013-MUS-100 No.8
2013/8/31

情報処理学会研究報告
IPSJ SIG Technical Reportγ

log p(X) ≥ E[log p(X|θ, H, V )]−11(ρ/τ ) 2
xγ−1 e− 2 (ρx+τ x )GIG(x|γ, ρ, τ ) =√2Kγ ( ρτ )+ E[log p(θ)] + E[log p(H)] + E[log p(V )]− E[log q(θ)] − E[log q(H)] − E[log q(V )] ≡ L

(46)

(51)

ここで，γ および ρ, τ ≥ 0 はパラメータであり，Kγ は第二種変形ベッセル関数を表す．式 (51) の対数は log x，x お

ここで，第一項では式 (42) に対する期待値計算が必要であ

よび x−1 の項を含んでいる．ガンマ分布は GIG 分布の特るが，解析的な計算は依然困難である．本研究では，L ≥ L

別な場合であるので，一種の共役性が成立している．このとなるさらなる変分下限 L を設計し，L の最大化を通し性質は LD-PSDTF の特別な場合である IS-NMF のベイズて L を間接的に最大化することを考える．これは，−L を

推定 [17] においても利用されている．また，期待値 E[x]−L の補助関数とみなした場合の補助関数法となっている．
具体的には，式 (30) および式 (31) を用いると，第一項の変分下限は次式で与えられる．

および E[x−1 ] は次式で計算可能である．
√
  √
√
√
τ Kγ+1 ( ρτ )
ρKγ−1 ( ρτ )1
E[x] = √
E
= √
(52)
√
√
ρKγ ( ρτ )x
τ Kγ ( ρτ )一方，q(V k ) は行列 GIG (Matrix GIG: MGIG) 分布 [21]E[log p(X|θ, H, V )]となることが分かる．MGIG 分布は半正定値行列上に定義

N


ν
c
E[log |Y n |] + E tr X n Y −1
(47)
=−
n
2 n=1

N
K



ν
≥−
log |Ωn | +−ME tr θk hkn V k Ω−1
n
2 n=1
k=1

K



T
−1 −1 −1+
E tr θk hkn V k Φkn X n Φkn
k=1

された分布であり，次式で与えられる．
M +12γM|X|γ− 2
|T |γ Bγ (RT /4)


1
−1
(53)
exp − tr RX + T X
2

MGIG(X|γ, R, T ) =ここで，γ は実数であり，R, T  0 は半正定値行列である．

M は行列 X のサイズであり，Bγ は第二種行列ベッセル関

ここで，Ωn は任意の半正定値行列，{Φkn }K
k=1 は足すと単数 [22] を表す．式 (53) の対数は log |X|，X および X −1

位行列となるような補助行列である．等号が成立するときの項を含んでいる．ウィシャート分布や逆ウィシャート分の条件は次式で与えられる．

布は，MGIG 分布の特別な場合である [23]．ただし，現在

Ωn =のところ期待値 E[X] および E[X −1 ] を解析的に計算するk E[Y kn ]
 −1 
 −1  −1Φkn = E Y −1knk E Y k n

−1

(48)

ことは困難であるため，モンテカルロ法に基づく近似計算(49)

法が提案されている [11, 24]．

このとき，L を単調増加させるには，各因子に関する変分
事後分布を次式に従って順番に更新すればよい．

算すると，最終的に以下を得る．

q(θk )=GIG(θk |γkθ , ρθk , τkθ )q(θ) ∝ p(θ) exp(Eq(H ,V ) [log q(X|θ, H, V )])q(H) ∝ p(H) exp(Eq(θ ,V ) [log q(X|θ, H, V )])上記議論を踏まえて，式 (50) に従って変分事後分布を計hh, ρhkn , τkn)
q(hkn )=GIG(hkn |γkn

(50)q(V ) ∝ p(V ) exp(Eq(θ ,H ) [log q(X|θ, H, V )])ここで，log q(X|θ, H, V ) は log p(X|θ, H, V ) の変分下限であり，式 (47) から期待値計算を省いたものである．

このときの変分パラメータは以下で与えられる．

γkθ =αc/Kρθk =2α + ν

ルに対して変分ベイズ法を用いて変分事後分布を求めるた
めには，事前分布と尤度関数との共役性を確立する必要がある．しかし，θk に着目すると，ガンマ事前分布の対数はθk および log θk に関する項を含んでいるが，式 (47) で与えられる対数尤度の変分下限は θk および θk−1 に関する項を含んでいる．このような不一致があると，事後分布の形
は事前分布と同じガンマ分布にはならない．
本研究では，q(θk ) および q(hkn ) は一般化逆ガウス (Gen-eralized Inverse Gaussian: GIG) 分布となることが分かる．
非負の実数上に定義された GIG 分布は次式で与えられる．

c 2013 Information Processing Society of JapanNtr E[hkn V k ] Ω−1nn=1

5.3 変分事後分布の計算
式 (40)，式 (41) および式 (43) で定義されたベイズモデ(54)

q(V k )=MGIG(V k |γkV , RVk , T Vk )τkθ =ν

Nn=1hγkn =a0−1
Φkn X n ΦTkn
tr E h−1
kn V kρhkn =2b0 + νtr E[θk V k ] Ω−1n −1 −1 hτkn =νtr E θk V k Φkn X n ΦTkn

(55)

γkV =ν0 /2
RVk =V −1
0 +ν

N


E[θk hkn ] Ω−1
n

n=1
N



T
T Vk =ν
E θk−1 h−1
kn Φkn X n Φkn
n=18Vol.2013-MUS-100 No.8
2013/8/31

情報処理学会研究報告
IPSJ SIG Technical Report
30 [dB]
3.7 dB25

4.1 dB

1.8 dB

反復STFT法を用いると音源分離品質は低下した

3.9 dB20

1.3 dB
1.4 dB

C4

SDR

SIR

E4

0

KL-NMF
KL-NMF with 反復STFT
IS-NMF
IS-NMF with 反復STFT
LD-PSDTF (提案法)1510

与えられた混合音

LD-PSDTF は KL-NMF や
IS-NMF よりも有意に高い性能を示した

G4

2

C4+E4

4

観測行列

C4+G4

6

8

観測行列

E4+G4

C4+E4+G4

10

12

観測行列

14[s]

観測行列

SAR

図 3 音源分離実験結果

6. 評価
LD-PSDTF を用いた音源分離実験について報告する．
本稿では，4 章で述べた乗法更新アルゴリズムに基づく最尤推定を用いるものとし，5 章で述べた変分ベイズ法に基づくベイズ推定については [11] を参考にされたい．

Log-Determinantダイバージェンスを最小化するような基底行列の凸結合を見つける
基底行列

基底行列

128

128

128

256

256

256

384

384

C4の周期
C4の周期

512

128

0

256

384

512

音量変化

0.10

6.1 実験条件

基底行列

0

2

4 6

実験には，RWC 研究用音楽データベース：楽器音 [25]図 4に収録されているピアノ (011PFNOM)，エレキギター(131EGLPM) およびクラリネット (311CLNOM) の単独音

512

128

0

256

384

512

音量変化

0.14

8 10 12 14

384

E4の周期
E4の周期

0

2

4 6

G4の周期
G4の周期

512

128

8 10 12 14

0

256

384

512

音量変化

0.35

0

2

4 6

2

8 10 12 14 x10

LD-PSDTF によるピアノ信号の分解結果

基底行列

基底行列

基底行列

を用いた．各楽器ごとに，異なる 3 つの音高 (C4, E4, G4)をもつ 2 秒間の音響信号を準備し，それらを 7 つの異なる組み合わせで重畳したもの (C4, E4, G4, C4+E4, C4+G4,

E4+G4, C4+E4+G4) を連結することで 14 秒の音響信号

128

128

128

256

256

256

384

384

384

512

512

128

を合成した．サンプリング周波数は 16[kHz] とした．
次に，LD-PSDTF を用いて，与えられた混合音を C4, E4,

G4 に対応する音源信号に分離することを試みた．まず，ガ

384

512

128

256

384

512

音量変化

音量変化

0.06

0

256

0.06

0

2

4 6

8 10 12 14

0

0

2

4 6

8 10 12 14

512

128

0.08

0

256

384

512

音量変化

0

2

4 6

2

8 10 12 14 x10

T
ウス窓を用いて局所信号 {on }N
n=1 を切り出し，X n = on onとすることで観測データ {X n }N
n=1 を得た．窓幅は 512 点，

図 5 LD-PSDTF によるクラリネット信号の分解結果

シフト長は 160 点とした (M = 512, N = 1400)．比較のたながらないことを意味する．図 4 および図 5 に示す通り，

め，振幅スペクトログラムに対する KL-NMF とパワースペ

LD-PSDTF を用いると，減衰音と持続音のいずれに対してクトログラムに対する IS-NMF も評価した（2.3 節およびも基底行列 V および音量変化 H を適切に推定することが2.4 節）．位相情報を復元するため，反復 STFT 法 [6] を用できた．ここで，各基底行列 V k における斜め縞の間隔はいる場合も評価した．各手法に対して，基底数は K = 3，反周期を表しており，V k の中心付近は循環行列に近くなっ

復回数は 100 回とし，乗法更新アルゴリズムを用いた．音ている．しかし，窓関数の影響で周辺部はそうはならない

源分離結果は，BSS Eval Toolbox [26] を用いて，source-to-ので，3.4 節で議論した通り，周波数領域においては周波

distortion ratio (SDR)，source-to-interferences ratio (SIR)

数成分間に相関が発生することは原理的に避けられない．

および sources-to-artifacts ratio (SAR) で評価した．

LD-PSDTF はこの影響を考慮することで優れた音源分離

6.2 実験結果
図 3 に示す通り，音源分離タスクにおいて LD-PSDTF

性能を達成している．また，観測行列である {X n }N
n=1 は格子状の分布をもっており，斜行型の分布をもつ基底行列{V k }K
k=1 の凸結合では十分な近似となっていないように思は NMF に対する明確な優位性を示した．SDR, SIR, SAR

える．しかし，LD ダイバージェンスは Y n が X n を過大の平均は，KL-NMF では 17.7 dB, 22.2 dB, 19.7 dB，IS-

評価しても小さなコストしか課さないことから，妥当な分NMF では 19.1 dB, 24.0 dB, 21.0 dB であったのに対し，解となっている．LD-PSDTF の主な課題は，計算コスト

LD-PSDTF では 23.0 dB, 27.7 dB, 25.1 dB であった*3 ．

が O(KN M 3 ) であり，NMF の O(KN M ) よりもはるかに

反復 STFT 法を用いると，音源信号の品質は低下した．こ

大きいことである．計算時間を短縮し，局所解を回避するのことは，[10] に示唆されている通り，スペクトログラム

ため，実際には IS-NMF で LD-PSDTF を初期値する，すの無矛盾性の向上が必ずしも音響信号の品質の向上につなわちある程度収束が進むまで基底行列を対角行列に制限*3

して反復更新を行う方法が推奨される．

著者の WEB サイトにおいてサンプルファイルが試聴可能.

c 2013 Information Processing Society of Japan

9

Vol.2013-MUS-100 No.8
2013/8/31

情報処理学会研究報告
IPSJ SIG Technical Report

7. おわりに[10]

本稿では，Log-Determinant 半正定値テンソル分解 (LD-PSDTF) と呼ぶ新しい因子分解法を提案した．LD-PSDTF
は IS-NMF の自然な拡張となっていることを示し，IS-NMF[11]と同様に収束性が保証された乗法更新則を用いて最尤推定が行えることを示した．また，ガンマ過程に基づくノンパ[12]

ラメトリックベイズモデルを提案し，観測データに合わせて基底数を自動調節する枠組みについて論じた．LD-PSDTF[13]

を音源分離に適用すると，KL-NMF や IS-NMF より優れた音源分離結果が得られることを実験で確認した．
今後の研究にはいくつかの興味深い方向性が考えられる．

[14]

まず，IS-NMF の拡張である複合自己回帰モデル [27,28] と同様に，LD-PSDTF の枠組みにソース・フィルタ理論を組み込むことで，音高ではなく音色に基づく高品質な音源分離の実現に取り組みたい．さらに，PSDTF の興味深いもう[15]

一つの変種である von Neumann (vN) ダイバージェンスに基づく PSDTF （式 (25) において φ(Z) = tr(Z log Z − Z)となるとき）に対して，最尤推定法およびベイズ推定法を[16]

確立したい．vN-PSDTF は KL-NMF の自然な拡張となっているため，幅広い分野への応用が期待できる．
謝辞: 本研究の一部は，JSPS 科研費 23700184, MEXT 科研
費 25870192, JST CREST OngaCREST の支援を受けた．

参考文献
[1]

[2]
[3]

[4]

[5]
[6]

[7]

[8]

[9]

K. Itoyama, M. Goto, K. Komatani, T. Ogata, and H. G.
Okuno. Query-by-example music information retrieval
by score-informed source separation and remixing technologies. EURASIP Journal on Advances in Signal Processing, Article ID 172961, 2010.
M. Goto. Active music listening interfaces based on signal processing. ICASSP, volume 4, pp. 1441–1444, 2007.
K. Yoshii, M. Goto, K. Komatani, T. Ogata, and H. G.
Okuno. Drumix: An audio player with real-time drumpart rearrangement functions for active music listening.
IPSJ Digital Courier, 3:134–144, 2007.
N. Sturmel, A. Liutkus, J. Pinel, L. Girin, S. Marchand,
G. Richard, R. Badeau, and L. Daudet. Linear mixing
models for active listening of music productions in realistic studio conditions. AES Convention, 2012.
D. Lee and H. Seung. Algorithms for non-negative matrix factorization. NIPS, pp. 556–562, 2000.
D. W. Griﬃn and J. S. Lim. Signal estimation frommodiﬁed short-time Fourier transform. IEEE Trans.
on Acoustics, Speech, and Signal Processing, 32(2):236–
243, 1984.
J. Le Roux, H. Kameoka, N. Ono, and S. Sagayama.
Explicit consistency constraints for STFT spectrogramsand their application to phase reconstruction. SAPA,
pp. 23–28, 2008.
J. Le Roux, H. Kameoka, N. Ono, and S. Sagayama.
Fast signal reconstruction from magnitude STFT spectrogram based on spectrogram consistency. DAFx, pp.
397–403, 2010.
H. Kameoka, T. Nishimoto, and S. Sagayama. Complex
NMF: A new sparse representation for acoustic signals.

c 2013 Information Processing Society of Japan


[17]

[18]

[19]

[20]

[21]

[22]
[23]

[24]

[25]

[26]

[27]

[28]

ICASSP, pp. 45–48, 2009.
J. Le Roux, E. Vincent, Y. Mizuno, H. Kameoka, N. Ono,and S. Sagayama. Consistent Wiener ﬁltering: Generalized time-frequency masking respecting spectrogram
consistency. LVA/ICA, pp. 89–96, 2010.
K. Yoshii, R. Tomioka, D. Mochihashi, and M. Goto. Inﬁnite positive semideﬁnite tensor factorization for source
separation of mixture signals. ICML, pp. 576–584, 2013.
P. Smaragdis and J. C. Brown. Non-negative matrix factorization for polyphonic music transcription. WASPAA,
pp. 177–180, 2003.
C. Févotte, N. Bertin, and J.-L. Durrieu. Nonnegative
matrix factorization with the Itakura-Saito divergence:With application to music analysis. Neural Computation, 21(3):793–830, 2009.
L. M. Bregman. The relaxation method of ﬁnding
the common points of convex sets and its application
to the solution of problems in convex programming.
USSR Computational Mathematics and Mathematical
Physics, 7(3):200–217, 1967.
M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono,and S. Sagayama. Convergence-guaranteed multiplicative algorithms for non-negative matrix factorization
with beta divergence. MLSP, pp. 283–288, 2010.
A. T. Cemgil. Bayesian inference for nonnegative matrix
factorisation models. Computational Intelligence and
Neuroscience, 2009:Article ID 785152, 2009.
M. Hoﬀman, D. Blei, and P. Cook. Bayesian nonparametric matrix factorization for recorded music. ICML,
pp. 439–446, 2010.
J. B. Allen and L. R. Rabiner. A uniﬁed approachto short-time Fourier analysis and synthesis. IEEE,
65(11):1558–1564, 1977.
B. Kulis, M. Sustik, and I. Dhillon. Low-rank kernel
learning with Bregman matrix divergences. Journal of
Machine Learning Research (JMLR), 10:341–376, 2009.
H. Sawada, H. Kameoka, S. Araki, and N. Ueda. Eﬃcientalgorithms for multichannel extensions of Itakura-Saito
nonnegative matrix factorization. ICASSP, pp. 261–264,
2012.
O. Barndorﬀ-Nielsen, P. Blæsild, J. L. Jensen, andB. Jørgensen. Exponential transformation models. Royal
Society of London, 379(1776):41–65, 1982.
C. S. Herz. Bessel functions of matrix argument. Annalsof Mathematics, 61(3):474–523, 1955.
R. W. Butler. Generalized inverse Gaussian distributions
and their Wishart connections. Scandinavian Journal of
Statistics, 25(1):69–75, 1998.
M. Yang, Y. Li, and Z. Zhang. Multi-task learning with
Gaussian matrix generalized inverse Gaussian model.
ICML, pp. 423–431, 2013.
M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka.
RWC music database: Music genre database and musical instrument sound database. ISMIR, pp. 229–230,
2003.
E. Vincent, R. Gribonval, and C. Févotte. Performancemeasurement in blind audio source separation. IEEE
Trans. on Acoustics, Speech, and Signal Processing,
14(4):1462–1469, 2006.
H. Kameoka and K. Kashino. Composite autoregressive
system for sparse source-ﬁlter representation of speech.
ISCAS, pp. 2477–2480, 2009.
K. Yoshii and M. Goto. Inﬁnite composite autoregressive models for music signal analysis. ISMIR, pp. 79–84,
2012.

10