Vol.2016-MUS-112 No.21
2016/8/1

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

NMF vs PLCA: 多重音生成過程のための無限因子モデルと無限混合モデル
吉井 和佳1,a)

中村 栄太1,b)

糸山 克寿1,c)

後藤 真孝2,d)

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

行列分解においては，入力行列 X や推定すべき行列 A お

1. はじめに

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

行列分解 (matrix factorization) は，現在の機械学習やを定式化することにより，様々な変種を得ることができる．

データマイニング分野を支える基礎技術である．推薦シス

例えば，主成分分析 (principal component analysis: PCA)

M ×N

テムや音源分離などさまざまな応用で，行列 X ∈ R

では，A が固有ベクトルの集合，B がそれらの重みに対応が入力として与えられたときに，X ≈ AB となるようなしており，X と AB の二乗誤差が最小化されるように学

二つの行列 A ∈ R

習を行う．一方，独立成分分析 (independent component

M ×K

および B ∈ R

K×N

を求める問題が現れる．ここで重要なのは，K をもとの行列 X のサイ

analysis) では，A が互いに独立なベクトルの集合，B がズである M や N よりもずっと小さくしておけば，コンパそれらの重みに対応しており，やはり X と AB の二乗誤

クトな表現が得られることである．実際，X の自由度（独差が最小化されるように学習を行うのが一般的である．

立に調節可能な変数の個数）は M N であるが，AB の自

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

由度は K(M + N ) となっており，自由度が大幅に制限さ非負値に制約した非負値行列因子分解 (nonnegative matrix

れている．この結果，X に含まれる冗長性が取り除かれ，

factorization: NMF) や確率的潜在要素解析 (probabilistic

X を構成する N 個の列ベクトルに含まれる典型的な K 個

latent component analysis: PLCA) と呼ばれる行列分解技のパターン（基底ベクトル）が A に，それらがどの程度含法が広く用いられている [1–5]．過去の文献では，両者はまれているかの重みが B として得られる．

（数式の上では）等価であると記述されたり，NMF を確率

1

2

a)
b)
c)
d)

京都大学 大学院情報学研究科 知能情報学専攻
Yoshida-honmachi, Sakyo, Kyoto, Kyoto 606-8501, Japan
産業技術総合研究所 情報技術研究部門
1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan
yoshii(at)i.kyoto-u.ac.jp
enakamura(at)sap.ist.i.kyoto-u.ac.jp
itoyama(at)kuis.kyoto-u.ac.jp
m.goto(at)aist.go.jp

c 2016 Information Processing Society of Japan
⃝

モデルとして解釈したものが PLCA であると紹介される場合があるが，このような理解は正しくない．確率モデルに基づく統計的推論の見地に立てば，NMF は因子モデル(factor model)，PLCA は混合モデル (mixture model) に対応しており，観測スペクトログラムの生成過程の捉え方に決定的な違いがある．

1

Vol.2016-MUS-112 No.21
2016/8/1

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

一般に，観測データの背後にある基本パターン（クラNMFの生成モデル

PLCAの生成モデル

ス）を発見しようと思えば，因子モデルと混合モデルのうちで適切な方を選択する必要がある．混合モデルでは，観＋測データ中の各サンプルはひとつのクラスに排他的に属す＋

振幅/パワースペクトル

ヒストグラムることを仮定している．したがって，混合モデルを用いると，標準的なクラスタリングの問題を解くことができる．
一方，因子モデルは，各サンプルは複数の（場合によってたがって，因子モデルを用いると，多重クラスタリングの問題を解くことができる．音源分離問題においては，観測される混合音スペクトログラムに含まれる基本パターンは

カ
ウ
ン
ト

/

は全ての）クラスに同時に属することを許容している．し

振
幅
パ
ワ
ー

周波数

図1

周波数

NMF と PLCA における観測スペクトルの生成モデル．NMF
は確率変数（振幅/パワースペクトル）の和に基づいて定式化

音源に対応していることを期待したうえで，何を「サンプされており，PLCA は確率分布（ヒストグラム）の和に基づル」と解釈するかが重要なポイントである．

いて定式化されている．

因子モデルである NMF では，観測スペクトログラムの各フレームにおけるスペクトルをサンプルとみなし，それ(gamma process: GaP) あるいはベータ過程 (beta process:を複数の基底スペクトルの重み付き和で近似する（図 1）．

BP) を用いることにより，GaP-NMF および BP-NMF が基底スペクトルは，ある音源の平均的なスペクトルに対定式化できる．特に，BP-NMF においては，各フレーム応しており，各フレームは複数の音源に割り当てられるにおいて各音源が存在するかどうかの二値変数が導入さことになる．このとき，観測スペクトルが従う，基底ス

れており，音源区間の推定が同時に行える利点がある．一

ペクトルの重み付き和をパラメータに持つ確率分布を定方，PLCA に対しては，ディリクレ過程 (Dirichlet process,

めることで，Euclidean NMF（EU-NMF, ガウス分布に対DP) を用いることにより，DP-PLCA が定式化できる．

応）[6]，Kullback-Leiber NMF（KL-NMF, ポアソン分布に対応）[7]，Itakura-Saito NMF（IS-NMF, 複素ガウス分

2. 関連研究布に対応）[8] など様々な変種を定式化することができる．

NMF や PLCA は，音楽音響信号の音源分離や自動採譜

本稿では特に，最も広く利用されている点と，PLCA とのに最も利用されている技術であり，性能を向上させるため関連から，KL-NMF について取り上げる．KL-NMF は，のさまざまな改良が試みられている．本章では，それぞれポアソン分布に基づいているため，観測スペクトルの振幅の確率モデルに関する関連技術を紹介する．

は整数値に制限されている（実用上は実数で問題ない）．
一方，混合モデルである PLCA では，時間・周波数平面

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

