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

信学技報
TECHNICAL REPORT OF IEICE.

［ポスター講演］音楽音響信号解析のための
ディリクレ過程に基づくベイズ潜在成分分析
吉井 和佳†

中村 栄太†

糸山 克寿†

後藤 真孝††† 京都大学 大学院情報学研究科 〒 606–8501 京都府京都市左京区吉田本町†† 産業技術総合研究所 情報技術研究部門 〒 305–8568 茨城県つくば市梅園 1–1–1
E-mail: †{yoshii,enakamura,itoyama}@sap.ist.i.kyoto-u.ac.jp, ††m.goto@aist.go.jp
あらまし 本稿では，音楽音響信号に対する音源分離を目的とした確率的潜在成分解析 (probabilistic latent compo-

nent analysis: PLCA) のノンパラメトリックベイズ拡張について述べる．最近よく利用される非負値行列因子分解(nonnegative matrix factorization: NMF) では，各フレームにおける混合音のスペクトルを，少数の基底スペクトルの重み付き和で近似する．すなわち，各フレームにおいて複数の音源が同時に生起することが許容されており，NMF
は因子モデルの一種である．一方，PLCA では，時間・周波数平面上のスペクトログラムをヒストグラムであるとみなし，その背後にある確率分布を推定する．このとき，各時間・周波数ビンの振幅値を量子化し，仮想的な音粒子の観測個数であるとみなしたうえで，各粒子をいずれかの音源に排他的に割り当てるため，PLCA は混合モデルの一種である．これまで，NMF の方が物理的に自然な解釈ができるにもかかわらず，実際には PLCA も盛んに利用され，優れた分離結果が得られている．本稿では，因子モデルである NMF に対してはガンマ過程あるいはベータ過程を，混合モデルである PLCA に対してはディリクレ過程を用いることにより，基底数を自動調節するためのノンパラメトリッ
クベイズモデルを構成できることを示し，変分ベイズ法あるいはギブスサンプリングを用いた推論方法を導出する．
キーワード 音楽音響信号解析，確率的潜在成分解析，非負値行列分解，音源分離，ノンパラメトリックベイズ

音楽音響信号の音源分離においては，X, A, B の要素を非

1. は じ め に

負値に制約した非負値行列因子分解 (nonnegative matrix fac-

行列分解 (matrix factorization) は，現在の機械学習やデー

torization: NMF) や確率的潜在成分解析 (probabilistic latent

タマイニング分野を支える基礎技術である．推薦システムやcomponent analysis: PLCA) と呼ばれる行列分解技法が広く

音源分離などさまざまな応用で，行列 X ∈ RM ×N が入力と用いられている [1, 2]．過去の文献では，両者は（数式の上でして与えられたときに，X ≈ AB となるような二つの行列は）等価であると記述されたり，NMF を確率モデルとして解

A∈R

M ×K

および B ∈ R

K×N

を求める問題が現れる．ここ釈したものが PLCA であると紹介される場合があるが，こので重要なのは，K をもとの行列 X のサイズである M や N よような理解は正しくない．確率モデルに基づく統計的推論の見りもずっと小さくしておけば，コンパクトな表現が得られるこ地に立てば，NMF は因子モデル (factor model)，PLCA は混とである．実際，X の自由度（独立に調節可能な変数の個数）

合モデル (mixture model) に対応しており，観測スペクトログ

は M N であるが，AB の自由度は K(M + N ) となっており，

ラムの生成過程の捉え方に決定的な違いがある．

自由度が大幅に制限されている．この結果，X に含まれる冗長

一般に，観測データの背後にある基本パターン（クラス）を性が取り除かれ，X を構成する N 個の列ベクトルに含まれる発見しようと思えば，因子モデルと混合モデルのうちで適切な

典型的な K 個のパターン（基底ベクトル）が A に，それらが方を選択する必要がある．混合モデルでは，観測データ中の各どの程度含まれているかの重みが B として得られる．

サンプルはひとつのクラスに排他的に属することを仮定してい行列分解において，入力行列 X や推定すべき行列 A およびる．したがって，混合モデルを用いると，標準的なクラスタリ

B が満たすべき性質を定め，X と AB との近似誤差を定式化

ングの問題を解くことができる．一方，因子モデルは，各サン

すれば，様々な変種が得られる．例えば，主成分分析 (principal

プルは複数の（場合によっては全ての）クラスに同時に属する

component analysis: PCA) では，A が固有ベクトルの集合，ことを許容している．したがって，因子モデルを用いると，多B がそれらの重みに対応しており，X と AB の二乗誤差が最小

重クラスタリングの問題を解くことができる．音源分離問題に化されるように学習を行う．一方，独立成分分析 (independent

おいては，観測される混合音スペクトログラムに含まれる基本

component analysis) では，A が互いに独立なベクトルの集合，パターンは音源に対応していることを期待したうえで，何をB がそれらの重みに対応しており，やはり X と AB の二乗誤「サンプル」と解釈するかが重要なポイントである．

差が最小化されるように学習を行う方法が提案されている．

因子モデルである NMF では，観測スペクトログラムの各フ—1—

レームにおけるスペクトルをサンプルとみなし，それを複数のNMFの生成モデル

PLCAの生成モデル．基底スペクト
基底スペクトルの重み付き和で近似する（図 1）
ルは，ある音源の平均的なスペクトルに対応しており，各フレー＋

ムは複数の音源に割り当てられることになる．このとき，観測＋

振幅/パワースペクトル

ヒストグラムスペクトルが従う，基底スペクトルの重み付き和をパラメータに持つ確率分布を定めることで，Euclidean NMF（EU-NMF,
カウント

アソン分布に対応）[4]，Itakura-Saito NMF（IS-NMF, 複素ガ

