社団法人 電子情報通信学会
THE INSTITUTE OF ELECTRONICS,
INFORMATION AND COMMUNICATION ENGINEERS

信学技報
TECHNICAL REPORT OF IEICE.

［ポスター講演］音楽音響信号解析のためのステューデント t 分布に基づく非負値行列分解と半正定値テンソル分解
吉井

和佳†

糸山

克寿†

後藤

真孝††

† 京都大学 大学院情報学研究科 〒 606–8501 京都府京都市左京区吉田本町†† 産業技術総合研究所 情報技術研究部門 〒 305–8568 茨城県つくば市梅園 1–1–1
E-mail: †{yoshii,itoyama}@kuis.kyoto-u.ac.jp, ††m.goto@aist.go.jp
あらまし 本稿では，裾が重い音響信号の音源分離に適した非負値行列分解 (NMF) と半正定値テンソル分解 (PSDTF)の新しい確率モデルについて述べる．従来，板倉・斎藤ダイバージェンスに基づく NMF (IS-NMF) を用いると，音源信号の複素スペクトル中の各要素が平均 0 の複素ガウス分布に従うという仮定のもとで，混合音のパワースペクト
ルを少数の音源信号のパワースペクトルの錐結合に分解できることが知られている．しかし，実際には，音源信号は裾が重い確率分布に従う場合が多いため，複素ガウス分布の代わりに複素コーシー分布に基づく Cauchy NMF が提案されている．本研究では，IS-NMF や Cauchy NMF を特殊形に含む，複素 t 分布に基づく NMF (t-NMF) を提案
し，統一的な学習アルゴリズムを導出する．t-NMF は対象となるデータの裾の重さに合わせて複素 t 分布の自由度を調節することで，IS-NMF と Cauchy NMF の中間的な性質をもつ NMF を実現できる．さらに，NMF の拡張である PSDTF に対しても同様の拡張を行い，多変量複素 t 分布に基づく PSDTF (t-PSDTF) を提案する．実験の結果，

t-NMF は IS-NMF より初期値依存性が低く，自由度を小さく抑えることで，安定して高精度な音源分離を行えること，t-PSDTF は t-NMF より精度の高い分離が行えることを確認した．
キーワード 音楽音響信号解析，非負値行列分解，半正定値テンソル分解，音源分離，t 分布[Poster Presentation] Nonnegative Matrix Factorization
and Positive Semidefinite Tensor Factorization
based on Student’s t Distributions for Music Signal Analysis
Kazuyoshi YOSHII† , Katsutoshi ITOYAMA† , and Masataka GOTO††† Kyoto University, Yoshida-honmachi, Sakyo-ku, Kyoto 606-8501, Japan†† National Institute of Advanced Industrial Science and Technology (AIST)
E-mail: †{yoshii,itoyama}@kuis.kyoto-u.ac.jp, ††m.goto@aist.go.jp
Abstract This paper presents new probabilistic models of nonnegative matrix factorization (NMF) and positive
semidefinite tensor factorization (PSDTF) for source separation of heavy-tailed audio signals. It is well known that
NMF based on the Itakura-Saito divergence (IS-NMF) is capable of decomposing the power spectrum of a mixture
signal into the weighted sum of source spectra. Real audio signals, however, tend to follow heavy-tailed probability
distributions such as the Cauchy distribution. In this paper we propose NMF based on the complex t distribution(t-NMF) that includes IS-NMF and Cauchy NMF. In addition, we propose PSDTF based on the multivariate complex t distribution (t-PSDTF). The experimental results showed that t-NMF was less sensitive to initialization than
IS-NMF and t-PSDTF significantly outperformed t-NMF in terms of source separation performance.
Key words Music signal analysis, NMF, PSDTF, source separation, t distribution

1. は じ め に
モノラル音響信号の音源分離においては，非負値行列分解とウィーナーフィルタを用いることが標準的になりつつある [1–9]．

NMF では，入力となる非負値行列（非負値ベクトルの集合）
を，二つの非負値行列（基底ベクトルの集合と音量ベクトルの—1—集合）の積で近似する，すなわち，各非負値ベクトルを少数の非負値基底ベクトルの重み付き和で近似する．NMF を音楽音基底行列の錐結合で近似: フレームの複素スペクトル

対角成分のみに着目（非負値のパワースペクトル）

響信号の非負値スペクトログラム（振幅あるいはパワースペク
トログラム）に適用する場合には，得られる基底ベクトルは異なる音高をもつ音源信号のスペクトルに対応することが期待さ[dB]

れる．最後に，NMF の結果を用いて，時間・周波数ビンごとに独立にウィナーフィルタリングを行うことにより，混合音ス
ペクトルを音源スペクトルの和に分解できる．しかし，音源ス
ペクトルの位相情報は混合音スペクトルのものをそのまま再利

半正定値基底行列の集合

基底ベクトルの錐結合で近似非負値基底ベクトルの集合

図 1 PSDTF と NMF との比較

この疑問を解消するため，最近，複素コーシー分布に基づく用することになり，逆フーリエ変換で得られる時間領域の音源

NMF (Cauchy NMF) が提案されている [17]．具体的には，音信号の品質に限界があることが指摘されていた [10]．

源信号の複素スペクトルが複素コーシー分布に従うと仮定する

近年，この問題を回避するため，半正定値テンソル分解 (PS-と，複素コーシー分布も再生性を持つことから，それらが重畳

DTF) が提案されている [11, 12]．PSDTF では，入力としてした混合音の複素スペクトルも複素コーシー分布に従う（注 1）．

半正定値行列の集合が与えられると，各半正定値行列を少数のこのとき，近似的に振幅スペクトルの加法性が成立する．また，

半正定値基底行列の重み付き和で近似する．図 1 で示す通り，

複素コーシー分布は裾が重い分布であるため（平均・分散とも行列の半正定値性はスカラやベクトルの非負値性の拡張概念でに定義されない），外れ値に頑健であるという性質を持っておあることから，PSDTF は NMF の自然な拡張となっている．

り，実在する音源信号をモデル化するのに都合が良い．