上の観測スペクトログラムをヒストグラム（カウントデー

無限混合モデルは，無限個の混合比と要素分布をパラタ）とみなし，そのヒストグラムを構成するひとつひとつ

メータに持ち，それら生成するための事前分布としてディのサンプルをある音源に排他的に割り当てていく（図 1）．

リクレ過程 (DP) を用いる．Sethuraman [9] は，無限個の混合モデルでは，複数のスペクトルの重畳を本質的に表現

混合比を明示的に生成するため，DP の構成法のひとつでできない．苦肉の策として，観測スペクトルを各々が周波ある棒折り過程 (stick-breaking process: SBP) を導出して数をもつ音粒子の集合，すなわちヒストグラムであるとみいる．Blei ら [10] は，モデルパラメータの事後分布を近似なすことで，音粒子は一つの音源にしか属せないが，スペ的に計算するため，変分ベイズ法 (variational Bayes: VB)

クトル全体としてみれば複数の音源に属していると解釈に基づく反復推定を提案している．しかし，実際の計算機できる．これは，自然言語処理分野におけるトピックモデでは無限個のパラメータを取り扱うことはできないため，

ルの考え方に着想を得ており，ある文書に着目すると，そあらかじめ定めておいた要素数で計算を打ち切る必要がこに含まれる単語は一つのトピックに排他的に割り当てあった．すなわち，最初に十分に大きな要素数で VB を初られるが，文書としてみれば複数のトピックに属すると解期化し，反復計算が進むにつれて，不要とされる要素を削釈するのと同じである．PLCA という名前も，代表的なト

除していく方法が取られていた．

ピックモデルある確率的潜在意味解析 (probabilistic latentsemantic analysis: PLSA) に倣って名づけられている．

DP の別の構成法として，中華料理店過程 (Chinese restaurant process: CRP) が知られている．CRP では，無限個のこのような違いを念頭におくと，観測スペクトルに含ま混合比を事前分布である DP で積分消去することで，無限れる音源数を自動的に推定するため，NMF や PLCA に対個のパラメータを陽に取り扱う必要がなくなる．Neal [11]するノンパラメトリックベイズモデルを自然に定式化するは，モデルパラメータの事後分布を近似的に計算するたことができる．具体的には，NMF に対しては，ガンマ過程

め，ギブスサンプリングを導出している．この方法では，

c 2016 Information Processing Society of Japan⃝2

Vol.2016-MUS-112 No.21
2016/8/1

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

観測データを表現するのに必要な要素数が適宜増減可能なため，VB のように必要以上の個数の要素を考慮することによる計算量の無駄がない．

2.3 ベータ過程に基づく無限因子モデル
無限因子モデルを定式化する別の方法は，各基底の生起確率（コイントスで表が出る確率）を考え，無限個ある生

PLCA [12] は，基本的なトピックモデルである PLSA [13]起確率の事前分布としてベータ過程 (BP) を導入することの拡張であるため，PLSA のベイズ拡張である潜在的ディである．各フレームでは，その生起確率に従って，無限個

リクレ配分法 (latent Dirichlet allocation: LDA) [14] と同の基底の中から複数の基底を同時に選ぶことができるが，様のベイズ拡張が可能である．Teh ら [15] は，トピック数

生起確率は非常にスパース（ほとんどの基底の生起確率は（PLCA における基底数に相当）を自動的に調節するため，

ほぼゼロ）であるため，実際には限られた少数の基底しか階層ディリクレ過程 (hierarchical Dirichlet process: HDP)

生起しない．Teh ら [20] は，無限個の生起確率を明示的に基づく LDA を提案している．LDA では，異なる文書間に生成するため，BP に対する SBP を導出している．さでトピックを共有するため，一階層の DP ではなく HDP

らに，スライスサンプリングとギブスサンプリングを併用を用いる必要がある．もし，各文書のトピックの事前分布することで，打ち切り近似なしで，モデルパラメータの事として文書ごとに独立な DP を仮定すると，推定されるト後分布からサンプリングを可能にしている．Gupta ら [21]ピック（単語分布）は各文書に独自のものとなってしまう．

は，SBP 表現を用いた BP-NMF を提案している．Paisley

これを防ぐため，各文書の DP を束ねる上位階層の DP がら [22] は，新たな形式を持つ SBP を導出し，VB を用い導入されている．

た事後分布推論方法を提案している．

PLCA は，調波構造を表現する基底をあらかじめ準備しBP の別の構成方法として，インド料理過程 (Indian buf-

たり，学習しておくことによって，優れた自動採譜精度を

fet process: IBP) が知られている [23]．DP に対する CRP

達成している [16–18]．しかし，基底数を自動推定するノと同様に，無限個の生起確率を積分消去することで，BP

ンパラメトリックベイズ拡張についてはこれまで提案されに対する IBP が得られる．また，BP の簡便な構成法としていなかった．後の章で述べるように，PLCA の無限モデて，弱極限近似も可能であり，各基底の生起確率に対してルを構成する上では，HDP ではなく DP を用いれば十分

ベータ事前分布を導入し，その第一パラメータを非常に小であり，VB あるいはギブスサンプリングのいずれを用いさく設定すればよい．Liang ら [24, 25] は，この近似に基ても事後分布の近似計算が可能である．

づく BP-NMF を提案しており，VB を用いた事後分布推論方法が示されている．

2.2 ガンマ過程に基づく無限因子モデル
無限因子モデルを定式化する方法の一つは，無限個の基

3. ノンパラメトリックベイズ NMF

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

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

トルをスパースに誘導する（限られた少数の重みのみがゼ

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

ロではない実効的な値を持つ）ように，事前分布としてガンベイズ拡張について紹介する．