振幅 パ
/ ワー

ガウス分布に対応）[3]，Kullback-Leiber NMF（KL-NMF, ポ
ウス分布に対応）[5] など様々な変種を定式化することができ

周波数

周波数

る．本稿では特に，最も広く利用されている点と，PLCA との関連から，KL-NMF について取り上げる．KL-NMF は，ポア
ソン分布に基づいているため，観測スペクトルの振幅は整数値に制限されている（実用上は実数で問題ない）．

図 1 観測スペクトルの生成モデル．NMF は確率変数（振幅/パワースペクトル）の和に基づいて定式化されており，PLCA は確率分布（ヒストグラム）の和に基づいて定式化されている．

一方，混合モデルである PLCA では，時間・周波数平面上の観測スペクトログラムをヒストグラム（カウントデータ）と程 (stick-breaking process: SBP) を導出している．Blei ら [7]みなし，そのヒストグラムを構成するひとつひとつのサンプル

は，モデルパラメータの事後分布を近似的に計算するため，変．しかし，混合モをある音源に排他的に割り当てていく（図 1）

分ベイズ法 (variational Bayes: VB) に基づく反復推定を提案

デルでは，複数のスペクトルの重畳を本質的に表現できない．

している．しかし，実際の計算機では無限個のパラメータを取苦肉の策として，観測スペクトルを各々が周波数をもつ音粒子

り扱うことはできないため，あらかじめ定めておいた要素数での集合，すなわちヒストグラムであるとみなすことで，音粒子計算を打ち切る必要があった．すなわち，最初に十分に大きなは一つの音源にしか属せないが，スペクトル全体でみれば複数要素数で VB を初期化し，反復計算が進むにつれて，不要とさの音源に属していると解釈できる．これは，自然言語処理分野

れる要素を削除していく方法が取られていた．

におけるトピックモデルの考え方に着想を得ており，ある文書

DP の別の構成法として，中華料理店過程 (Chinese restau-に着目すると，そこに含まれる単語は一つのトピックに排他的

rant process: CRP) が知られている．CRP では，無限個の混に割り当てられるが，文書としてみれば複数のトピックに属す

合比を事前分布である DP で積分消去することで，無限個のパると解釈するのと同じである．PLCA という名前も，代表的なラメータを陽に取り扱う必要がなくなる．Neal [8] は，モデル

トピックモデルある確率的潜在意味解析 (probabilistic latent

パラメータの事後分布を近似的に計算するため，ギブスサンプ

semantic analysis: PLSA) に倣って名づけられている．

リングを導出している．この方法では，観測データを表現する

このような違いを念頭におくと，観測スペクトルに含まれるのに必要な要素数が適宜増減可能なため，VB のように必要以音源数を自動的に推定するため，NMF や PLCA に対するノン上の個数の要素を考慮することによる計算量の無駄がない．

パラメトリックベイズモデルを自然に定式化することができる．

PLCA [9] は，基本的なトピックモデルである PLSA [10] の

具体的には，NMF に対しては，ガンマ過程 (gamma process:拡張であるため，PLSA のベイズ拡張である潜在的ディリクレ

GaP) あるいはベータ過程 (beta process: BP) を用いること配分法 (latent Dirichlet allocation: LDA) [11] と同様のベイにより，GaP-NMF および BP-NMF が定式化できることが既ズ拡張が可能である．Teh ら [12] は，トピック数（PLCA におに知られている．特に，BP-NMF においては，各フレームにおける基底数に相当）を自動的に調節するため，階層ディリクレ

いて各音源が存在するかどうかの二値変数が導入されており，過程 (hierarchical Dirichlet process: HDP) に基づく LDA を

音源区間の推定が同時に行える利点がある．一方，PLCA に対提案している．LDA では，異なる文書間でトピックを共有す

しては，ディリクレ過程 (Dirichlet process: DP) を用いるこるため，一階層の DP ではなく HDP を用いる必要がある．もとにより，DP-PLCA が定式化できる．

し，各文書のトピックの事前分布として文書ごとに独立な DP

2. 関 連 研 究を仮定すると，推定されるトピック（単語分布）は各文書に独NMF や PLCA は，音楽音響信号の音源分離や自動採譜に最ねる上位階層の DP が導入されている．

自のものとなってしまう．これを防ぐため，各文書の DP を束も利用されている技術であり，性能を向上させるためのさまざ

PLCA は，調波構造を表現する基底をあらかじめ準備した

まな改良が試みられている．本章では，それぞれの確率モデル

り，学習しておくことによって，優れた自動採譜精度を達成しに関する関連技術を紹介する．

ている [13, 14]．しかし，基底数を自動推定するノンパラメト

2. 1 ディリクレ過程に基づく無限混合モデル

リックベイズ拡張についてはこれまで提案されていなかった．

無限混合モデルは，無限個の混合比と要素分布をパラメー

後の章で述べるように，PLCA の無限モデルを構成する上でタに持ち，それら生成するための事前分布としてディリクレ

は，HDP ではなく DP を用いれば十分であり，VB あるいはギ

過程 (DP) を用いる．Sethuraman [6] は，無限個の混合比を

ブスサンプリングのいずれを用いても事後分布の近似計算が可

明示的に生成するため，DP の構成法のひとつである棒折り過

能である．

—2—

2. 2 ガンマ過程に基づく無限因子モデル
無限因子モデルを定式化する方法の一つは，無限個の基底にパースに誘導する（限られた少数の重みのみがゼロではない実周波数

対する非負の重みを導入し，その無限次元の重みベクトルをス

観測 と再構成との距離を最小化するよう と を反復最適化