PSDTF を音源分離に用いる場合は，入力となる半正定値行列

より一般には，音源信号の複素スペクトルが複素対称 α 安定

は，各フレームにおける混合音の複素スペクトルとその共役ス(SαS) 分布に従うと仮定すると，混合音の複素スペクトルも複ペクトルとの直積となる．NMF では半正定値行列の対角成分，

素 SαS 分布に従う，すなわち Fractional スペクトル（注 2） のすなわち非負値ベクトルのみを対象としており，周波数ビン間加法性が成立する [18]．ここで，α は特性指数と呼ばれる分布の相関が無視されていることに注意する．一方，PSDTF ではの裾の重さを調節するパラメータ (0 < α <= 2) で，α = 1, 2

位相情報を取り扱うことができるため，各フレームごとに独立のときそれぞれ複素コーシー分布および複素ガウス分布に帰着にウィナーフィルタリングを行う（当該フレームの周波数ビンする．この生成モデルのもとで，事後分布推論の観点から，混は相関を考慮しながら一挙に処理される）ことにより，音源信

合音スペクトルを音源スペクトルに分解するための一般化ウィ号の複素スペクトルを直接推定可能である．

ナーフィルタリングが導ける．しかし，複素 SαS 分布の確率密

NMF には様々な変種が存在するが，ある仮定のもとで，板度関数は α = 1, 2 のときを除いて解析的に表現することができ

倉・斎藤 (IS) ダイバージェンスに基づく NMF (IS-NMF) [13]ないため，一般の α について，複素 SαS 分布に基づく NMF（あるいはその拡張である LogDet ダイバージェンスに基づくの最尤推定アルゴリズムは導出されていない．現在のところ，

PSDTF (LD-PSDTF) [11]）を用いて音源分離を行うことが理

その特殊形である IS-NMF [13] と Cauchy NMF [17] に対する論的に妥当となる．具体的には，音源信号の複素スペクトルが最尤推定アルゴリズムが個別に提案されているのみである．

複素ガウス分布に従うと仮定すると，複素ガウス分布は再生性

本稿では，確率密度関数が解析的に表現できる複素 t 分布に基を持つことから，それらが重畳した混合音の複素スペクトルもづく NMF (t-NMF) を提案する．複素 t 分布も自由度と呼ばれ複素ガウス分布に従う．したがって，パワースペクトルの加法る裾の重さを調節するパラメータ ν (0 < ν) を持ち，ν = 1, ∞性が成立し，混合音スペクトルに対する複素ガウス尤度の最大のときそれぞれ複素コーシー分布および複素ガウス分布に帰着化が，IS ダイバージェンスの最小化と等価となる．

する．ただし，それ以外の一般の ν に対しては，分布の再生性

実際には，カルバック・ライブラー (KL) ダイバージェンス

は成立しない．したがって，スペクトルの加法性も一般に成立に基づく NMF (KL-NMF) を用いて音源分離を行うことが一しないが，複素 SαS 分布の軸とは異なるものの，Cauchy NMF

般的である [14]．具体的には，音源信号の振幅スペクトルがポと IS-NMF の中間的な性質を持つ NMF を実現することは実用

アソン分布に従うことを仮定すると，ポアソン分布は再生性を面と学術面で有意義である．まず，理論的にスペクトルの加法持つことから，混合音の振幅スペクトルもポアソン分布に従う．

性が保証されなくても，実験的に音源分離精度を最大化するよ

したがって，振幅スペクトルの加法性が成立し，混合音スペク

う自由度 ν を最適化できる．また，Cauchy NMF と IS-NMF

トルに対するポアソン尤度の最大化が，KL ダイバージェンス

を含む統一的な最尤推定アルゴリズムを導出することにより，の最小化と等価となる．このモデルでは，振幅スペクトルは度

学習時の振る舞いを明らかにすることができる．

数表，すなわち音粒子のヒストグラムであるという物理的に意（注 1）：中心極限定理によると，独立同分布の確率変数を考えると，その確率分

味を持たない解釈がなされていることに注意する．にもかかわ布がなんであれ，確率変数の和の確率分布は変数が多くなるにつれて，ガウス分

らず，実験的には KL-NMF が IS-NMF よりも良い結果を与え

布に収束する．ただし，これが成立するのは確率変数の分散が有限の場合のみでることが多いことが知られている [15, 16]．振幅スペクトルのあり，確率変数が従う分布の裾が |x|−α−1 (0 < α < 2) の冪乗で減衰（分散は無限大に発散）する場合は，特性指数 α の安定分布に収束する．

加法性が成立するように，混合音に対して物理的に意味のある（注 2）：複素スペクトルに対して，各要素ごとに絶対値を取り，α 乗したもの．

生成モデルを考えられるだろうか？

α = 1 のとき振幅スペクトル，α = 2 のときパワースペクトルとなる．

—2—

本稿ではさらに，多変量 t 分布に基づく PSDTF (t-PSDTF)を提案する．t-PSDTF も t-NMF と同様に，ν = 1, ∞ のときが，多変量 t 分布の再生性は ν = 1 のときも成立しなくなるこ

周波数

Cauchy PSDTF（本稿で提案）と LD-PSDTF [12] に帰着する

観測 と再構成との距離を最小化するように と を反復最適化と注意する．スペクトルの加法性が成立するのは ν = ∞ である LD-PSDTF のみであるが，実際のデータに対して ν を実験
時間的に最適化できることは実用上重要である．t-PSDTF に対する最尤推定アルゴリズムは，Cauchy NMF [17], IS-NMF [13],図 2 IS-NMF による混合音スペクトログラムの低ランク分解

Cauchy PSDTF, LD-PSDTF [12] をすべて特殊形に含む統一的な乗法更新則となっている．Cauchy NMF については，[17]2. 2 コスト関数最小化としての LD-PSDTF

で提案されているものとは異なり，見通しが良く，収束性が保PSDTF の目標は，半正定値行列 （注 3） の集合（テンソル）
X = [X 1 , · · · , X N ] ∈ CM ×M ×N に対して，少数の半正定値証された新しい学習アルゴリズムを与える．