マ過程 (GaP) を導入することである．ここで，無限混合モ
デルにおける混合比とは異なり，基底の重みの総和は 1 と3.1 NMF に基づく音源分離なるように正規化されている必要はない．Roychowdhury

混合音に対して短時間フーリエ変換を行って得られる複素ら [19] は，無限個の重みを明示的に生成するため，GaP に

スペクトログラムを X̃ = [x̃1 , · · · , x̃N ] ∈ CM ×N ，k 番目の対する SBP を提案し，行列分解モデルに対する VB を導出している．しかし，この手法の実装は複雑であり，今だ

音源信号の複素スペクトログラムを X̃ k = [x̃k1 , · · · , x̃kN ] ∈CM ×N とする．M は周波数ビン数，N はフレーム数であほとんど利用されていない．一方，DP に対する CRP にる．観測した混合音が K 個の音源信号の周波数領域瞬時相当する表現はまだ確立されておらず，効率的なサンプリ

混合であると仮定すると，以下が成立する．

ング方法は研究の途上にある．

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

X̃ =

K
!

k=1

X̃ k

⇔ x̃n =

K
!

x̃kn

(1)

k=1

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

いま，観測変数 X̃ を潜在変数 X̃ k に分解したい．しかばよい．Hoﬀman ら [5] は，この近似に基づく GaP-NMF

し，これは不良設定問題であるので，X̃ k に対応する振幅スを提案しており，VB を用いることで，観測データに合わ×N(xknm =
ペクトログラム X k = [xk1 , · · · , xkN ] ∈ RM+

せて適切に基底数が調節できることを報告している．ただし，最初に十分多くの基底を準備しておき，VB の反復に

×N
|x̃knm |) は，ランク 1 の行列 Y k = [y k1 , · · · , y kN ] ∈ RM+

で近似できるものと仮定する（図 2）．

合わせて不要な基底を削除していく方式のため，計算量はX k ≈ w k hTk = Y kdef

比較的多くなることに注意が必要である．

c 2016 Information Processing Society of Japan
⃝

(2)

3

Vol.2016-MUS-112 No.21
2016/8/1

情報処理学会研究報告
IPSJ SIG Technical Report
観測 と再構成との距離を最小化するよう と を反復最適化

4.0

周
波
数

2.0

3.0

1.0
0

時間

図 2

2.0

図 3 形状母数が異なる幾つかのガンマ分布．

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

"K
ここで，振幅スペクトルの加法性，すなわち X = k=1 X k"Kおよび Y = k=1 Y k が成立することを仮定すると，
X≈K!

wk hTk = Y
def

(3)

k=1

×N
が得られる．ここで，X = [x1 , · · · , xN ] ∈ RM
Y =+×N
[y 1 , · · · , y N ] ∈ RM
とすると，
+

xn ≈

1.0

K
!

def

hkn wk = y n

(4)

k=1

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

xknm ∼ Poisson(yknm )(5)

ここで，複数の音源信号の重畳における振幅スペクトルの加法性を仮定すると（実際には成立しないことに注意），ポ

xnm ∼ Poisson(ynm )(6)を得る．ここで，尤度関数である式 (6) の対数をとって符号を反転させると，c

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

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

K
!

def

θk hkn wk = y n

(9)

k=1

ここで，θk ≥ 0 は基底 k の大域的な重みである．この θ に対し，観測データ X を表現するのに必要な基底 k 以外の要素 θk がゼロとなるようなスパースな学習を行いたい．
次に，θ ，W ，H に対して事前分布を導入する．まず，

W 及び H の各要素は非負値であるので，ガンマ事前分布を用いると都合が良い．
w
wkm ∼ Gamma(aw0 , b0 )hkn ∼ Gamma(ah0 , bh0 )

(10)
(11)

ここで，a∗0 > 0 及び b∗0 > 0 はそれぞれ，ガンマ分布の形

アソン分布の再生性から− log(xnm |ynm ) = xnm log

3.2 ガンマ過程に基づく NMF

xnm− xnm + ynm
ynm= DKL (xnm ||ynm )ても同様にガンマ事前分布を仮定する．
# αc $,αθk ∼ GammaK(12)

ここで，α > 0 及び c > 0 は超パラメータである．形状パラ(7)となり，負の対数尤度は xnm と ynm の KL ダイバージェ
ンスと定数項を除いて等しいことが分かる．したがって，
式 (6) の最大化（最尤推定）は式 (7) の最小化と等価である．KL-NMF の目標は，全時間・周波数平面にわたる対"数尤度関数の和 nm log p(xnm |ynm ) を最大化するようなW および H を求めることである．

NMF の実行結果を用いて，混合音の音源分離を行うに．
メータが小さくなるほど 0 が出る確率が大きくなる（図 3）"c
ただし，Eprior [θk ] = K ，Eprior [ k θk ] = c である．
ここで，式 (10)，式 (11) 及び式 (12) で構成される有限モデルに対して，K → ∞ となる極限を考えると，以下のガンマ過程が得られる．

G ∼ GaP(α, G0 )(13)N
ここで，G0 は空間 U（w ∈ RM
+ と h ∈ R+ の直積空間）上．
に定義された基底測度であり，G0 (U ) = c を満たす（図 4）
このとき，G は U 上の離散測度となり，空間 U の任意のは，ウィナーフィルタを用いる．

yknm
E[x̃knm |x̃nm ] =
x̃nm
ynm

状パラメータと逆尺度パラメータである．更に，θ に対し(8)

分割 {Ui }Ii=1 に対して

この処理では，X̃ k の位相は X̃ の位相と同一であるというG(Ui ) ∼ Gamma(αG0 (Ui ), α)

仮定がおかれている．最後に，逆フーリエ変換を用いて，が成立している．ただし，E[G] = G0 である．無限小区間

E[X̃ k |X̃] から k 番目の音源信号を復元すればよい．