効的な値を持つ）ように，事前分布としてガンマ過程 (GaP) を導入することである．ここで，無限混合モデルにおける混合比
時間とは異なり，基底の重みの総和は 1 となるように正規化されている必要はない．Roychowdhury ら [15] は，無限個の重みを

図 2 混合音スペクトログラムに対する NMF．

明示的に生成するため，GaP に対する SBP を提案し，行列分解モデルに対する VB を導出している．しかし，この手法の実装は複雑であり，今だほとんど利用されていない．一方，DPに対する CRP に相当する表現はまだ確立されておらず，効率的なサンプリング方法は研究の途上にある．

GaP の簡便な構成法には，弱極限近似 (weak-limit approximation) が知られており，各重みに対してガンマ事前分布を

3. 1 NMF に基づく音源分離
混合音に対して短時間フーリエ変換を行って得られる振幅ス
×N，音源 k の振幅
ペクトログラムを X = [x1 , · · · , xN ] ∈ RM
+
M ×N
とする．
スペクトログラムを X k = [xk1 , · · · , xkN ] ∈ R+

M は周波数ビン数，N はフレーム数である．便宜的に，振幅スペクトルの加法性を仮定すると，以下が成立する．

導入し，その形状パラメータを非常に小さく設定すればよい．

Hoﬀman ら [2] は，この近似に基づく GaP-NMF を提案しておxn =り，VB を用いることで，観測データに合わせて適切に基底数

K


xkn

(1)

k=1

が調節できることを報告している．ただし，最初に十分多くの

ここで，K は音源数である．いま，観測変数 xn を潜在変数

基底を準備しておき，VB の反復に合わせて不要な基底を削除

xkn に分解したい．しかし，これは不良設定問題であるので，していく方式のため，計算量は比較的多くなることに注意が必

潜在変数 xkn は，音源 k に対応する基底スペクトル wk のス

要である．

ケール hkn の相似形 y kn で近似できると仮定する．

2. 3 ベータ過程に基づく無限因子モデル
無限因子モデルを定式化する別の方法は，各基底の生起確率（コイントスで表が出る確率）を考え，無限個ある生起確率の事前分布としてベータ過程 (BP) を導入することである．各
フレームでは，その生起確率に従って，無限個の基底の中から複数の基底を同時に選ぶことができるが，生起確率は非常にスパース（ほとんどの基底の生起確率はほぼゼロ）であるた
め，実際には限られた少数の基底しか生起しない．Teh ら [16]は，無限個の生起確率を明示的に生成するため，BP に対する

SBP を導出している．さらに，スライスサンプリングとギブ
スサンプリングを併用することで，打ち切り近似なしで，モデ
ルパラメータの事後分布からサンプリングを可能にしている．

Gupta ら [17] は，SBP 表現を用いた BP-NMF を提案していdef

xkn ≈ hkn wk = y kn

(2)

したがって，式 (2) を式 (1) へ代入すると，

xn ≈

K


def

hkn wk = y n

(3)

k=1

を得る．すなわち，混合音のスペクトル xn は，K 個の基底ス
ペクトル wk の線形和 y n で近似できることを意味する．
音源分離によく利用される KL-NMF に対応する確率モデル
では，潜在変数 xknm が yknm をパラメータに持つポアソン分
布（E[xknm ] = yknm ）に従うことを仮定する．

xknm ∼ Poisson(yknm )(4)

る．Paisley ら [18] は，新たな形式を持つ SBP を導出し，VBを用いた事後分布推論方法を提案している．

BP の別の構成方法として，インド料理過程 (Indian buﬀetprocess: IBP) が知られている [19]．DP に対する CRP と同

ここで，ポアソン分布の再生性から

xnm ∼ Poisson(ynm )(5)様に，無限個の生起確率を積分消去することで，BP に対するを得る．ここで，尤度関数である式 (5) の対数をとって符号をIBP が得られる．また，BP の簡便な構成法として，弱極限近反転させると，似も可能であり，各基底の生起確率に対してベータ事前分布を導入し，その第一パラメータを非常に小さく設定すればよい．

Liang ら [20, 21] は，この近似に基づく BP-NMF を提案して

c− log(xnm |ynm ) = xnm log

xnm− xnm + ynm
ynm= DKL (xnm ||ynm )(6)おり，VB を用いた事後分布推論方法が示されている．

3. ノンパラメトリックベイズ NMFとなり，負の対数尤度は xnm と ynm の KL ダイバージェンスと定数項を除いて等しい．したがって，式 (5) の最大化（最尤推定）

ンマ過程およびベータ過程に基づくノンパラメトリックベイズ

は式 (6) の最小化と等価である．KL-NMF の目標は，全時間・
周波数平面にわたる対数尤度関数の和 nm log p(xnm |ynm ) を

拡張について紹介する．

最大化するような W および H を求めることである．

本章では，NMF の音源分離への適用方法について述べ，ガ—3—

3. 2 ガンマ過程に基づく NMF

基底測度

本節では，GaP-KL-NMF について説明する．まず，GaP の
弱極限近似に基づくベイズモデルの定式化を説明し，VB およびギブスサンプリングの適用について述べる．

3. 2. 1 ベイズモデルの定式化
ま ず ，式 (3) に 対 し ，K 次 元 の 非 負 値 ベ ク ト ル θ =[θ1 , θ2 , · · · , θK ] を導入する．
xn ≈

K


def

θk hkn wk = y n

図 3 NMF のためのガンマ過程事前分布．

(7)を計算したい．しかし，分母の周辺尤度 p(X) は解析的に計算できないため，対数周辺尤度 log p(X) の変分下限 L(q) を構成

k=1

ここで，θk  0 は基底 k の大域的な重みである．この θ に対し，逐次最大化を行うことで log p(X) を近似したい．具体的