2. 複素ガウス分布に基づく NMF と PSDTF
本章では，t-PSDTF を定式化するうえで必要となる数学的な背景について説明する．まず，コスト関数最小化の立場から，

NMF [13, 19] と PSDTF [11, 12] の一般的な定式化を説明し，行列の集合（テンソル）W = [W 1 , · · · , W K ] ∈ CM ×M ×K お
N ×Kの積 Y = [Y 1 , · · · , Y N ] ∈よび H = [h1 , · · · , hK ] ∈ R+

CM ×M ×N で近似することである．
def

X ≈ W ×3 H = Y(4)PSDTF が NMF の自然な拡張になっていることを示す．NMF

ここで，×3 は 3 モード積 [20] であり，W k ∈ CM ×M  0 おと同様に PSDTF にも様々なコスト関数が考えうるが，LogDet

よび hk ∈ RN
+ はそれぞれ基底行列およびアクティベーション(LD) ダイバージェンスに基づく PSDTF (LD-PSDTF) は，板ベクトルである．式 (4) は次式で書き直せる．

倉・斉藤 (IS) ダイバージェンスに基づく NMF (IS-NMF) の数学的拡張となっている．IS-NMF や LD-PSDTF を音源分離に適用する例を通して，コスト関数の最小化は，（多変量）複素

Xn ≈

K


def

hkn W k = Y n

(5)

k=1

ガウス分布で与えられる尤度関数の最大化（確率モデルの最尤

すなわち，PSDTF の目標は，半正定値行列を少数の半正定値推定）と等価であることを示す．

行列の錐結合で近似することである．
観測行列 X n と再構成行列 Y n との間の誤差 D(X n |Y n ) を

2. 1 コスト関数最小化としての IS-NMF
NMF の目標は，与えられた非負値行列 X = [x1 , · · · , xN ] ∈×N
M ×K
を，二つの非負値行列 W = [w1 , · · · , wK ] ∈ R+
RM
+
N ×K
および H = [h1 , · · · , hK ] ∈ R+の積 Y = [y 1 , · · · , y N ] ∈M ×Nで近似することである（図 2）．
R+

X ≈ W ×2 H = W H

T def

=Y

(1)

ここで，×2 は 2 モード積 [20] であり，wk ∈ RM
+ および

h k ∈ RN
+ はそれぞれ基底ベクトルおよび対応するアクティベー

評価する尺度として，von-Neumann ダイバージェンス [21] やLogDet ダイバージェンス [22] を利用することが考えられる．
本稿では，特に後者に着目する．





−M
DLD (X n |Y n ) = tr X n Y −1
− log X n Y −1
n
n

(6)

この値は常に非負であり，xn = y n のときのみ 0 をとる．通常の距離尺度と異なり，対称性が成り立たない．
式 (1) と式 (4)，式 (2) と式 (5)，式 (3) と式 (6) を比較する

ションベクトルである．K  min(M, N ) は基底数であり，予と，PSDTF は NMF の自然な拡張であることが明らかである

め定めておく必要がある．式 (1) は次式で書き直せる．

（図 1）．半正定値の対角行列では，対角成分が固有値に対応

xn ≈

K


し，それらはすべて非負値であることから，X n = diag(xn )，def

hkn wk = y n

(2)

k=1

Y n = diag(y n )，W k = diag(wk ) とすれば（注 4），PSDTF は
NMF に帰着する．すなわち，NMF では PSDTF のように，M

すなわち，NMF の目標は，非負値ベクトルを少数の非負値ベ個の要素間の相関を考慮することができない．

2. 3 尤度最大化としての IS-NMF と LD-PSDTF

クトルの錐結合で近似することである．
観測ベクトル xn と再構成ベクトル y n との間の誤差 D(x|y)

モノラル音響信号の音源分離を例に，確率モデルの最尤を評価する尺度には，KL ダイバージェンス [19] や IS ダイバー推定の観点から IS-NMF と LD-PSDTF を捉え直す．まず，ジェンス [13] などが知られている．本稿では，後者に着目する．

観測される音響信号（混合音）の複素スペクトログラムを

DIS (xn |y n ) =

M 

m=1

xnm
xnm
− log
−1
ynm
ynm

S = [s1 , · · · , sN ] ∈ CM ×N ，k 番目の音源信号の複素スペクト(3)

この値は常に非負であり，xn = y n のときのみ 0 をとる．
通常の距離尺度と異なり，対称性が成り立たない，すなわちDIS (xn |y n ) =| DIS (y n |nn ) であることに注意する．

ログラムを S k = [sk1 , · · · , skN ] ∈ CM ×N とする．ここで，M（注 3）：半正定値行列と非負値行列とは異なる概念である．ある複素エルミート
行列 A = AH ∈ CM ×M が半正定値行列である ⇔ A = V V H となる複素行列 V が存在すること ⇔ A の M 個の固有値が全て非負であること．
（注 4）：diag(a) はベクトル a を対角成分に持つ対角行列を意味する．

—3—

は周波数ビン数，N はフレーム数である．観測した混合音が K

ここで，P k および Qk は次式で求まる．

個の音源信号の瞬時混合であると仮定すると，以下が成り立つ．

S=

K


S k ⇔ sn =

k=1

K


Pk =
skn

(7)

N


hkn Y −1
Qk =
n

n=1

N


−1
hkn Y −1n X nY n

(13)

n=1

hkn の非負性と W k の半正定値性は自然に保たれているが，k=1

我々の目標は，観測変数 S を潜在変数 S k の和に分解すること

tr(W k ) = 1 を満たすよう，反復のたびに W k および hk のスである．しかし，これは不良設定問題であるので，何らかの制

ケールを調整しておく．

約を導入して最適解を同定する必要がある．

すべての半正定値行列を対角行列に限定すれば，よく知られ我々は，各音源信号は局所的に定常なガウス過程に従う，すた IS-NMF の乗法更新アルゴリズム [23, 24] に帰着する．具

なわち，各フレーム n における局所連続信号は定常であり，時体的には，式 (11) および式 (12) において，X n = diag(xn )，不変な共分散関数（カーネル）をもつガウス過程に従うと仮定