への分割を {Uk }∞k=1 とすると，G(Uk ) = θk である．α は

c 2016 Information Processing Society of Japan
⃝

(14)

4

Vol.2016-MUS-112 No.21
2016/8/1

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

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

基底測度

図 4

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

集中度と呼ばれ，α が小さくなるほど θ はよりスパースになる．計算機上では K → ∞ は扱えないが，K を α に比べて十分大きな値に設定すれば，式 (12) はガンマ過程の良い近似となる（弱極限近似）．

3.2.2 変分ベイズ法

×N
Require: 非負値行列 X ∈ RM
, 最大基底数 K ，ガンマ過程の+wh h
集中度 α，ガンマ事前分布のパラメータ aw
0 , b0 , a0 , b0
1: 変分事後分布 q(θ), q(W ), q(H) をランダムに初期化
2: while not converged do
3:
λknm ∝ exp(log Eq [yknm ])!
!
+ nm λknm xnm , α + nm Eq [wkm hkn ])4:
q(θk ) = G( αcK
5:λknm ∝ exp(log Eq [yknm ])!
!
w
6:q(wkm ) = G(aw0 +n λknm xnm , b0 +n Eq [θk hkn ])7:
λknm ∝ exp(log Eq [yknm ])!
!
h8:q(hkn ) = G(ah0 +m λknm xnm , b0 +m Eq [θk wkm ])9: end while
10: return 変分事後分布 q(θ), q(W ), q(H)

式 (17) で与えられる変分下限 L(q) の第一項は式 (6) で式 (6), (10), (11), (12) で定義されるノンパラメトリック計算できる対数ポアソン尤度の期待値であるが，依然としベイズ KL-NMF (GaP-KL-NMF) に対する VB についてて解析的に計算できない．そのため，凹関数 log(x) に対し述べる [5]．今，観測データ x が与えられたときに，ベイて Jensen の不等式を用いると更なる変分下限

ズの定理を用いて未知パラメータ θ, W , H の事後分布p(θ, W , H|X) =p(X, θ, W , H)
p(X)

(15)

を計算したい．しかし，分母の周辺尤度 p(X) は解析的に計算できないため，対数周辺尤度 log p(X) の変分下限 L(q)を構成し，逐次最大化を行うことで log p(X) を近似した
い．具体的には，任意の分布 q(θ, W , H) を導入し，凹関
数 log(x) に対して Jensen の不等式を用いると%p(X, θ, W , H)dθdW dHlog p(X) = log q(θ, W , H)q(θ, W , H)%p(X, θ, W , H)
≥ q(θ, W , H) log
dθdW dHq(θ, W , H)

= Eq(θ ,W ,H ) [log p(X, θ, W , H)]def− Eq(θ ,W ,H ) [log q(θ, W , H)] = L(q) (16)を得る．等号成立条件は q(θ, W , H) = p(θ, W , H|X) であり，このとき L(q) が最大値をとる．しかし，真の事後分布 p(θ, W , H|X) は計算困難であるため，変分事後分布を因子分解可能な形 q(θ, W , H) = q(θ)q(W )q(H) に限定
し，その中でも以下の通り計算できる変分下限

nm

&

k

xnm Eq log

!

k

)

!
yknm
−
Eq [yknm ]
λknm
nm
k
knm
*
+ !
!
!
yknm
≥
−
xnm
λknm Eq log
Eq [yknm ]
λknm
nm
=

!

λknm

k

knm

def

= Eq [log q(X|θ, W , H)](21)"を得る．ここで，λknm は k λknm = 1 を満たす補助変数である．変分下限を最大化する λknm は，ラグランジュの未定乗数法を用いて，λknm ∝ exp(log Eq [yknm ]) と求まる．
最後に，各パラメータの変分事後分布を導出する．実際には式 (17) で与えられる元の変分下限 L(q) ではなく，式 (21)を用いて得られた更なる変分下限を最大化することになる．
すなわち，式 (18), (19), (20) において，log p(X, θ, W , H)の代わりに次式を用いればよい．

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

L(q) = Eq [log p(X|θ, W , H)] + Eq [log p(θ)](22)

具体的には，最適な変分事後分布 q ∗ (θ) は，θ に関連する+ Eq [log p(W )] + Eq [log p(H)]+ H(q(θ)) + H(q(W )) + H(q(H))

Eq [log p(X|θ, W , H)]
&
'
()
!
!
!
c
= Eq
xnm logyknm −
yknm(17)を最大化するものを求めたい．ここで，H(·) はエントロ
ピーを表す．これは，変分事後分布 q(θ)q(W )q(H) の真の事後分布 p(θ, W , H|X) に対する KL ダイバージェンス

項のみを取り出すと以下の通り計算できる．
!
!
clog q ∗ (θ) =xnm λknm log θk −
θk Eq [wkm ] Eq [hkn ]knm+! ## αc
k

K

$

knm

− 1 log θk − αθk

$

(23)

を最小化することと等価である．式 (17) を逐次最大化す

従って，θk の事後分布はガンマ分布となることが分かる．

るには，以下の更新式を収束するまで繰り返せば良い．

同様に，最適な q ∗ (W ) や q ∗ (H) もガンマ分布として求まq(θ) ∝ exp(Eq(H ,W ) [log p(X, θ, W , H)])(18)

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

c 2016 Information Processing Society of Japan⃝

る．Algorithm 1 に更新則を示す．反復ごとに，E[θk ] が十分に小さい基底 k を削除していけば，実効的な基底数

K+ を自動的に推定できる．最終的に，音源分離を行う際には，各パラメータの期待値を用いるのが一般的である．

5

Vol.2016-MUS-112 No.21
2016/8/1

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