し，観測データ X を表現するのに必要な基底 k 以外の要素 θkには，任意の分布 q(θ, W , H) を導入し，凹関数 log(x) に対しがゼロとなるようなスパースな学習を行いたい．

w
wkm ∼ Gamma(aw0 , b0 )(8)て Jensen の不等式を用いると
p(X, θ, W , H)
log p(X) = log q(θ, W , H)dθdW dHq(θ, W , H)p(X, θ, W , H)dθdW dH q(θ, W , H) log
q(θ, W , H)

hkn ∼ Gamma(ah0 , bh0 )(9)= Eq(θ,W ,H) [log p(X, θ, W , H)]

次に，θ ，W ，H に対して事前分布を導入する．まず，W 及
び H の各要素は非負値であるので，ガンマ事前分布を用いる．

def

ここで，a∗0 > 0 及び b∗0 > 0 はそれぞれ，ガンマ分布の形状パラメータと逆尺度パラメータである．更に，θ に対しても同様にガンマ事前分布を仮定する．
 αc 
,αθk ∼ GammaK− Eq(θ,W ,H) [log q(θ, W , H)] = L(q)(14)

を得る．等号成立条件は q(θ, W , H) = p(θ, W , H|X) であり，このとき L(q) が最大値をとる．しかし，真の事後分布(10)

ここで，α > 0 及び c > 0 は超パラメータである．形状パラメータが小さくなるほど 0 が出る確率が大きくなる．ただし，
cEprior [θk ] = K，Eprior [ k θk ] = c である．

p(θ, W , H|X) は計算困難であるため，変分事後分布を因子分解可能な形 q(θ, W , H) = q(θ)q(W )q(H) に限定し，その中で L(q) を最大化するものを求めたい．具体的には，以下の更新式を収束するまで繰り返せば良い．

ここで，式 (8)，式 (9) 及び式 (10) で構成される有限モデル

q(θ) ∝ exp(Eq(H,W ) [log p(X, θ, W , H)])

(15)に対して，K → ∞ となる極限を考えると，以下のガンマ過程

q(H) ∝ exp(Eq(θ,W ) [log p(X, θ, W , H)])(16)が得られる．

q(W ) ∝ exp(Eq(θ,H) [log p(X, θ, W , H)])(17)

G ∼ GaP(α, G0 )(11)

変分下限 L(q) の計算では，対数ポアソン尤度の期待値が必要

N
ここで，G0 は空間 U（w ∈ RM
+ と h ∈ R+ の直積空間）上に定となるが，依然として解析的に計算できない．そのため，凹関．このと義された基底測度であり，G0 (U ) = c を満たす（図 3）

数 log(x) に対して Jensen の不等式を用いると更なる変分下限

き，G は U 上の離散測度となり，空間 U の任意の分割 {Ui }Ii=1
に対してG(Ui ) ∼ Gamma(αG0 (Ui ), α)(12)が成立している．ただし，E[G] = G0 である．無限小区間への分割を {Uk }∞k=1 とすると，G(Uk ) = θk である．α は集中度と呼ばれ，α が小さくなるほど θ はよりスパースになる．計算機
上では K → ∞ は扱えないが，K を α に比べて十分大きな値に設定すれば，式 (10) はガンマ過程の良い近似となる（弱極限Eq [log p(X|θ, W , H)]





c
= Eq
yknm −
yknm
xnm lognm=

xnm Eq log

nm




nm

k



xnm





k

λknm

k

λknm Eq log

k

yknm
λknm

−



Eq [yknm ]

knm


yknm
−
Eq [yknm ]
λknm
knm

def

= Eq [log q(X|θ, W , H)]

近似）．



(18)

3. 2. 2 変分ベイズ法を得る．ここで，λknm は式 (5), (8), (9), (10) で定義されるノンパラメトリックベイズ

ある．変分下限を最大化する λknm は，ラグランジュの未定乗

KL-NMF (GaP-KL-NMF) に対する VB について述べる [2]．
今，観測データ X が与えられたときに，ベイズの定理を用い

p(X, θ, W , H)p(X)

数法を用いて，λknm ∝ exp(log Eq [yknm ]) と求まる．
最後に，各パラメータの変分事後分布を導出する．実際には，

L(q) ではなく，式 (18) を用いて得られた更なる変分下限を最て未知パラメータ θ, W , H の事後分布p(θ, W , H|X) =k λknm = 1 を満たす補助変数で大化することになる．すなわち，式 (15), (16), (17) において，(13)

log p(X, θ, W , H) の代わりに次式を用いればよい．
—4—

Algorithm 1 GaP-KL-NMF に対する変分ベイズ法

Algorithm 2 GaP-KL-NMF に対するギブスサンプリング

×N
Require: 非負値行列 X ∈ RM
, 最大基底数 K ，ガンマ過程の集+
w h h
中度 α，ガンマ事前分布のパラメータ aw
0 , b0 , a0 , b0

×N
Require: 非負値行列 X ∈ RM
, 最大基底数 K ，ガンマ過程の集+

1: 変分事後分布 q(θ), q(W ), q(H) をランダムに初期化

1: パラメータ θ, W , H をランダムに初期化

2: while not converged do

2: while not converged do

λknm ∝ exp(log Eq [yknm ])q(θk ) = G( αc+ nm λknm xnm , α + nm Eq [wkm hkn ])
K

3:
4:
5:
6:
7:
8:

w h h
中度 α，ガンマ事前分布のパラメータ aw
0 , b0 , a0 , b0

3:
4:

λknm ∝ yknmθk ∼ G( αc+ nm λknm xnm , α + nm wkm hkn )K