Y n = diag(y n )，W k = diag(wk ) とすると以下を得る．

する．このとき，ガウス過程の任意の周辺分布は多変量ガウ
ス分布になることから，あるサンプリング周波数に従ってサンプルされた離散信号は多変量ガウス分布に従う．したがって，
フーリエ変換はアフィン変換であることから，各フレームにおける複素スペクトル skn は多変量複素ガウス分布に従う．

skn | Y kn ∼ N c (skn | 0, Y kn )

(8)

1

2 2
xnm hkn /ynmn
wkm ← wkm
n hkn /ynm
1

2 2
m xnm wkm /ynmhkn ← hkn
m wkm /ynm

(14)

(15)

2. 5 ウィナーフィルタリング
最終的に，式 (8) および式 (9) から，混合音スペクトル sn が

ここで，共分散行列を Y kn = hkn W k とした．すなわち，各与えられたときの音源スペクトル skn の事後分布は多変量複素音源信号において，W k で定まる基本的なダイナミクスは定常

ガウス分布となり，その平均と分散は次式で与えられる．

であるが，そのスケール hkn のみが時変化すると考えることで，音源信号の非定常性を制限する．
このとき，音源信号が重畳して得られる観測信号も局所的にE[skn |sn ] = Y kn Y −1
n sn

(16)

V[skn |sn ] = Y kn − Y kn Y −1
n Y kn(17)

定常なガウス過程に従う．式 (7) に着目すると，複素ガウス分

この処理では，sk の位相は s の位相とは異なる点に注意する．

布の再生性から以下を得る．

各フレーム n ごとに周波数ビン間の相関を考慮しながら一挙に

sn | Y n ∼ N c (sn | 0, Y n )ただし，Y n =

(9)

分離を行うことで，高品質な分離が可能となる．一方，IS-NMF
では式 (16) や式 (17) における半正定値行列はすべて対角行列k Y kn である．したがって，観測される混合

音スペクトル sn に対する対数尤度関数は以下で与えられる．

log p(sn |Y n ) = −M log(π) − log |Y n | − sn Y −1n sn= − log |Y n | − tr X n Y −1nc

ここで，X n = sn sHnに制限されているため，各周波数ビン n, m ごとに独立に分離が行われることになり，sk の位相は s の位相と同一となる．

3. 複素 t 分布に基づく NMF と PSDTF(10)

 0 とした．式 (10) を式 (6) と比較す

本章では，多変量複素ガウス分布の代わりに，より裾の重い
多変量複素 t 分布を尤度関数に用いることで，外れ値に対してると，X n 一定であるから，Y n に関する対数尤度の最大化は，頑健な PSDTF (t-PSDTF) を提案する．多変量複素 t 分布は，

LogDet ダイバージェンスの最小化と等価である．したがって，

多変量複素コーシー分布や多変量複素ガウス分布を特殊形とし式 (10) が LD-PSDTF の確率モデルとなっている．

て含む．複素 t 分布に基づく NMF (t-NMF) も同様に定式化

すべての半正定値行列を対角行列に限定すれば，LD-PSDTFすることができ，すべての半正定値行列を対角行列に限定すれは IS-NMF に帰着する．すなわち，式 (8) や式 (9) で与えられ

ば，t-PSDTF は t-NMF に帰着する．ただし，確率分布の再生

る多変量複素ガウス分布の分散行列は対角行列となり，音源ス性が成立するのは，一変量あるいは多変量複素ガウス分布と一ペクトル skn や混合音スペクトル sn の各要素が独立であり，変量複素コーシー分布のみである．この時に限り，すなわち，それぞれ異なる複素ガウス分布に従うことを意味する．

LD-PSDTF [12], IS-NMF [13], Cauchy NMF [17] に関しての

2. 4 乗法更新アルゴリズム

み，混合音スペクトルに対する理論的に妥当な生成モデルを考LD-PSDTF のパラメータ W および H を求めるため，補助えることができる．提案する t-PSDTF はこれらをすべて包含関数法に基づく収束性の保証された乗法更新アルゴリズムが提する一般的なモデルである．

案されている [11, 12]．本稿では結果のみ記すと，乗法更新則

3. 1 複素 α 安定分布と複素 t 分布は以下で与えられる．

1. 章で述べた通り，本来の目標は複素対称 α 安定 (SαS) 分

1

1

1

1

1

W k ← W k Qk2 (Qk2 W k P k W k Qk2 )− 2 Qk2 W k

(11)



−1
tr X n Y −1
n W kY ntr W k Y −1n(12)hkn ← hkn

12

布に基づく NMF や PSDTF の定式化であることを念頭に置
き，複素 SαS 分布と複素 t 分布の性質を整理しておく．一変量の場合と多変量の場合で，確率分布の再生性が成立する条件が異なることに注意する．

—4—

0.4

0.4

0.3

0.3

0.2

0.2

3. 1. 3 複素 α 安定分布と複素 t 分布との比較
複素 α 安定分布では α = 1, 2 の場合に，複素 t 分布では

ν = 1, ∞ の場合に，それぞれ複素コーシー分布と複素ガウス

0.1

0.1
-5

-4

-3

-2

-1

0

1

2

3

4

-5

5

-4

-3

-2

-1

0

1

2

3

4

5

Standard distribution:

S S distribution:

図 3 実数の確率変数に対する一変量対称 α 安定分布と t 分布

一変量複素 α 安定分布の確率密度関数を Sαc (μ, σ α ) と表す．
ここで，α は特性指数 (0 < α <= 2)，μ は位置 (location) パラ
メータ (μ ∈ C)，σ α は分散 (dispersion) パラメータ (σ > 0)である．特に μ = 0 のとき複素対称 α 安定 (complex symmet-

ric α-stable: complex SαS) 分布と呼ぶ．一般の α について，
Sαc (μ, σ α ) は解析的に書けないが，特性関数は簡潔に書ける．
ϕ(t) = E ec(μ, 2σ 2 ) =T1c (μ, σ 2 ) = C c (μ, σ 2 ) お よ び S2c (μ, σ 2 ) = T∞N c (μ, 2σ 2 )，多変量の場合は S1c (μ, Σ) = T1c (μ, Σ) = C c (μ, Σ)c
および S2c (μ, Σ) = T∞(μ, 2Σ) = N c (μ, 2Σ) が成立する．