Algorithm 2 GaP-KL-NMF に対するギブスサンプリング
×N
Require: 非負値行列 X ∈ RM
, 最大基底数 K ，ガンマ過程の+wh h
集中度 α，ガンマ事前分布のパラメータ aw
0 , b0 , a0 , b0

1: パラメータ θ, W , H をランダムに初期化
2: while not converged do
3:
λknm ∝ yknm!
!
4:θk ∼ G( αc+ nm λknm xnm , α + nm wkm hkn )K
5:λknm ∝ yknm!
!
w
6:
wkm ∼ G(aw0 +n λknm xnm , b0 +n θk hkn )7:
λknm ∝ yknm!
!
h8:
hkn ∼ G(ah0 +m λknm xnm , b0 +m θk wkm )9: end while
10: return パラメータ θ, W , H

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

3.3.1 ベイズモデルの定式化
ま ず ，式 (4) に 対 し ，K 次 元 の 二 値 ベ ク ト ル z n =[z1n , z2n , · · · , zKn ] を導入する．
xn ≈

K
!

def

znk hkn wk = y n

(28)

k=1

ここで，znk ∈ {0, 1} はフレーム n における基底 k のオン・オフを制御する二値変数である．すなわち，znk はア

3.2.3 ギブスサンプリング
決定的な近似推論方法である VB とは対称的に，確率的な近似推論法であるギブスサンプリングを用いれば，パラメータの事後分布からのサンプルを得ることができる．
式 (15) で与えられる真の事後分布は解析的に計算できないとしても，あるパラメータ（例えば θ ）の事後分布が，他のパラメータ（例えば H および W ）が既知であれば，解析的に計算できる場合がある．このとき，そのような「条件クティベーション hkn に対するマスクとして作用する．この Z ∈ {0, 1}N ×K に対し，観測データ X を表現するのに必要な基底 k 以外の変数 znk がゼロとなるようなスパース

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

znk ∼ Bernoulli(πk )(29)付き」事後分布からパラメータをサンプリングすることを次に，π ，W ，H に対して事前分布を導入する．まず，各パラメータについて順番に繰り返せば，真の事後分布か

W 及び H に対しては，GaP-KL-NMF と同様に，式 (10)らのサンプルが得られる．GaP-KL-NMF においては，以および式 (11) を用いる．コイン投げの確率 πk に対しては，下の手順を反復すればよい．

ベータ事前分布を仮定する．
,
αc α(K − c),πk ∼ Beta
KKθ ∼ p(θ|H, W , X)

(24)

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

(25)

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

具体的には，θ に着目すると，前節の VB と同様の変分下限を考えることにより

+

k

K

$

knm

− 1 log θk − αθkを得る．ここで，λknm は"$が小さくなるほど 0 が出る確率が大きくなる．ただし，"Eprior [znk ] = Kc ，Eprior [ k znk ] = c である．
ここで，式 (10)，式 (11) 及び式 (29) で構成される有限ベータ過程が得られる．

≥ log q(θ, H, W , X)!
!
c=xnm λknm log θk −
θk wkm hkn
knm

ここで，α > 0 及び c > 0 は超パラメータであり， αc
K

モデルに対して，K → ∞ となる極限を考えると，以下のlog p(θ|H, W , X) ∝ log p(θ, H, W , X)

! ## αc(30)

G ∼ BP(α, G0 )(31)N
ここで，G0 は空間 U （w ∈ RM
+ と h ∈ R+ の直積空間）(27)

k λknm = 1 を満たす補助変数である．変分下限を最大化する λknm は，ラグランジュの未定乗数法を用いて，λknm ∝ yknm と求まる．VB とは異上に定義された基底測度であり，G0 (U ) = c を満たす．このとき，G は U 上の離散測度となり，空間 U の無限小区間への分割 {Uk }∞k=1 に対して.
/
G(Uk ) ∼ Beta αG0 (Uk ), α(1 − G0 (Uk ))(32)

なり，λknm を最適化すれば常に等号を成立させることがが成立している．ただし，E[G] = G0 であり，G(Uk ) = πk可能であり，変分下限のタイトな評価が可能になっている．

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

最終的に，θk の条件付き事後分布はガンマ分布となるこ

いため，ディリクレ過程やガンマ過程と異なり，任意の分とが分かる．Algorithm 2 に更新則を示す．VB とは異割に対する周辺分布が陽に求まらないことに注意する．αなり，θk が非常に小さな値となった場合であっても，確率は集中度と呼ばれ，α が小さくなるほど θ はよりスパース

的に復活する可能性があるので，安易に削除してはいけなになる．計算機上では K → ∞ は扱えないが，K を α にい．音源分離を行う際には，各パラメータの条件付き事後比べて十分大きな値に設定すれば，式 (30) はベータ過程の分布の期待値を用いるのが一般的である．

良い近似となる（弱極限近似）．

c 2016 Information Processing Society of Japan⃝

6

Vol.2016-MUS-112 No.21
2016/8/1

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

3.3.2 変分ベイズ法

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

式 (6), (10), (11), (29), (30) で定義される BP-KL-NMF
では，通常の VB の適用は困難である．一見，GaP-KL-

NMF と同様に，Jensen の不等式を用いて対数尤度の変分下限が得られるように思われる．実際，計算を行うとlog p(X|π, Z, H, W )
(33)
!
c !
≥
xnm λknm log(znk wkm hkn ) −znk wkm hkn
knm

knm