λknm ∝ exp(log Eq [yknm ])wq(wkm ) = G(aw0 +n λknm xnm , b0 +n Eq [θk hkn ])5:

λknm ∝ yknm

6:

wkm ∼ G(aw
0 +λknm ∝ exp(log Eq [yknm ])hq(hkn ) = G(ah0 +m λknm xnm , b0 +m Eq [θk wkm ])7:

λknm ∝ yknm

8:

hkn ∼ G(ah
0 +





w
n λknm xnm , b0 +hm λknm xnm , b0 +9: end while

9: end while

10: return 変分事後分布 q(θ), q(W ), q(H)10: return パラメータ θ, W , Hn θk hkn )m θk wkm )3. 3 ベータ過程に基づく NMF

log q(X, θ, W , H) = log q(X|θ, W , H) + log p(θ)+ log p(W ) + log p(H)



(19)

Algorithm 1 に更新則を示す．反復ごとに，E[θk ] が十分に

本節では，BP-KL-NMF について説明する．まず，BP の弱極限近似に基づくベイズモデルの定式化を説明し，VB および
ギブスサンプリングの適用について述べる．

小さい基底 k を削除していけば，実効的な基底数 K+ を自動的

3. 3. 1 ベイズモデルの定式化に推定できる．最終的に，音源分離を行う際には，各パラメー

ま ず ，式 (3) に 対 し ，K 次 元 の 二 値 ベ ク ト ル z n =[zn1 , zn2 , · · · , znK ] を導入する．

タの期待値を用いるのが一般的である．

3. 2. 3 ギブスサンプリング
式 (13) で与えられる真の事後分布は解析的に計算できないが，あるパラメータ（例えば θ ）の事後分布が，他のパラメータ（例えば H および W ）が既知であれば，解析的に計算できる場合がある．このとき，そのような「条件付き」事後分布か
らパラメータをサンプリングすることを各パラメータについて

xn ≈

K


def

znk hkn wk = y n

(24)

k=1

ここで，znk ∈ {0, 1} はフレーム n における基底 k のオン・オフを制御する二値変数である．すなわち，znk はアクティベーショ
ン hkn に対するマスクとして作用する．この Z ∈ {0, 1}N ×K

順番に繰り返せば，真の事後分布からのサンプルが得られる．

に対し，観測データ X を表現するのに必要な基底 k 以外の変GaP-KL-NMF においては，以下の手順を反復すればよい．

数 znk がゼロとなるようなスパースな学習を行うため，znk は
非常に裏が出やすいコイン投げで確率的に定まるものとする．

θ ∼ p(θ|H, W , X)

(20)

H ∼ p(H|θ, W , X)

(21)

W ∼ p(W |θ, H, X)(22)を考えることによりlog p(θ|H, W , X) ∝ log p(θ, H, W , X)

+

k

K

knm



− 1 log θk − αθkを得る．ここで，λknm は次に，π ，W ，H に対して事前分布を導入する．まず，Wc，
なるほど 0 が出る確率が大きくなる．ただし，Eprior [znk ] = K
Eprior [ k znk ] = c である．

(23)k λknm = 1 を満たす補助変数である．変分下限を最大化する λknm は，ラグランジュの未定乗数法を用いて，λknm ∝ yknm と求まる．VB とは異なり，λknmを最適化すれば常に等号を成立させることが可能であり，変分下限のタイトな評価が可能になっている．

Algorithm 2 に更新則を示す．VB とは異なり，θk が非常に小さな値となった場合であっても，確率的に復活する可能性があるので，安易に削除してはいけない．

げの確率 πk に対しては，ベータ事前分布を仮定する．
αc α(K − c),(26)
πk ∼ BetaK
K
ここで，α > 0 及び c > 0 は超パラメータであり，αcが小さくK log q(θ, H, W , X)c=xnm λknm log θk −
θk wkm hkn
  αc(25)及び H に対しては，式 (8) および式 (9) を用いる．コイン投

具体的には，θ に着目すると，前節の VB と同様の変分下限

knm

znk ∼ Bernoulli(πk )ここで，式 (8)，式 (9) 及び式 (25) で構成される有限モデルに対して，K → ∞ となる極限を考えると，以下のベータ過程が得られる．

G ∼ BP(α, G0 )(27)

N
ここで，G0 は空間 U （w ∈ RM
+ と h ∈ R+ の直積空間）上に定義された基底測度であり，G0 (U ) = c を満たす．このとき，G は U 上の離散測度となり，空間 U の無限小区間への分
割 {Uk }∞
k=1 に対して


G(Uk ) ∼ Beta αG0 (Uk ), α(1 − G0 (Uk ))

(28)
—5—

が成立している．ただし，E[G] = G0 であり，G(Uk ) = πk でAlgorithm 3 BP-KL-NMF に対するギブスサンプリング

ある．ベータ過程はコルモゴロフの拡張定理を満たさないため，

×N
Require: 非負値行列 X ∈ RM
, 最大基底数 K ，ベータ過程のパ+w h h
ラメータ α，c，ガンマ事前分布のパラメータ aw
0 , b0 , a0 , b0

ディリクレ過程やガンマ過程と異なり，任意の分割に対する周辺分布が陽に求まらない．α は集中度と呼ばれ，α が小さくな

1: パラメータ Z, W , H をランダムに初期化るほど θ はよりスパースになる．計算機上では K → ∞ は扱

2: while not converged do
3:

λknm ∝ yknm

4:

znk ∼ 式 (32) および (33)5:λknm ∝ yknm

6:

式 (5), (8), (9), (25), (26) で定義される BP-KL-NMF では，