3. 1. 1 複素 α 安定分布

iRe(t∗ X)

分布に帰着する．具体的には，一変量の場合は S1c (μ, σ) =  α= exp − 12 σ |t|α + iRe(t∗ μ) (18)参考のため，実数の確率変数に対する一変量の対称 α 安定分布と t 分布の確率密度関数を図 3 に示す．ここで，α 安定分布の分散パラメータは S2 (0, 0.5) = T∞ (0, 1) かつ S1 (0, 1) = T1 (0, 1)となるよう適当に調整されている. 裾の減衰の仕方は異なるものの，両分布ともに，ガウス分布とコーシー分布の中間的な特徴を持つ裾の重い分布を表現することができる．本研究では，
確率密度関数の取り扱いの容易さから複素 t 分布に着目する．

3. 2 尤度最大化としての定式化
PSDTF を頑健化するため，式 (9) で与えられる尤度関数に特性関数と確率密度関数はフーリエ変換の関係にあるので，二おいて，多変量複素ガウス分布の代わりに，自由度 ν の多変量

つの確率密度関数の畳込み（二つの確率変数の和の確率分布）

複素 t 分布を用いることを考える．

は，対応する特性関数の積と等価である．したがって，同じ特性指数の安定分布に従う確率変数の和は，同じ特性指数の安定分布に従う．すなわち，全ての α について再生性を持つ．

X1 ∼ Sαc (μ1 , σ1α )X2 ∼ Sαc (μ2 , σ2α )⇒ X1 +X2 ∼ Sαc (μ1 +μ2 , σ1α +σ2α ) (19)多変量複素 α 安定分布のうち，楕円型 (elliptically contoured)分布の確率密度関数は Sαc (μ, Σ) と表す．ここで，μ は位置ベ
クトル (μ ∈ Cd )，Σ は分散行列 (Σ  0 ∈ Cd×d ) である．この分布の特性関数は以下で与えられる．

 

α
H
ϕ(t) = E eiRe(t X) = exp − 12 tH Σt 2 + iRe(tHμ) (20)
このことから，多変量の場合には，α = 2 の多変量複素ガウス分布の場合に限り再生性を持つ．

sn | Y n ∼ Tνc (sn | 0, Y n )ここで，Y n =

(23)

K

k=1 hkn W k である．このとき，観測データ

X 全体に渡る対数尤度は次式で与えられる．
log p(X |Y) =
c

=

N 

n=1

N


log p(sn |Y n )

(24)

n=1





2 
ν
log 1 + tr X n Y −1
− log |Y n | − M +
n
2
ν

c

ここで，= は定数項を除いて等号が成立することを示す．
式 (24) を Y に関して直接最大化することは困難であるため，補助関数法（付録）を用いて式 (24) の下限関数を設計し，その下限関数を逐次最大化することで，式 (24) を間接的に最大化することを考える．

3. 1. 2 複素 t 分布
一変量複素 t 分布の確率密度関数を Tνc (μ, σ 2 ) と表す．ここで，ν は自由度 (ν > 0)，μ は位置パラメータ (μ ∈ C)，σ 2 は
分散パラメータ (σ > 0) である．一般の ν について，Tνc (μ, σ 2 )3. 2. 1 下限関数導出のための不等式
下 限 関 数 導 出 の た め に 必 要 な 不 等 式 を 説 明 す る ．ま ず，f (Z) = log |Z| が凹関数であることに着目すると，f (Z) に対して一次のテイラー展開を行うことで，次式を得る．

は解析的に書ける．

Tνc (x|μ, σ 2 ) =2Γ( 2+ν)2
πνσ 2 Γ( ν2 )



2 |x − μ|2
1+
ν
σ2

−1
− log |Z| >= − log |Ω| − tr(Ω Z) + M

− 2+ν
2(21)

したがって，ν → ∞ の複素ガウス分布あるいは ν = 1 の複素コーシー分布のときに限り再生性を持つ．
多変量複素 t 分布の確率密度関数を Tνc (μ, Σ) と表す．ここで，

ここで，Ω は任意の半正定値行列（展開点）であり，M は正方行列 Z の行数あるいは列数である．等号成立条件は，

Ω = Z で与えられる．次に，任意の半正定値行列 A に対して
g(Z) = tr(Z −1 A) は凸関数であることに着目すると，澤田ら

μ は位置ベクトル (μ ∈ Cd )，Σ は分散行列 (Σ  0 ∈ Cd×d )の提案する不等式 [25] を適用可能である．

であり，z = x − μ とすると，Tνc (μ, Σ) は以下で与えられる．


K
−1 



K
H
−tr
Z
A >
tr Z −1kk Φk AΦk=−
k=1

2d Γ( 2d+ν)2Tνc (x|μ, Σ) =νd(πν) Γ( 2 )|Σ|
1+2 H −1z Σ zν− 2d+ν
2(22)

多変量の場合は，ν → ∞ の多変量複素ガウス分布ときに限り再生性を持つ．

(25)

(26)

k=1
K
ここで，{Z k }K
k=1 は任意の半正定値行列であり，{Φk }k=1 は

足すと単位行列になるような補助行列である ( k Φk = I)．等
号成立条件は，Φk = Z k ( k Z k )−1 で与えられる．

—5—

3. 2. 2 下限関数の導出

t-PSDTF では，仮想的な観測行列 πn X n に対して，低ラン

これら不等式を用いると，式 (24) に対する補助関数は以下

でに πn X n には Y n の性質が反映されており，ある種の平滑の通り導出できる．

化が行われていると見ることができる．この結果，観測データ

log p(X |Y)
c

>
=

N


に対して過剰にフィットすることが抑制され，LD-PSDTF よりも外れ値に対して頑健な分解ができると考えられる．

− log |Ωn | − tr(Ω−1n Y n) + M