が得られるが，znk は二値変数であるため，第一項の対数の中の値がゼロとなってしまう場合が存在する．したがっ×N
Require: 非負値行列 X ∈ RM
, 最大基底数 K ，ベータ過程の+w
h h
パラメータ α，c，ガンマ事前分布のパラメータ aw
0 , b0 , a0 , b0
1: パラメータ Z, W , H をランダムに初期化
2: while not converged do
3:
λknm ∝ yknm
4:
znk ∼ 式 (36) および (37)5:λknm ∝ yknm!
!
w
6:
wkm ∼ G(aw0 +n λknm xnm , b0 +n znk hkn )7:
λknm ∝ yknm!
!
h8:
hkn ∼ G(ah0 +m λknm xnm , b0 +m znk wkm )9: end while
10: return パラメータ Z, W , Hて，k, n, m に関する和に分解することができない．

4. ノンパラメトリックベイズ PLCA

3.3.3 ギブスサンプリング
本節では，ギブスサンプリングに基づく近似推論法について説明する [25]．本手法の肝は，二値行列 Z のサンプリングにある．具体的には，π ，Z ¬nk ，W , H が与えられたもとで（¬ はそのインデクス以外の変数の集合を表す），
あるフレーム n における基底 k の状態を表す二値変数 znk

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

= γ πk

4.1 PLCA に基づく分離
PLCA に お い て は ，混 合 音 の 振 幅 ス ペ ク ト ロ グラム

化しておく必要がある．経験的には，時間・周波数平面上における振幅の平均が 1 となるように混合音のスペクトロ¬k(ynm+ wkm hkn )xnm e−wkm hkn

(34)

m

グラム正規化し，四捨五入をしておくことで，PLCA が期
待通り動作することが多かった．ここで，量子化されたひp(znk = 0|π, Z ¬nk , W , H, X)
0¬kPoisson(xnm |ynm))
= γ(1 − πk )とつひとつの「サンプル」を xi (1 ≤ i ≤ I) で表すものとm= γ ′ (1 − πk )¬k
ここで，ynm=ついて提案する．

みなす．したがって，あらかじめ整数値になるように量子m0

ディリクレ過程に基づくノンパラメトリックベイズ拡張に

X = [x1 , · · · xn ] ∈ RM ×N をヒストグラムデータであるとp(znk = 1|π, Z ¬nk , W , H, X)
0¬kPoisson(xnm |ynm
+ wkm hkn )= γπk′

本章では，PLCA の音源分離への適用方法について述べ，0¬k xnm
(ynm
)

(35)

m

"

i̸=k wim hin であり，γ および γ′は正規化定数である．また，式 (29) 及び (30) には共役性が成立しているため，パラメータ π を積分消去することで，効率的な周辺化ギブスサンプリングを行うこともできる．

p(znk = 1|Z ¬nk , W , H, X)
(36)
⎛
⎞
!
αc ⎠ 0 ¬k∝⎝
zn′ k +
(ynm + wkm hkn )xnm e−wkm hkn
K
′
m
n ̸=n

p(znk = 0|Z ¬nk , W , H, X)
⎛
⎞
!
0
α(K
−
c)
¬k xnm⎠ (ynm∝⎝(1 − zn′ k ) +
)
K
′
m

(37)

n ̸=n

その他のパラメータの推論方法は GaP-KL-NMF と同様する．I はサンプル数であり，ヒストグラムの総和である"ので I = nm xnm となる．

xi は確率変数であり，その取りうる値は，フレームイ

ンデクスおよび周波数インデクスのペアである．具体的には，時間・周波数平面に大量の粒子がばらまかれ，積み重なってたくさんの山ができている場面を想像するとわかりやすい．ある粒子 xi が時間・周波数平面に投げ込まれたとき，その粒子がたまたまフレーム n・周波数 m で止まったとすると，xi = (n, m) となる．このとき，xi が従う時
間・周波数平面上の確率分布を p(n, m) とする．時間・周波数平面上のヒストグラム（混合音のスペクトログラム）
は，p(n, m) から得られたサンプル群である．

PLCA では，p(n, m) が因子分解できることを仮定する．
p(n, m) = p(n)K!

p(k|n)p(m|k)

(38)

k=1

すなわち，xi の生成過程では，まず，p(n) に従ってフレーであり，Algorithm 3 に更新則を示す．GaP-KL-NMF とム n が選択され，次に，p(k|n) に従って基底 k が選択さ同じく，πk が非常に小さな値となった場合であっても，確れ，最後に，p(m|k) に従って周波数 m が選択される．通

率的に復活する可能性があるので，安易に削除してはいけ常，p(n) は観測データから経験的に求めることができるのない．音源分離を行う際には，各パラメータの条件付き事

で，p(k|n) および p(m|k) を EM アルゴリズムを用いて推後分布の期待値を用いるのが一般的である．

定することが行われている [12]．

c 2016 Information Processing Society of Japan⃝7

Vol.2016-MUS-112 No.21
2016/8/1

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

確率分布が推定されれば，NMF と同じくウィナーフィ
ルタを用いて音源分離ができる．

E[x̃knm |x̃nm ] =p(n, m, k)x̃nmp(n, m)(39)

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

ここで，p(n, m, k) = p(n)p(k|n)p(m|k) である．

z i ∼ Categorical(π)

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

xi ∼ま ず ，観 測 変 数 X = {x1 , · · · , xI } お よ び 潜 在 変 数xz(θkm φkn ) inm ik(42)

n=1 m=1

分布を用いると都合がよい．

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

M ×N

として 1-of-M N 表現で表す．すなわち，サンプル i がフレーム n・周波数 m を選択する場合，xinm = 1 であり，xi
中のそれ以外の要素はすべて 0 である．また，各サンプルはある基底 k をひとつ選択するため，サンプル i に対応すφk ∼ Dir(γ)(44)とする．さらに，π に対しては，DP の構成法の一つである棒折り過程 (SBP) に基づく事前分布を仮定する．

πk = vkk−1
0(1 − vk′ )

(45)

k′ =1

中のそれ以外の要素はすべて 0 である．
を考慮して，以下の等価な表現に変換しておく．

(43)

ディリクレ分布について，全ての次元で同じ値を取るもの