wkm ∼ G(aw0 +7:

λknm ∝ yknm

通常の VB は困難である．一見，Jensen の不等式を用いて対8:hkn ∼ G(ah0 +数尤度の変分下限が得られるように思われる．実際には，

9: end while

えないが，K を α に比べて十分大きな値に設定すれば，式 (26)はベータ過程の良い近似となる（弱極限近似）．

3. 3. 2 変分ベイズ法wn λknm xnm , b0 +hm λknm xnm , b0 +

n znk hkn )m znk wkm )10: return パラメータ Z, W , H

log p(X|π, Z, H, W )
(29)
c 


xnm λknm log(znk wkm hkn ) −znk wkm hkn
knm

knm

4. ノンパラメトリックベイズ PLCAが得られるが，znk は二値変数であるため，第一項の対数の中本章では，PLCA の音源分離への適用方法と，ディリクレ過の値がゼロとなってしまう場合が存在する．したがって，k, n,

程に基づくノンパラメトリックベイズ拡張について提案する．

4. 1 PLCA に基づく分離

m に関する和に分解することができない．

PLCA においては，混合音の振幅スペクトログラム X =3. 3. 3 ギブスサンプリング
ギブスサンプリングにおいては，二値行列 Z の取り扱いが重M ×Nをヒストグラムデータであるとみなす．
[x1 , · · · xn ] ∈ R+

要である [21]．具体的には，π ，Z ¬nk ，W , H が与えられたもしたがって，あらかじめ整数値になるように量子化しておく必，あるフレーム nとで（¬ はそのインデクス以外の変数の集合）

要がある．経験的には，時間・周波数平面上における振幅の平の基底 k の状態を表す二値変数 znk をサンプリングする際に，均が 1 となるように混合音のスペクトログラム正規化し，四捨

基底 k がオンの場合とオフの場合とに分けて，尤度を計算する．

五入をしておくことで，PLCA が期待通り動作することが多かった．ここで，量子化されたひとつひとつの「サンプル」をp(znk = 1|π, Z ¬nk , W , H, X)¬kPoisson(xnm |ynm
+ wkm hkn )= γπkxi (1  i  I) で表すものとする．I はサンプル数であり，ヒnm xnm となる．

ストグラムの総和であるので I =m= γ πk ¬k(ynm + wkm hkn )xnm e−wkm hkn

xi は確率変数であり，その取りうる値は，フレームインデク(30)m

周波数平面に大量の粒子がばらまかれ，積み重なってたくさんp(znk = 0|π, Z ¬nk , W , H, X)¬kPoisson(xnm |ynm))
= γ(1 − πk )の山ができている場面を想像するとわかりやすい．ある粒子 xi
が時間・周波数平面に投げ込まれたとき，その粒子がたまたまフm= γ (1 − πk ) ¬k x(ynm ) nm

(31)

m
¬k
ここで，ynm=i=k| wim hin であり，γ および γ

ため，パラメータ π を積分消去することで，効率的な周辺化ギ
ブスサンプリングを行うこともできる．

このとき，xi が従う時間・周波数平面上の確率分布を p(n, m)
トログラム）は，p(n, m) から得られたサンプル群である．

PLCA では，p(n, m) が因子分解できることを仮定する．
p(n, m) = p(n)p(znk = 1|Z ¬nk , W , H, X)(32)
αc  ¬k∝zn k +
(ynm + wkm hkn )xnm e−wkm hkn
K

m
n =n
|

n =n
|

レーム n・周波数 m で止まったとすると，xi = (n, m) となる．とする．時間・周波数平面上のヒストグラム（混合音のスペクは正規化定数である．また，式 (25) 及び (26) には共役性が成立している

p(znk = 0|Z ¬nk , W , H, X)α(K − c)∝(1 − zn k ) +Kスおよび周波数インデクスのペアである．具体的には，時間・Kp(k|n)p(m|k)

(34)

k=1

すなわち，xi の生成過程では，まず，p(n) に従ってフレーム

n が選択され，次に，p(k|n) に従って基底 k が選択され，最後に，p(m|k) に従って周波数 m が選択される．通常，p(n) は観測データから経験的に求めることができるので，p(k|n) および ¬k x(ynm ) nmp(m|k) を EM アルゴリズムを用いて推定することが行われて(33)m

Algorithm 3 に更新則を示す．GaP-KL-NMF と同じく，
πk が非常に小さな値となった場合であっても，確率的に復活する可能性があるので，安易に削除してはいけない．

いる [9]．

4. 2 ディリクレ過程に基づく PLCA
本節では，DP-PLCA について説明する．まず，ベイズモデ
ルを定式化し，VB および周辺化ギブスサンプリングの適用について述べる．

—6—

4. 2. 1 ベイズモデルの定式化

Algorithm 4 DP-PLCA に対する変分ベイズ法

まず，観測変数 X = {x1 , · · · , xI } および潜在変数 Z =Require: 観測データ X = {x1 , · · · , xI }, 最大基底数 K ，ディリク{z 1 , · · · , z I } を定義する．ヒストグラム中の各サンプルは，離レ過程のパラメータ α，ディリクレ事前分布のパラメータ β, γ

散的なフレームインデクスと周波数インデクスを取るため，サンプル i の値を xi ∈ {0, 1}

M ×N

として 1-of-M N 表現で表す．

すなわち，サンプル i がフレーム n・周波数 m を選択する場合，

xinm = 1 であり，xi 中のそれ以外の要素はすべて 0 である．
また，各サンプルはある基底 k をひとつ選択するため，サンプ
ル i に対応する潜在変数 z i ∈ {0, 1}K を 1-of-K 表現で表す．
さらに，式 (34) における確率変数 n, m, k の間の依存性を