3. 4 t-NMF・IS-NMF・Cauchy NMF の比較

n=1







2 
ν
ψn + ψn−1 1 + tr X n Y −1
−1
− M+
n
2
ν
c

=

N


−1
−tr(Ω−1n Y n ) − ψn



n=1

>
=

N

n=1

− ψn−1多変量複素 t 分布に基づく t-PSDTF の特殊形として，複素t 分布に基づく NMF が得られる．2. 2 節で説明した通り，対角




2M
+ 1 tr X n Y −1nν

行列 X n = diag(xn )，Y n = diag(y n )，W k = diag(wk ) を考えると，以下の更新式を得る．

−tr(Ω−1n Y n)


クな再構成行列 Y n をフィットさせることが目的であるが，す2 2
nm xnm )hkn /ynm
n (π


(35)

2 2
nm xnm )wkm /ynm
m (π


(36)

wkm ← wkm


K


2M
−1
H
+1
tr h−1
kn W k Φkn X n Φkn
ν


(27)

hkn ← hkn

件（下限関数の最大化条件）は以下で与えられる．

(28)


2 
ψn = 1 + tr X n Y −1
n
ν

(29)

Φkn = Y kn Y −1n(30)

式 (27) の下限関数を W k および hkn に関してそれぞれ偏微
分してゼロとおくことで，以下の更新式を得る．
1

1

N




−1
tr πn X n Y −1
n W kY ntr W k Y −1nhkn Y −1
Qk =
n

n=1

N


ν
2
−1
ynm
x−1
+
nm
2+ν
2+ν

−1
(38)

−1
hkn Y −1n (πn X n )Y n

ν : 2 で調和平均をとることを意味し，自由度 ν が大きくなるほど，xnm を重視するようになり，ν が小さくなるほど，ynm を(32)

タ xnm と再構成データ ynm を 1 : 2 の重みで考慮する．

Cauchy NMF は ν = 1 の t-NMF と等価であるため，観測デー
t-NMF では，IS-NMF と Cauchy NMF に対して，見通しのよい統一的な乗法更新アルゴリズムを与える．特に，Cauchy

NMF に対する乗法更新アルゴリズムは，[17] で提案されてい
(33)

n=1

2M + ν
2tr(X n Y −1n )+νπnm xnm =と等価であるため，観測データ xnm のみに着目する．一方，

12

る二種類のアルゴリズムのいずれとも異なる．

4. 評 価 実 験

また，πn は次式で求められる．

πn =

このとき，仮想的な観測データ πnm xnm は次式となる．

(31)

ここで，P k および Qk は次式で求まる．

Pk =

(37)

重視するようになる．例えば，IS-NMF は ν → ∞ の t-NMF

1

1

W k ← W k Qk2 (Qk2 W k P k W k Qk2 )− 2 Qk2 W k
hkn ← hkn

2+ν
−1
2xnm ynm+ν

これは，観測データ xnm と再構成データ ynm に対して，重み3. 3 乗法更新アルゴリズム

11

ここで，πnm は次式で求まる．

πnm =Ωn = Y nn hkn /ynm

m wkm /ynm

k=1

ここで，Ωn ，ψn および Φkn は補助変数であり，等号成立条

1

(34)

t-NMF および t-PSDTF の音源分離精度や初期値依存性を確認するために行った評価実験について報告する．

したがって，hkn の非負性と V k の半正定値性は自然に保たれ4. 1 実 験 条 件ているが，tr(W k ) = 1 を満たすよう，反復ごとに W k および実験には，RWC 研究用音楽データベース：楽器音 [26] に収録

hk をスケーリングしておく．

されているピアノ (011PFNOM)，エレキギター (131EGLPM)

式 (31) および式 (32) は，式 (11) および式 (12) の自然な拡およびクラリネット (311CLNOM) の単独音を用いた．各楽器張となっている．実際，多変量複素 t 分布の自由度を ν → ∞ごとに，異なる 3 つの音高 (C4, E4, G4) をもつ 2 秒間の音響とすると，πn → 1 となり，多変量複素ガウス分布に基づく信号を準備し，それらを 7 つの異なる組み合わせで重畳したも

PSDTF (LD-PSDTF) に帰着する．一方，ν = 1 とすると，多の (C4, E4, G4, C4+E4, C4+G4, E4+G4, C4+E4+G4) を連

変量複素コーシー分布に基づく PSDTF が得られる．すなわち，結することで 14 秒の音響信号を合成した．サンプリング周波

提案する多変量複素 t 分布に基づく t-PSDTF では，πn X n を

数は 16[kHz] とした．まず，ガウス窓を用いて各フレームの複

仮想的な観測行列であるとみなして LD-PSDTF を行っているT
素スペクトル {sn }Nn=1 を計算し，半正定値行列 X n = sn snと解釈できる．ただし，πn は再構成行列 Y n に依存しているとすることで観測データ {X n }N
n=1 を得た．窓幅は 512 点，シ

ため，Y n が更新されるたびに観測行列 πn X n も変化する．

フト長は 160 点とした (M = 512, N = 1400)．

—6—

30

SDR

[dB]

SIR

観測される音響信号（混合音）

SAR
0

20

2

4

0.5 1 2 5 10 20 ∞ 1 ∞t -NMFt -PSDTF

0.5 1 2 5 10 20 ∞ 1 ∞t -NMF0.5 1 2 5 10 20 ∞ 1 ∞t -PSDTF

t -NMF

6

観測行列

8

観測行列

10

12

観測行列

14[s]

観測行列

Gaussian

Gaussian
CauchyCauchyPoissonCauchyGaussian

Gaussian

CauchyPoissonCauchyCauchyGaussian0Poisson10

Gaussian

観測行列

多変量複素t分布に基づく尤度を最大化するような基底行列を推定

t -PSDTF

基底行列

図 4 自由度を変化させた t-NMF と t-PSDTF による音源分離結果

次に，t-NMF および t-PSDTF を用いて，混合音を 3 つの音高に対応する音源信号に分離することを試みた．t-NMF で

基底行列

基底行列

128

128

128

256

256

256

384

384

Period of C4

512

128 256 384 512

Period of E4

512

128 256 384 512

アクティベーション

アクティベーション

データに用い，ν = 0.5, 1 (Cauchy NMF [17]), 2, 5, 10, 20,

00∞ (IS-NMF [13]) のときを評価した．t-PSDTF に関しては，

図5は X n の対角成分（非負値のパワースペクトル）のみを観測

2

4

6

8 10 12 14

2

4

6

384

Period of G4

512

128 256 384 512

アクティベーション

8 10 12 14 0

2

4

6

8 10 12 14

時間領域 t-PSDTF を用いたピアノ信号の分解結果

ν = 1 (Cauchy PSDTF), ∞ (LD-PSDTF [12]) のときを評価した．ここで，t-PSDTF は，多変量複素 t 分布がフーリエ変
換（線形変換）に対して閉じていることから，時間領域における等価な定式化が可能である [12]．このとき，ウィナーフィル
タリングを行うと，各音源信号の各フレームの複素スペクトルではなく，局所信号が直接求まる．また，振幅スペクトログラ
ムがポアソン分布に従うと仮定する KL-NMF とも比較した．
各手法はパラメータをランダムに初期化して，100 回試行した．計算量の削減のため，t-PSDTF は対応する t-NMF の平均的な結果を用いて初期化を行った（ν = 1, ∞ のとき）．音源分離結果は，BSS Eval Toolbox [27] を用いて，source-to-

distortion ratio (SDR)，source-to-interferences ratio (SIR)
および sources-to-artifacts ratio (SAR) で評価した．

4. 2 実 験 結 果
実験結果から，自由度を ν = 1, 2 程度に設定した t-NMF は，
従来の IS-NMF や KL-NMF より優れた分離精度を達成した．
どの自由度であっても t-NMF の分離精度の最大値は IS-NMFとほぼ同等であるが，自由度が ν = 1, 2 程度の t-NMF は初期値依存性が低く，安定して精度の良い音源分離が可能である
ことが分かった．ここで，自由度が ν = 1 の t-NMF (CauchyNMF) では，振幅スペクトルの加法性が成立するものの，自由度が ν = 2 の t-NMF の方がわずかに優れている傾向がみられることが興味深い．このことは，スペクトルの加法性が理論的に成立しなくても，実用上は最適な自由度は異なる場合があることを示している．また，ν の値を逆温度パラメータとみなし，少しずつ上昇させていくことで，よりよい局所解に誘導する一種のアニーリングも可能であると考えられる．

t-PSDTF は t-NMF に対する明確な優位性を示した．図 5 に示す通り，t-PSDTF を用いると，減衰音と持続音のいずれに対しても基底行列 W および音量変化 H を適切に推定できた．
ここで，各基底行列 W k における斜め縞の間隔は周期を表しており，W k の中心付近は巡回行列に近くなっている．しかし，5. お わ り に
本稿では，モノラル音響信号の音源分離を行うための新しい確率モデルである t-PSDTF およびその特殊形である t-NMFを提案した．実際の音響信号は優ガウス性を持つことが多いため，裾が重い確率分布である多変量複素 t 分布（特殊形として多変量複素ガウス/コーシー分布を含む）を尤度関数に仮定した．
このとき，パラメータを最尤推定するための効率的な乗法更新
アルゴリズムを導出した．結果的に，t-PSDTF は，従来知られていた LD-PSDTF [12]，IS-NMF [13]，Cauchy NMF [17]をその特殊形として含み，統一的な乗法更新アルゴリズムで学習が可能となった．特に，パワースペクトルの加法性に基づ
く IS-NMF と振幅スペクトルの加法性に基づく Cauchy NMFの中間的な性質をもつ NMF を実現したことの意義は大きい．
実験の結果，t-NMF は初期値依存性が低く，t-NMF および

t-PSDTF ともに高品質な音源分離ができることを確かめた．
今後は，多変量複素対称 α 安定分布を尤度関数とした SαS-PSDTF およびその特殊形である SαS-NMF の導出に取り組
みたい．特に，一変量の複素 SαS 分布は，複素ガウス/コー
シー分布を特殊な場合として含み，パワースペクトル (α = 2)や振幅スペクトル (α = 1) に限らず，あらゆる α に関して

Fractional スペクトルの加法性を理論的に保証できる点で画期的である [18]．ただし，確率密度関数は解析的に書くことがで
きず，最尤推定には特別な方法が必要になると考えられる．

付録

目的関数 F (θ) を変数 θ に関して最大化する問題を考える．
補助関数法 [28, 29] を用いると，F (θ) を「間接的に」最大化することができる．まず，F (θ) に対する補助関数として，++F (θ) >= F (θ, φ) となる下限関数 F (θ, φ) を考える．ここで，φ は補助変数である．このとき，以下の反復更新則

窓関数の影響で周辺部はそうはならないので，周波数成分間にφnew ← argmaxφ F + (θ old , φ)(A·1)

相関が発生することは原理的に避けられない．t-PSDTF はこθ new ← argmaxθ F + (θ, φnew )(A·2)の影響を考慮することで優れた音源分離性能を達成している．
また，t-PSDTF に対する乗法更新アルゴリズムが，ν の値にを用いると，F (θ) は単調非減少となる．このアルゴリズムの関わらず，安定して収束することを確認した．

収束性は保証されている.

—7—

多変量複素ガウス分布 （Q(x) = (x − μ)H Σ−1 (x − μ)）

N c (x|μ, Σ) =1
|Σ|−1 exp (−Q(x))
πd

[12]

多変量複素 t 分布


−
2d+ν
) 2d
2
−1
2
1
+
Q(x)
|Σ|
ν
Γ( 2 ) (πν)dνΓ(Tνc (x|μ, Σ) =2d+ν2

多変量複素コーシー分布

C c (x|μ, Σ) =

d 

i=1

[13]

[14]

i−

1
2



2d+1
2d
|Σ|−1 (1 + 2Q(x))− 2
πd[15]
謝辞: 本研究の一部は，JSPS 科研費 No. 24220006, 26700020,

26280089，JST CREST，栢森情報科学振興財団の支援を受けた．

文

[16]

献

[1] P. Smaragdis, C. Févotte, G. Mysore, N. Mohammadiha,and M. Hoﬀman. Dynamic source separation using nonnegative factorizations: A uniﬁed view. IEEE Signal Processing
Magazine, 31(3):66–75, 2014.
[2] A. T. Cemgil. Bayesian inference for nonnegative matrix
factorisation models. Computational Intelligence and Neuroscience, 2009:1–17, 2009.
[3] A. Ozerov, A. Liutkus, R. Badeau, and G. Richard. Codingbased informed source separation: Nonnegative tensor factorization approach. IEEE Transactions on Audio, Speech,
and Language Processing, 21(8):1699–1712, 2013.
[4] D. Liang, M. Hoﬀman, and D. Ellis. Beta process sparse
nonnegative matrix factorization for music. In International
Society for Music Information Retrieval Conference (ISMIR), pages 375–380, 2013.
[5] F. J. Canadas-Quesada, P. Vera-Candeas, N. Ruiz-Reyes,
J. J. Carabias-Orti, and P. Cabanas-Molero.
Percussive/harmonic sound separation by non-negative matrix factorization with smoothness/sparseness constraints.
EURASIP Journal on Audio, Speech, and Music Processing, 26(1):1–17, 2014.
[6] N. Souviraà-Labastie, E. Vincent, and F. Bimbot. Music
separation guided by cover tracks: Designing the joint NMF
model. In International Conference on Acoustics, Speech,
and Signal Processing (ICASSP), pages 484–488, 2015.
[7] N. Q. K. Duong D. El Badawy, A. Ozerov. Relative group
sparsity for non-negative matrix factorization with application to on-the-ﬂy audio source separation. In International Conference on Acoustics, Speech, and Signal Processing (ICASSP), pages 256–260, 2015.
[8] K. O’Hanlon, M. Sandler, and M. D. Plumbley. Matrix factorisation incorporating greedy Hellinger sparse coding applied to polyphonic music transcription. In International
Conference on Acoustics, Speech, and Signal Processing(ICASSP), pages 3112–3116, 2015.
[9] C.-W. Wu and A. Lerch. Drum transcription using partially ﬁxed non-negative matrix factorization with template
adaptation. In International Society for Music Information
Retrieval Conference (ISMIR), 2015.
[10] J. Le Roux, E. Vincent, Y. Mizuno, H. Kameoka, N. Ono,and S. Sagayama. Consistent Wiener ﬁltering: Generalized time-frequency masking respecting spectrogram consistency. In International Conference on Latent Variable
Analysis and Signal Separation (LVA/ICA), pages 89–96,
2010.
[11] K. Yoshii, R. Tomioka, D. Mochihashi, and M. Goto. Beyond NMF: Time-domain audio source separation withoutphase reconstruction. In International Society for Music

[17]

[18]

[19]

[20]
[21]

[22]

[23]
[24]
[25]

[26]

[27]

[28]

[29]

Information Retrieval Conference (ISMIR), pages 369–374,
2013.
K. Yoshii, R. Tomioka, D. Mochihashi, and M. Goto. Inﬁnite positive semideﬁnite tensor factorization for source
separation of mixture signals. In International Conference
on Machine Learning (ICML), pages 576–584, 2013.
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.
P. Smaragdis and J. C. Brown. Non-negative matrix factorization for polyphonic music transcription. In IEEE Workshop on Applications of Signal Processing to Audio andAcoustics (WASPAA), pages 177–180, 2003.
D. Fitzgerald, M. Cranitch, and E. Coyle. On the use of
the beta divergence for musical source separation. In Irish
Signals and Systems Conference (ISSC), pages 1–6, 2009.
B. King, C. Févotte, and P. Smaragdis. Optimal cost function and magnitude power for NMF-based speech separation and music interpolation. In International Workshop
on Machine Learning for Signal Processing (MLSP), pages
1–6, 2012.
A. Liutkus, D. Fitzgerald, and R. Badeau. Cauchy nonnegative matrix factorization. In IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA), 2015.
A. Liutkus and R. Badeau. Generalized Wiener ﬁltering with fractional power spectrograms. In International
Conference on Acoustics, Speech, and Signal Processing(ICASSP), pages 266–270, 2015.
D. Lee and H. Seung. Algorithms for non-negative matrix
factorization. In Neural Information Processing Systems(NIPS), pages 556–562, 2000.
T. G. Kolda and B. W. Bader. Tensor decompositions andapplications. SIAM Review, 51(3):455–500, 2009.
K. Tsuda, G. Rätsch, and M. K. Warmuth. Matrix exponentiated gradient updates for on-line learning and Bregman
projection. Journal of Machine Learning Research (JMLR),
6:995–1018, 2005.
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.
亀岡弘和. 非負値行列因子分解. 計測と制御, 51(9):835–844,
2012.
亀岡弘和. 非負値行列因子分解の音響信号処理への応用. 日本音響学会誌, 68(11):559–565, 2012.
H. Sawada, H. Kameoka, S. Araki, and N. Ueda. Eﬃcientalgorithms for multichannel extensions of Itakura-Saito nonnegative matrix factorization. In International Conferenceon Acoustics, Speech, and Signal Processing (ICASSP),
pages 261–264, 2012.
M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka. RWC
music database: Music genre database and musical instrument sound database. In International Conference on Music Information Retrieval (ISMIR), pages 229–230, 2003.
E. Vincent, R. Gribonval, and C. Févotte. Performancemeasurement in blind audio source separation. IEEE
Transactions on Acoustics, Speech, and Signal Processing,
14(4):1462–1469, 2006.
M. Hoﬀman, D. Blei, and P. Cook. Bayesian nonparametricmatrix factorization for recorded music. In International
Conference on Machine Learning (ICML), pages 439–446,
2010.
M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono,and S. Sagayama. Convergence-guaranteed multiplicativealgorithms for non-negative matrix factorization with beta
divergence. In International Workshop on Machine Learning for Signal Processing (MLSP), pages 283–288, 2010.

—8—