サンプル i が基底 k を選択する場合，zik = 1 であり，z i

さらに，式 (38) における確率変数 n, m, k の間の依存性θ k ∼ Dir(β)

ここで，β および γ は超パラメータであり，それぞれのる潜在変数 z i ∈ {0, 1}K を 1-of-K 表現で表す．すなわち，

p(k)p(n|k)p(m|k)(41)

次に，パラメータ π, θ, φ に対する事前分布を導入する．

Z = {z 1 , · · · , z I } を定義しておく．ヒストグラム中の各サ

p(n, m) =M
0まず， θ, φ に対しては，共役事前分布であるディリクレ

4.2.1 ベイズモデルの定式化

K
!

N
0

vk ∼ Beta(1, α)(46)

ここで，α は超パラメータであり，π が v に変数変換されていることに注意する．あるいは，より簡便にディリクレ

(40)

k=1

これは，k について対称な形となっており，p(k) を混合
比，p(n|k)p(m|k) を要素分布とみなせば，p(n, m) が混合
分布として定式化されていることが分かる．したがって，分布を導入する方法もある．

π ∼ Dir

#α$

(47)
K
ここで，式 (43)，式 (44) 及び式 (47) で構成される有限モデルに対して，K → ∞ となる極限を考えると，以下のディリクレ過程が得られる．

K → ∞ として無限混合分布を考える際に，ディリクレ過G ∼ DP(α, G0 )(48)

程を利用することができる．一方，式 (38) においては，フ"K
レーム n ごとに混合分布 p(m|n) = k=1 p(k|n)p(m|k) をN
ここで，G0 は空間 U（θ ∈ RM
+ と φ ∈ R+ の直積空間）上

HDP を用いる必要が出てきてしまい，複雑になる．

とき，G は U 上の離散測度となり，空間 U の任意の分割考えているので，無限化する際には，LDA と同じ理由でLDA には HDP が必要であるのに，PLCA では DP で十分である理由は，PLCA では PLSA と同様にフレームの生起確率 p(n) を考慮しているからである．この結果，式 (38)から式 (40) への変換が可能になっている．しかし，p(n)に定義された確率測度であり，G0 (U ) = 1 を満たす．この{Ui }Ii=1 に対して[G(U1 ), G(U2 ), · · · , G(UI )]∼ Dir(αG0 (U1 ), αG0 (U2 ), · · · , αG0 (UI ))(49)

はすでに観測済みのフレーム N 個に対する確率分布であが成立している．ただし，E[G] = G0 である．無限小区間

り，未知の観測データ（新たなフレーム）の生成機構を有への分割を {Uk }∞k=1 とすると，G(Uk ) = πk である．α はしておらず，「生成モデル」としては完全とは言えない．一

集中度と呼ばれ，α が小さくなるほど π はよりスパースに方，LDA では，p(n) を考慮しないことにより，新たに観なる．計算機上では K → ∞ は扱えないが，K を α に比測される文書に対する汎化能力を獲得している点に注意がべて十分大きな値に設定すれば，式 (47) は DP の良い近似必要である．PLCA のベイズ化においても同様のアプローとなる（弱極限近似）．一般には，SBP に基づく事前分布

チが可能であるが，予備実験においては正しく音源分離をである式 (45) を用いた方が近似誤差が小さく，良い結果を行うことができなかった．

与えることが多い．あとで触れるように，無限次元の π をいま，p(k), p(n|k), p(m|k) はすべて標準的な離散分布で

c 2016 Information Processing Society of Japan⃝

積分消去した CRP 表現に基づく推論も可能である．

8

Vol.2016-MUS-112 No.21
2016/8/1

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

4.2.2 変分ベイズ法

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

DP-PLCA においても，真の事後分布 p(Z, v, θ, φ|X) を求めることは困難であるが，VB を用いて近似的に事後分布を推定することができる．具体的には，因子分解できる形に限定した変分事後分布 q(Z, v, θ, φ) = q(Z)q(v)q(θ)q(φ)
を考え，真の事後分布との KL ダイバージェンスを最小化するように反復最適化を行う．

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

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

Require: 観測データ X = {x1 , · · · , xI }, 最大基底数 K ，ディリ
クレ過程のパラメータ α，ディリクレ事前分布のパラメータ β, γ1: 変分事後分布 q(Z), q(v), q(θ), q(φ) をランダムに初期化
2: while not converged do
3:q(z i ) = Categorical(η
i)
"
#
!
! !
′
4:
q(vk ) = Beta 1 + i E[zik ], α + i K
k′ =k+1 E[zik ]5:q(θ k ) = Dir(ω k )6:
q(φk ) = Dir(λk )7: end while
8: return 変分事後分布 q(Z), q(v), q(θ), q(φ)

機上での取り扱いに都合がよい．また，ディリクレ事前分布の共役性から，他のパラメータである θ および φ も同時

これに従って計算を進めると，Z の変分事後分布は以下のに積分消去することができ，推論すべき確率変数は Z のみ

通り求まる．

になる．この結果，収束が早く，局所解にも陥りにくい性
質のよいサンプラーが構成できる．推論時には，観測デー

q(z i ) = Categorical(η i )(54)

ここで，ηik = !Kρik ρ ′ であり，ρik は次式で求まる．
k′ =1

ik

log ρik = E[log vk ] +
+

!

k−1
!

k′ =1

xinm E[log θkm ] +nmE[log(1 − vk′ )]!

ら適切な個数を探索することができる．

5. おわりに
本稿では，音楽音響信号に対する音源分離のための主要xinm E[log φkn ]

(55)

nm

な二つの行列分解技法である NMF と PLCA について，対応する確率モデルの性質を明らかにし，比較検討を行った．
具体的には，ガンマ過程およびベータ過程に基づく NMF