1: 変分事後分布 q(Z), q(v), q(θ), q(φ) をランダムに初期化
2: while not converged do
4:q(z i ) = Categorical(η
i)



 
q(vk ) = Beta 1 + i E[zik ], α + i Kk =k+1 E[zik ]5:q(φk ) = Dir(λk )6:q(θ k ) = Dir(ω k )3:

7: end while
8: return 変分事後分布 q(Z), q(v), q(θ), q(φ)

考慮して，以下の等価な表現に変換しておく．

p(n, m) =Kp(k)p(n|k)p(m|k)

(35)

k=1

いは，より簡便にディリクレ分布を導入する方法もある．

π ∼ Dir

これは，k について対称な形となっており，p(k) を混合比，p(n|k)p(m|k) を要素分布とみなせば，p(n, m) が混合分布として定式化されていることが分かる．したがって，K → ∞として無限混合分布を考える際に，ディリクレ過程を利用で

α

(40)

K

ここで，式 (38) 及び式 (40) で構成される有限モデルに対し，

K → ∞ となる極限を考えると，以下のディリクレ過程を得る．
G ∼ DP(α, G0 )(41)

きる．一方，式 (34) においては，フレーム n ごとに混合分布p(m|n) = K
k=1 p(k|n)p(m|k) を考えているので，無限化する

N
ここで，G0 は空間 U （θ ∈ RM
+ と φ ∈ R+ の直積空間）上に際には，LDA と同じ理由で HDP を用いる必要がある．

定義された確率測度であり，G0 (U ) = 1 を満たす．このとき，PLCA では HDP ではなく DP で十分である理由は，PLSAと同様にフレームの生起確率 p(n) を考慮しているからである．
この結果，式 (34) から式 (35) への変換が可能になる．しかし，p(n) はすでに観測済みのフレーム N 個に対する確率分布であり，未知の観測データ（新たなフレーム）の生成機構を有してG は U 上の離散的な確率測度となり，空間 U の任意の分割{Ui }Ii=1 に対して[G(U1 ), G(U2 ), · · · , G(UI )]∼ Dir(αG0 (U1 ), αG0 (U2 ), · · · , αG0 (UI ))

(42)

おらず，
「生成モデル」としては完全ではない．一方，LDA でが成立する．ただし，E[G] = G0 である．無限小区間への分割

は，p(n) を考慮しないことで，新たに観測される文書に対するを {Uk }∞
k=1 とすると，G(Uk ) = πk である．α は集中度と呼ば汎化能力を獲得している．PLCA のベイズ化においても同様のアプローチが可能かもしれない．

K を α に比べて十分大きな値に設定すれば，式 (40) は DP の

いま，p(k), p(n|k), p(m|k) はすべて離散分布であるとし，そのパラメータを p(k) = πk ，p(n|k) = φkn ，p(m|k) = θkm と
する．ただし， k πk = 1， n φkn = 1， m θkm = 1 を満たす．これらを用いると，PLCA のベイズモデルを定式化できる．まず，X と Z に対する尤度関数は以下の通り与えられる．

z i ∼ Categorical(π)
xi ∼

M
N 


れ，α が小さくなるほど π はよりスパースになる．実際には，良い近似となる（弱極限近似）．一般には，SBP に基づく事前分布である式 (39) を用いた方が近似誤差が小さく，良い結果を与えることが多い．あとで触れるように，無限次元の π を積分消去した CRP 表現に基づく推論も可能である．

4. 2. 2 変分ベイズ法(36)

DP-PLCA においても，真の事後分布 p(Z, v, φ, θ|X) を解(37)

後分布を推定できる．具体的には，因子分解できる形に限定し析的に計算することは困難であるが，VB を用いて近似的に事(φkn θkm )xinm zik

n=1 m=1

た変分事後分布 q(Z, v, φ, θ) = q(Z)q(v)q(φ)q(θ) を考え（潜次に，パラメータ π, φ, θ に対する事前分布を導入する．まず，

在変数とパラメータの独立性 q(Z, v, φ, θ) = q(Z)q(v, φ, θ) をφ, θ に対しては，共役事前分布のディリクレ分布を用いる．

仮定するだけで自動的に導かれる），真の事後分布との KL ダφk ∼ Dir(β)θ k ∼ Dir(γ)(38)

イバージェンスを最小化するように反復最適化を行う．

q(Z) ∝ exp(Eq(v,φ,θ) [log p(X, Z, v, φ, θ)])(43)ては，DP の構成法の一つである棒折り過程 (SBP) に基づく事

q(v) ∝ exp(Eq(Z,φ,θ) [log p(X, Z, v, φ, θ)])(44)

前分布を仮定する．

q(φ) ∝ exp(Eq(Z,v,θ) [log p(X, Z, v, φ, θ)])(45)

q(θ) ∝ exp(Eq(Z,v,φ) [log p(X, Z, v, φ, θ)])(46)

ここで，β および γ は超パラメータである．さらに，π に対し

πk = v kk−1
(1 − vk )vk ∼ Beta(1, α)

(39)

k =1

ここで，α は超パラメータで，π は v の変数変換である．ある

これらの更新式に従って計算を進めると，Z の変分事後分布は
以下の通り求まる．

—7—

Algorithm 5 DP-PLCA に対するギブスサンプリング
Require: 観測データ X = {x1 , · · · , xI }, ディリクレ過程のパラ
メータ α，ディリクレ事前分布のパラメータ β, γ

本稿では，音楽音響信号に対する音源分離のための主要な二

1: パラメータ Z をランダムに初期化

つの行列分解技法である NMF と PLCA について，対応する

2: while not converged do

確率モデルの性質を明らかにし，比較検討を行った．具体的に

p(z i |Z ¬i , X) ∼ 式 (53) および (54)3:

5. お わ り には，ガンマ過程およびベータ過程に基づく NMF のノンパラ

4: end while

メトリックベイズモデルを紹介し，ディリクレ過程に基づく

5: return パラメータ Z

PLCA を提案した．今後は大規模なデータで比較実験を行うとともに，より効率的で局所解に対して頑健な変分ベイズ法およq(z i ) = Categorical(η i )ここで，ηik = K ρik ρ
k =1

ik

log ρik = E[log vk ] +
+



(47)

であり，ρik は次式で求まる．

k−1
xinm E[log φkn ] +16H01744, JST CREST OngaCREST，および栢森情報科学振興財

E[log(1 − vk )]

文xinm E[log θkm ]

(48)

nm

次に，v の変分事後分布は次式で求まる．
⎛
⎞
K

 
q(vk ) = Beta ⎝1 +E[zik ], α +E[zik ]⎠ (49)iik =k+1

最後に，φ および θ の変分事後分布は以下の通り求まる．

q(φk ) = Dir(λk )q(θ k ) = Dir(ω k )(50)

ここで，λk および ω k は以下で求められる．

λkn = β +

謝辞: 本研究の一部は，JSPS 科研費 24220006, 26700020, 26280089，

団の支援を受けた．

k =1

nm

びギブスサンプリングの開発に取り組みたい．

xinm E[zik ] ωkm = β +mxinm E[zik ] (51)nAlgorithm 4 に更新則を示す．GaP-KL-NMF と同じく，
E[πk ] が十分に小さい基底 k を削除していけばよい．
4. 2. 3 ギブスサンプリング
中華料理店過程 (CRP) に基づく効率的な周辺化ギブスサン．無限次元の π が積分
プリングも可能である（Algorithm 5）
消去されており，有限しか扱えない計算機上での取り扱いに都合がよい．また，ディリクレ事前分布の共役性から，θ および

φ も同時に積分消去することができ，推論すべき確率変数は Zのみになる．具体的には，潜在変数 z i の条件付き事後分布はp(z i |Z ¬i , X) ∝ p(z i , xi |Z ¬i , X ¬i )(52)

で与えられる．したがって，潜在変数 z i が既存の基底 k に割り当てられ（zik = 1），観測変数 xi がフレーム n および周波数 m に割り当てられる確率は（xinm = 1），p(zik = 1, xinm = 1|Z ¬i , X ¬i )
(53)

 
 





i =i
| zi km i =i
| xi nm zi k + βn i =i
| xi nm zi k + γ∝
I −1+αi =i| zi k + βNi =i
| zi k + γMで与えられ，z i が新たな基底 knew に割り当てられる確率は，

p(ziknew = 1, xinm = 1|Z ¬i , X ¬i ) ∝α1
I − 1 + α MN

(54)となる．このとき，実行的な音源数 K+ が増減しうる．

献[1] 亀岡弘和. 非負値行列因子分解の音響信号処理への応用. 日本音響学会誌, 68(11):559–565, 2012.
[2] M. Hoﬀman et al. Bayesian nonparametric matrix factorization for recorded music. ICML, 2010.
[3] D. Lee and H. Seung. Algorithms for non-negative matrix
factorization. NIPS, 2000.
[4] P. Smaragdis and J. Brown. Non-negative matrix factorization for polyphonic music transcription. WASPAA, 2003.
[5] C. Févotte et al. Nonnegative matrix factorization with the
Itakura-Saito divergence: With application to music analysis. Neural Computation, 21(3):793–830, 2009.
[6] J. Sethuraman. A constructive deﬁnition of Dirichlet priors.
Statistica Sinica, 4:639–650, 1994.
[7] D. Blei and M. Jordan. Variational inference for Dirichlet
process mixtures. Bayesian Analysis, 1(1):121–144, 2006.
[8] R. Neal. Markov chain sampling methods for Dirichlet process mixture models. J. of Com. and Grap. Stat., 2000.
[9] M. Shashanka et al. Probabilistic latent variable models as
nonnegative factorizations. Computational Intelligence and
Neuroscience, 2008:1–8, 2008.
[10] T. Hofmann. Learning the similarity of documents: An
information-geometric approach to document retrieval andcategorization. NIPS, 2000.
[11] D. M. Blei et al. Latent Dirichlet allocation. JMLR, 2003.
[12] Y. W. Teh et al. Hierarchical Dirichlet processes. JASA,
101:1566–1581, 2006.
[13] B. Fuentes et al. Harmonic adaptive latent component analysis of audio and application to music transcription. IEEETASLP, 21(9):1854–1866, 2013.
[14] E. Benetos and S. Dixon. Multiple-instrument polyphonicmusic transcription using a temporally constrained shiftinvariant model. JASA, 133(3):1727–1741, 2013.
[15] A. Roychowdhury and B. Kulis. Gamma processes, stickbreaking, and variational inference. AISTATS, 2015.
[16] Y. W. Teh et al. Stick-breaking construction for the Indianbuﬀet process. AISTATS, 2007.
[17] S. Gupta et al. A nonparametric Bayesian Poisson gamma
model for count data. ICPR, 2012.
[18] J. Paisley et al. Variational inference for stick-breaking beta
process priors. ICML, 2011.
[19] T. L. Griﬃths and Z. Ghahramani. Inﬁnite latent featuremodels and the Indian buﬀet process. NIPS, 2006.
[20] D. Liang et al. Beta process sparse nonnegative matrix factorization for music. ISMIR, 2013.
[21] D. Liang and M. D. Hoﬀman. Beta process non-negative
matrix factorization with stochastic structured mean-ﬁeld
variational inference. In NIPS Workshop on Advances in
Variational Inference, 2014.

—8—