次に，v の変分事後分布は以下の通り求まる．

q(vk )

タを表現するのに必要な基底数 K を自動的に増減しながのノンパラメトリックベイズモデルを紹介し，ディリクレ
過程に基づく PLCA を提案した．今後は大規模なデータ'= Beta 1 +!

E[zik ], α +

i

K
! !
i

(

E[zik′ ]k′ =k+1

で比較実験を行っていきたい．

(56)

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

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

26280089，16H01744, JST CREST OngaCREST，および栢
森情報科学振興財団の支援を受けた．

q(θk ) = Dir(ω k )

(57)

参考文献

q(φk ) = Dir(λk )

(58)

[1]
[2]

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

ωkm = β +!

xinm E[zik ](59)xinm E[zik ]

(60)

n

λkn = β +

!

[3]

m

Algorithm 4 に更新則を示す．GaP-KL-NMF と同じ
く，反復ごとに，E[πk ] が十分に小さい基底 k を削除して[4]

いけば，実効的な基底数 K+ を自動的に推定できる．最終的に，音源分離を行う際には，各パラメータの期待値を用[5]

いるのが一般的である．

4.2.3 ギブスサンプリング
DP-PLCA に対しては，中華料理店過程 (CRP) に基づ[6]く効率的な周辺化ギブスサンプリングを行うことができる．この手法では，無限個の混合比である π が積分消去されており，明示的に実体化する必要がなくなるため，計算

c 2016 Information Processing Society of Japan⃝[7]

亀岡弘和. 非負値行列因子分解の音響信号処理への応用.
日本音響学会誌, 68(11):559–565, 2012.
P. Smaragdis, C. Févotte, G. Mysore, N. Mohammadiha,and M. Hoﬀman. Dynamic source separation using nonnegative factorizations: A unified view. IEEE Signal
Processing Magazine, 31(3):66–75, 2014.
H. Kameoka and K. Kashino. Composite autoregressive
system for sparse source-filter representation of speech.
In IEEE International Symposium on Circuits and Systems (ISCAS), pages 2477–2480, 2009.
K. Yoshii and M. Goto. Infinite composite autoregressive models for music signal analysis. In International
Conference on Music Information Retrieval (ISMIR),
pages 79–84, 2012.
M. Hoﬀman, D. Blei, and P. Cook. Bayesian nonparametric matrix factorization for recorded music. In International Conference on Machine Learning (ICML),
pages 439–446, 2010.
D. Lee and H. Seung. Algorithms for non-negative matrix
factorization. In Neural Information Processing Systems (NIPS), pages 556–562, 2000.
P. Smaragdis and J. C. Brown. Non-negative matrix factorization for polyphonic music transcription. In IEEE

9

Vol.2016-MUS-112 No.21
2016/8/1

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

[8]

[9]
[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA), pages 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.
J. Sethuraman. A constructive definition of Dirichlet priors. Statistica Sinica, 4:639–650, 1994.
D. M. Blei and M. I. Jordan. Variational inference for
Dirichlet process mixtures. Bayesian Analysis, 1(1):121–
144, 2006.
R. M. Neal. Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational
and Graphical Statistics, 9(2):249–265, 2000.
M. Shashanka, B. Raj, and P. Smaragdis. Probabilistic latent variable models as nonnegative factorizations.
Computational Intelligence and Neuroscience, 2008:1–
8, 2008.
T. Hofmann. Learning the similarity of documents: An
information-geometric approach to document retrieval
and categorization. In Advances in Neural Information
Processing Systems (NIPS), pages 914–920, 2000.
D. M. Blei, A. Y. Ng, and M. I. Jordan. Latent Dirichlet allocation. Journal of Machine Learning Research,
3:993–1022, 2003.
Y. W. Teh, M. I. Jordan, M. J. Beal, and D. M. Blei. Hierarchical Dirichlet processes. Journal of the American
Statistical Association, 101:1566–1581, 2006.
G. Grindlay and D. P. W. Ellis. Transcribing multiinstrument polyphonic music with hierarchical eigeninstruments. Journal Selected Topics Signal Processing,
5(6):1159–1169, 2011.
B. Fuentes, R. Badeau, and G. Richard. Harmonic adaptive latent component analysis of audio and application to music transcription. IEEE Transactions on Audio, Speech, and Language Processing, 21(9):1854–1866,
2013.
E. Benetos and S. Dixon. Multiple-instrument polyphonic music transcription using a temporally constrained shift-invariant model. Journal of Acoustical Society of America, 133(3):1727–1741, 2013.
A. Roychowdhury and B. Kulis. Gamma processes,
stick-breaking, and variational inference. In International Conference on Artificial Intelligence and Statistics (AISTATS), pages 800–808, 2015.
Y. W. Teh, D. Görür, and Z. Ghahramani. Stickbreaking construction for the Indian buﬀet process. In
International Conference on Artificial Intelligence andStatistics (AISTATS), pages 556–563, 2007.
S. K. Gupta, D. Phung, and S. Venkatesh. A nonparametric Bayesian Poisson gamma model for countdata. In International Conference on Pattern Recognition (ICPR), pages 1815–1818, 2012.
J. Paisley, L. Carin, and D. Blei. Variational inference
for stick-breaking beta process priors. In International
Conference on Machine Learning (ICML), pages 889–
896, 2011.
T. L. Griﬃths and Z. Ghahramani. Infinite latent featuremodels and the Indian buﬀet process. In Advances in
Neural Information Processing Systems (NIPS), pages
475–482, 2006.
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.

c 2016 Information Processing Society of Japan⃝[25]

D. Liang and M. D. Hoﬀman. Beta process non-negative
matrix factorization with stochastic structured meanfield variational inference. In NIPS Workshop on Advances in Variational Inference, pages 1–6, 2014.

10