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

信学技報
TECHNICAL REPORT OF IEICE.

非負値行列因子分解とその音響信号処理応用
亀岡弘和†† 東京大学 大学院情報理工学系研究科 / 日本電信電話株式会社 コミュニケーション科学基礎研究所〒 113-8656 東京都文京区本郷 7-3-1 / 〒 243-0198 神奈川県厚木市森の里若宮 3-1
E-mail: †kameoka@hil.t.u-tokyo.ac.jp / kameoka.hirokazu@lab.ntt.co.jp
あらまし

実世界にはパワースペクトル，画素値，頻度など，非負値で表されるデータが多い。このように非負値で表されたデータを有意な加法的な構成成分に分解することを目的とした多変量解析手法を非負値行列因子分解(non-negative matrix factorization: NMF) といい，近年音響信号処理分野をはじめ様々な分野で注目を集めている。
本稿では，NMF の基本性質，アルゴリズムの導出方法，生成モデルとしての解釈，音響信号処理への応用とそのための拡張モデルについて解説する。
キーワード

多変量解析，非負値データ，行列分解，スパース性，補助関数法，生成モデル，音響信号処理

Non-negative matrix factorization and its applications toaudio signal processing
Hirokazu KAMEOKA†† Graduate School of Information Science and Technology, The University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan
E-mail: †kameoka@hil.t.u-tokyo.ac.jp / kameoka.hirokazu@lab.ntt.co.jp
Abstract In this paper, I will give a brief introduction to a data analysis technique called non-negative matrix
factorization (NMF), which has attracted a lot of attention in the ﬁeld of audio signal processing in recent years.
I will mention some basic properties of NMF, eﬀects induced by the non-negative constraints, how to derive an
iterative algorithm for NMF, and some attempts that have been made to apply NMF to audio processing problems.
Key words Multivariate analysis, non-negative data, matrix factorization, sparsity, auxiliary function method,
generative model, audio signal processing

1. は じ め に
実世界には，パワースペクトル，画素値，頻度など，非負値で表されるデータが多い。主成分分析や独立成分分析などの多といい，様々な分野で近年注目を集めている。本稿では，NMFの定式化，基本的な性質，アルゴリズムの導出方法，音響信号処理の問題に焦点を当てた改良・拡張のアイディアについて解説する。

変量解析では，所与のデータを複数の加法的な成分に分解することを目的とするが，これと同様に上述のような非負値のデー

2. NMF とはタから構成成分を抽出することが役立つ場面が多い。例えば，

以後，データをベクトルで表記することとする。例えば，複数の音源の音響信号が混在する多重音のパワースペクトルか

画像データであれば各ピクセルの画素値がベクトル要素とら個々の音源のパワースペクトルをうまく取り出すことができなり，パワースペクトルであれば各周波数におけるパワー値

れば，雑音除去や音源分離などに役立てることができるし，顔がベクトル要素となる。今，N 個の非負値データベクトル

画像データを目や鼻などの顔のパーツに該当する画像データに=0,K が与えられたとしよう。これらを観測ベ
y 1 , . . . , y N ∈ R>

うまく分解することができれば，顔認証や顔画像合成などに役

クトルと呼ぶ。ただし，R=0,K は K 次元の非負値ベクトル全立てることができる。以上のように，非負値のデータを加法的体の集合を表す。NMF では，これら観測ベクトルが，いずれな構成成分に分解することを目的とした多変量解析手法を非負

も M 個の基底ベクトルの適当な重みつき和によって表された値行列因子分解 (non-negative matrix factorization: NMF) [1]

ものと見なされ，すべての観測ベクトルを最も良く説明するよ>—1—

基底行列と係数行列が重要な意味をもつ。

NMF では，観測データの中で共起する成分をひとまとめにしたものが基底ベクトルの推定結果になる傾向がある。例えば，
どのデータも成分 a，b，c から成っているとき，a，b，c を基底にすれば全データを完全に表現することは可能である。ここ
でもし，成分 a や成分 b が生起するときいつも成分 c も揃って生起するようなら，成分 a と成分 c，および成分 b と成分 c をひとまとめにしたものを基底としても当該データを同様に表現可能である。この方が，少ない基底でデータを表現可能な分，
より「節約的」である。このように，共起する成分をひとまとめにして基底と置いた方が，基底を節約でき，節約した分を他の成分にフィットするのに充てられるようになるわけである。
図 1 NMF によるスペクトログラムの分解

先に述べたように，NMF では少ない基底で観測データをでき

うな M 個の基底ベクトルおよび重み係数を推定することが目るだけ良く表す必要があるため，上記のような基底が自然に得的である。よって，NMF では加法性が成り立つ量のみが対象られる傾向にある。

となる。画素値やパワースペクトルはいずれも厳密には加法性

NMF では，係数行列の非負性の制約により，基底同士を加が成り立つ量ではないが，NMF を適用する場面では近似的に算することしかできない (減算ができない) ため，係数行列の要加法性が成り立つと仮定されているということに注意が必要で素がスパース (疎) になる傾向がある。これはベクトルの合成をある。パワースペクトルの非加法性については 6. で詳しく述べ

考えるとイメージしやすい。観測ベクトルを近似するには，観る。以上の加法性に関する仮定の他に，基底ベクトルと重み係

測ベクトルと近い方向を向いた基底ベクトル以外の係数はでき数がいずれも非負値であるという仮定が置かれる点も NMF にるだけ小さくした方が良い。スパース性は統計的独立性と概念おける重要なポイントである。すなわち，NMF は，観測ベク

的に関連が深く，上述の共起する成分をひとまとめにしたもの>
トル y n を基底ベクトル h1 , . . . , hM ∈ R=0,K の非負結合 (結を基底にする働きは，なるべく各基底の係数が互いにコヒーレ
ントにならないように基底を決定しようとする働きと見なすこ合係数 u1,n , . . . , uM,n が非負値の線形結合)とができる。つまり，NMF では非負制約による副次的な効果

yn ≃

M
X

hm um,n (n = 1, . . . , N )

(1)

m=1

で近似する問題と見なせる。ここで，観測ベクトルを並べた行
列を Y = [y 1 , . . . , y N ] = (yk,n )K×N ，基底ベクトルを並べた
行列を H = [h1 , . . . , hM ] = (hk,m )K×M ，結合係数 um,n を

m 行 n 列の要素とした行列を U = (um,n )M ×N とすると，式(1) は Y ≃ HU と同じ意味である。このように NMF は，観測ベクトルを並べた行列を二つの非負値行列の積に分解する問題と捉えることができる。これが，非負値行列因子分解と呼ばれる所以である。NMF のイメージを掴んでもらうために，（注 1）

スペクトログラムを行列と見なして NMF を適用した例を

図 1 に示す。

として係数が互いに独立になるように基底が求まる傾向にある。

4. NMF アルゴリズム
4. 1 正値行列因子分解と NMF
観測行列を二つの非負値行列の積で表そうという NMF の基
本概念自体は，Paatero らによって最初に提案されている [2]。

Paatero らは誤差行列 Y − HU の Frobenius ノルム (行列の要素の二乗和) で HU の Y からの乖離度 DEU (H, U ) を定義
し，さらに，H と U の各行列要素が非負であることを保証する目的で，各行列要素が 0 になろうとすると無限大のペナル
ティを課す対数障壁関数を定義し，これに適当な係数を掛けて乖離度 DEU (H, U ) に足したものを目的関数とした最適化アルゴリズムを提案している。Paatero らの方法は，最終的に得ら

3. NMF の基本性質

れる H と U が上記の障壁関数の性質より正値行列に限られる

NMF では通常，基底数 M は観測ベクトルの次元 K やデー

ことから，positive matrix factorization (PMF) と呼ばれる。

タの個数 N より小さく設定される。例えば M = K の場合，

一方，上述のような障壁関数を用いずとも行列要素の非負性をH = I (I は単位行列) であるような分解表現 Y = HU を得

保証しながら Y ≃ HU となる H と U を効率的に得る反復アることができるが，明らかにこの分解は意味をなさない。まルゴリズムが Lee らによって考案され [1]，これをきっかけにた，M = N の場合，U = I であるような分解表現 Y = HU

NMF が広く知れ渡るようになった。

を得ることができるが，この分解からも意味を見出せない。

4. 2 行列の乖離度

M < min(K, N ) のとき，NMF は観測行列 Y を低いランクの

NMF は，HU の Y からの乖離度を表す規準の定義に応じ行列で近似しようとしていることに相当し，その場合に求まるて異なる最適化問題に帰着する。Lee らは，前述の二乗誤差規準に加え，I ダイバージェンス規準の NMF アルゴリズムを導（注 1）：各時刻周辺における信号のスペクトルを時間順に並べて表示したものをスペクトログラムと呼ぶ。

出している [1]。当然ながら異なる乖離度規準を用いた場合では異なる解が最適解として得られる。このため，最適解が所望の—2—りである。
定義 1. θ = {θi }1<をパラメータとする目的関数 D(θ) に=i<=I対し，D(θ) = minθ̄ G(θ, θ̄) が成り立つとき，G(θ, θ̄) を D(θ)
の補助関数と定義する。
定理 1 (補助関数法). 補助関数 G(θ, θ̄) を，θ̄ に関して最小化するステップと，θ1 , · · · , θI に関して最小化するステップθ̄ ← argmin G(θ, θ̄)

(5)

θ̄

θi ← argmin G(θ, θ̄) (i = 1, . . . , I)

(6)

θi

を繰り返すと，目的関数 D(θ) の値は単調減少する。
図2

DEU/KL/IS (y|x) を x の関数と見たときのグラフ

小野氏による補助関数法に関する解説 [5] があるので，詳細解であるためには，背後にあるデータの生成プロセスに合った
適切な規準を設定することが重要である。例えば，パワースペ
クトルを対象とした NMF では板倉斎藤擬距離が乖離度規準として用いられることがある [3] が，これは 6. 2 で述べる信号波
形の生成プロセスに関する仮定に基づいている。

y, x ∈ R とすると，y の x からの二乗誤差，I ダイバージェ
ンス，板倉斎藤擬距離はそれぞれ

はそちらを参照して頂きたい。余談だが，著者は大学院生時代，不完全データに対する確率モデルのパラメータ推定手法である Expectation-Maximization アルゴリズムの拡張解釈を通し，Lee らとは独立に補助関数法を考案した。当時，研究室ス
タッフであった小野氏とその性質や補助関数の設計方法について多くの議論を交わし，ともにこれまで補助関数法の音響信号処理における様々な最適化問題への応用に取り組んできた。小野氏の解説記事にはその取り組みの一部が紹介されているので，DEU (y|x) = (y − x)y
DKL (y|x) = y log − y + xxyyDIS (y|x) = − log − 1
x
x
2

(2)

補助関数法の NMF 以外への適用例に興味のある読者は是非参(3)

照して頂きたい。

4. 4 DEU=KL (H; U ) を規準とした NMF(4)

上 述 の 補 助 関 数 法 の 原 理 に 基 づ き ，ま ず Lee ら に よ るで与えられる。いずれも，x = y のときにのみ 0 となり，x がDEU (H, U ) を規準とした NMF アルゴリズム [1] の導出のy から離れるほど増大する関数である。図 2 にそれぞれを x の概要を示す。DEU (H, U ) を展開し，H と U に依らない項を関数と見たときのグラフを示す。二乗誤差は y を中心に対称で省略すると，

あるのに対し，I ダイバージェンスと板倉斎藤擬距離は非対称

H,UDEU (H, U ) =であり，x が y を下回る場合に，より過大なペナルティを課す関数であることが分かる。また，板倉斎藤擬距離は y と x の比と書ける。ただし，

xk,n =

ある。これらを用いれば，HU の Y からの乖離度を，

X
k,n

˛X
“
”
˛
D· yk,n ˛hk,m um,nmで測ることができる (· は EU/KL/IS を表す)。

(7)

k,n

のみで表される関数であるため，y と x のスケールに非依存でD· (H, U ) =X(−2yk,n xk,n + x2k,n )X

hk,m um,n

(8)

m
z

で あ る 。ま た ，= は z に 関 係 す る 項 の み に 関 す る 等 号 を表 す。DEU (H, U ) に お い て x2k,n が ，H と U の 行 列 要 素

hk,1 , . . . , hk,M , u1,n , . . . , uM,n を含んだ非線形関数項である。

4. 3 補助関数法

ここで，この項に対し，行列要素 hk,m , um,n ごとの関数の和

これらを最小化する非負値の H と U を求めるのが目的であに分離した形をした上限関数を設けたい。上述のような上限関るが，いずれも非負制約つき非線形最適化問題であり，解析的数は，2 次関数が凸関数であることを利用し，以下の Jensen のに解くことはできない。もし，目的関数が行列要素 hk,m , um,n

不等式を用いることで設計可能である。Jensen の不等式とは，P任意の凸関数 g ，I 個の実数 x1 , . . . , xI ， i λi = 1 を満たす正ごとの関数の和に分離した形をしていれば，各行列要素ごとに非負値制約の中で解を探索することができて好都合であるが，残念ながらこれらの目的関数ではそうはなっていない。Lee ら
は，目的関数の上限となる補助関数を反復的に降下させること
で目的関数を間接的に降下させていく，補助関数法と呼ぶ方法
をベースに，補助関数として行列要素 hk,m , um,n ごとの関数の和に分離した形をとるものをうまく設計することで，当該制
約付き非線形最適化問題の解を見通し良く探索することができると考えた。補助関数の定義と補助関数法の原理は以下のとおの重み係数 λ1 , . . . , λ1 のもとで成り立つ“X” X
λi z i <λi g(zi )
g
=
i

(9)

i

(等号成立は z1 = · · · = zI ) のような不等式のことをいう。これを用いると，x2k,n に対し，

x2k,n <=

Xm

λk,m,n„hk,m um,nλk,m,n

«2

(10)

—3—

のような不等式を立てることができる。ただし，λk,m,n > 0,Pm λk,m,n = 1 であり，等号は

hk,m um,nλk,m,n =
xk,n

(11)

のとき成立する。以上より，DEU (H, U ) における x2k,n を式 (10)の右辺に置き換えることにより得られる関数 GEU (H, U , –)作れない。ここで，任意の微分可能な凹関数 g に対し，′
g(x) <= g(α) + (x − α)g (α)(17)(等号成立は x = α) が成り立つことを利用すると，
log xk,n <= log αk,n +1(xk,n − αk,n )
αk,n

(18)

は DEU (H, U ) の 補 助 関 数 の 要 件 を 満 た す。た だ し ，– =のような不等式が立てられ，行列要素の 1 次式で表された上{λk,m,n }K×M ×N とする。補助関数が得られれば，限関数を設計できる。以上の二つの不等式の等号成立条件は– ← argmin GEU (H, U , –)(12)

λk,m,n と αk,n がそれぞれ–H ← argmin GEU (H, U , –)(13)

λk,m,n =hk,m um,n, αk,n = xk,n
xk,n

(19)

H

U ← argmin GEU (H, U , –)

(14)

U

のときである。以上より，DIS (H, U ) における 1/xk,n の項を
式 (16) の右辺に，log xk,n の項を式 (18) の右辺に置き換えるを反復的に行えば DEU (H, U ) を単調に降下させていくこと

ことにより得られる関数 GIS (H, U , –, ¸) は DIS (H, U ) の補ができる。式 (12) については既に述べたとおり式 (11) で与え助関数の要件を満たす。あとは定理 1 に従って各更新式を求められる。式 (13) と式 (14) は hk,m >= 0, um,n >= 0 の制約のもとで解かれなければならない点に注意が必要だが，有り難いれば板倉斎藤擬距離規準の NMF アルゴリズムが導かれる。

ことに GEU (H, U , –) は行列要素ごとの関数の和に分かれた

バージェンスと呼ぶ規準で統一的に記述することができ，筆者

形をしているので，これらはいずれも行列要素ごとの一変数

らはこれを規準とした一般化 NMF アルゴリズムを [6] で導い関数の非負制約つき最小化問題に帰着し，∂GEU /∂hk,m = 0,ている。導出のポイントは前述の筆者らのアイディア [4] と同∂GEU /∂um,n = 0 を解くことで更新式を容易に得ることがで様であり，凸関数の項に対しては式 (9) を，凹関数の項に対しきる。

ては式 (17) の不等式を用いて上限関数を設計する，という方二乗誤差，I ダイバージェンス，板倉斎藤擬距離は，β ダイ

DKL (H, U ) を規準とした NMF アルゴリズムも以上と同様

針で閉形式の更新式からなる反復アルゴリズムを補助関数法ににして導くことができる。負の対数関数が凸関数であることに基づき導くことができる。以上のアイディアで導出される具体注意すると，− log xk,n の項に対し Jensen の不等式„«
Xhk,m um,n− log xk,n <
λk,m,n log=−
λk,m,nm

的なアルゴリズムは以下のとおりである。β = 0, 1, 2 のとき，

4. 4 と 4. 5 で示した板倉斎藤擬距離規準，I ダイバージェンス
規準，二乗誤差規準の NMF アルゴリズムと等価となる。

が立てられる。これにより行列要素 hk,m , um,n ごとの関数のβ ダイバージェンス規準の NMF アルゴリズム [6]

和に分離した形の補助関数を設けることができ，導出は以下同（ 1 ） H, U を非負値行列に初期設定する。

様である。

（ 2 ） 以下の更新を収束するまで繰り返す。

4. 5 DIS (H; U ) を規準とした NMF
次に，筆者らによる，板倉斎藤擬距離 DIS (H, U ) を規準とした NMF アルゴリズムの導出 [4] を以下に示す。H, U に依らない項を省略すると，DIS (H, U ) は«
X „ yk,n
H,UDIS (H, U ) =+ log xk,n
xk,n

0X

yk,n xβ−2k,n um,n

B nhk,m ← hk,m B
@ X

xβ−1
k,n um,n

1φ(β)
C
C
A

n

(15)

k,n

と書ける。先と同様の方針で補助関数法により当該非負制約つき最適化問題を解決するためには，DIS (H, U ) に対し，行列要
素 hk,m , um,n ごとの関数の和に分離した形をした上限関数を設計できるかどうかがやはり鍵となる。まず，逆数関数は正領域において凸関数であるため，1/xk,n の項に関しては，Jensen

0X

yk,n xβ−2k,n hk,m

B kum,n ← um,n B
@ X

xβ−1
k,n hk,m

1φ(β)
C
C
A

k

8
>
1/(2 − β)>><φ(β) = 1>>>:1/(β − 1)(β < 1)
(1 <
=β<
= 2)
(β > 2)

の不等式

1
xk,n

<
=

„ ﬃ
«
hk,m um,nλk,m,n 1
λk,m,nm

X(16)

ところで，文書データのクラスタリングを目的として開発された probabilistic Latent Semantic Analysis (pLSA) [7] は，が 成 り 立 つ 。た だ し ，λk,m,n は 先 と 同 様 ，λk,m,n > 0,Pm λk,m,n = 1 を満たす。次に，log xk,n の項に関しては，

NMF と類似した考え方をベースにしたデータ解析手法である。

正の対数関数は凹関数なので Jensen の不等式では上限関数がH と U を正規化するかどうかの違いを除けば pLSA は I ダイ

両者の関連については紙面の都合上割愛せざるを得なかったが，—4—バージェンス規準の NMF とアルゴリズム的に等価である。

形変換であるため，時間周波数分解により得られる複素スペク
トログラムも加法的である。しかし，複素スペクトログラムか5. 観測行列の生成プロセスら振幅 (またはパワー) スペクトログラムへの変換は非線形で以上のような規準で NMF を行うことが，観測データの背後あるため，(振幅およびパワー) スペクトログラムは実際には非にどのような生成プロセスを仮定した場合に正当化されるか，加法的である。端的に言えば，波形 A と波形 B の和のスペクということについて考えよう。二乗距離，I ダイバージェンス，

トルは波形 A のスペクトルと波形 B のスペクトルの和とは等板倉斎藤距離を規準とした NMF は，観測行列の要素 yk,n がしくなると限らない。従って，振幅 (またはパワー) スペクトロ

xk,n を平均とした正規分布，Poisson 分布，指数分布

グラムを加法的な成分に分解したところで，それぞれの成分がyk,n ∼ N (yk,n ; xk,n , σ 2 )(20)yk,n ∼ Poisson(yk,n ; xk,n )(21)

yk,n ∼ Exponential(yk,n ; xk,n )(22)

物理的に何に対応するかは定かではない。

NMF をスペクトログラムに適用するアプローチでは，周波に従って独立に生成されたと仮定した場合の H, U の最尤推定問題と各々等価である。ただし，

数成分比さえ同じであれば波形 (ゲインや位相スペクトル) が異なっていたとしても同一音と見なそう，という考え方がベースとなっているが，著者らはこの考え方を基に複素スペクトロ
グラムをモデル化し，観測複素スペクトログラムを加法的な成分に分解する「複素 NMF」と呼ぶ方法を提案している [9]。

N (z; µ, σ2−(z−µ)2 /2σ 2
1) = √2πσ2e(23)Poisson(z; µ) = µz e−µ /z! (z >
= 0)

(24)

Exponential(z; µ) = µ1 e−z/µ (z >
= 0)

(25)

6. 2 板倉斎藤擬距離規準の NMF [3]
先述のとおりパワースペクトル自体は加法的ではないが，各構成音の信号を互いに統計的に独立な確率変数と見なしたとき，それぞれのパワースペクトルの期待値は加法的となる。これは，である。このことは，以下により確かめられる。式 (20)〜(22) に

NMF において置かれるスペクトルの加法性の仮定が，期待値より立てられる xk,n の対数尤度 L(xk,n ) = log p(yk,n |xk,n ) はの意味では正当化される場合があるということを示唆している。

いずれも xk,n = yk,n のときに最大となる。すなわち，L(yk,n ) >=多重音信号中の各構成音の信号が短時間区間毎に平均が 0

L(xk,n ) である。よって，対数尤度差 L(yk,n ) − L(xk,n ) は，の (巡回) 定常 Gauss 過程に従って生成されたと仮定すると，

xk,n = yk,n のときにのみ 0 になる，yk,n と xk,n の近さを各区間における離散 Fourier 変換の各周波数成分は平均が 0

表す非負の尺度と見なせ，式 (20)〜(22) の場合の対数尤度差の複素正規分布に従う。すなわち，区間 n における m 番目

L(yk,n ) − L(xk,n ) は，それぞれ式 (2)〜(4) と等しい。

の構成音の周波数 k の成分 (構成音 m の複素スペクトログ(m)6. スペクトログラムモデルとしての NMF
Smaragdis らは，音楽音響信号の自動採譜を目的とし，音楽音響信号の (振幅またはパワー) スペクトログラムを観測行列と見なして NMF を適用し，楽音ごとのスペクトログラムに分

(m)

(m)

(m)

(m)

ラム) を sk,n とすると，sk,n は sk,n ∼ NC (sk,n ; 0, νk,n ) に2(m)
1 −|z−µ| /ν従う。ただし，NC (z; µ, ν) = πνeである。νk,n は構(m)
成音 m のパワースペクトログラムの期待値 E[|sk,n |2 ] を表す

パラメータである。ここで，観測信号の複素スペクトログラP (m)ム yk,n が yk,n = m sk,n のように構成音の複素スペクトロ(m′ )ルは加法的である)，という仮定と，各構成音は周波数成分比

グラムの和で与えられ，sk,n と sk,n (m =| m′ ) が互いに独P(m)立と仮定できるなら，yk,n ∼ NC (yk,n ; 0, m νk,n ) に従う。
P(m)
よって，xk,n = m νk,n と置くと，yk,n が与えられたもとでが時不変でゲインのみが時間変化する，という仮定の下で，観の xk,n の対数尤度は L(xk,n ) = − log πxk,n − |yk,n |2 /xk,n

測スペクトログラムから個々の構成音のスペクトル (周波数成と な る 。こ の 対 数 尤 度 は xk,n = |yk,n |2 の と き 最 大 と な分比) と各時刻におけるゲインを推定する問題は，NMF と同形の問題と見なせることが本アプローチのポイントである。例

る の で L(|yk,n |2 ) >= L(xk,n ) で あ る 。実 は ，対 数 尤 度 差
2
2
L(|yk,n | ) − L(xk,n ) >= 0 は，|yk,n | と xk,n の板倉齋藤擬

えば，音楽は平均律音階に該当する限られた種類の音高の音で距離 DIS (|yk,n |2 |xk,n ) と等しい。従って，各構成音のパワー構成されるため，各音高の音が時刻に依らずいつも同じ周波数

スペクトログラムの期待値を表す νk,n に関し，周波数成分比

成分比を有すると仮定できる場合，音楽音響信号のスペクトロが時不変であるような構造 νk,n = hk,m um,n を仮定すれば，グラムに NMF を適用すると各基底ベクトルが一つの音階の音

H = (hk,m )K×M と U = (um,n )M ×N の最尤推定問題は，観のスペクトルとなり，結合係数が各時刻におけるそれぞれの音

測パワースペクトログラム |yk,n |2 を要素にもつ行列 Y に対し階の音のゲインとなるように得られる。以上のアプローチの有板倉斎藤擬距離規準の NMF を行うことと等価となる [3]。

解する方法 [8] を提案している。瞬時瞬時の観測スペクトルが構成音のスペクトルの重みつき和によって表される (スペクト

(m)

(m)

(m)

効性は経験的に知られているが，上述の，核となる 2 つの仮定

6. 3 可変基底 NMF [10]〜[12]は実際には成り立たない。本章では，これらの仮定をより現実

NMF をスペクトログラムに適用するアプローチでは基本的に即したものに改変した改良モデルを紹介する。

に，分解したい各構成音の周波数成分比は時不変であるという

6. 1 複素 NMF [9]仮定が置かれるが，音楽音響信号を楽音ごとに分解する問題を

音の信号波形は加法的である。短時間 Fourier 変換やウェー扱う上では必ずしもこの仮定は成り立たない。歌声や音声はも

ブレット変換などの時間周波数分解は基本的には信号波形の線

ちろん，ピアノやヴァイオリンの各音符に対応する音響信号の

—5—

文

スペクトルは時々刻々と変化する。このように同一音源の音が持続していたとしてもスペクトルが時間変化する場合，通常[1]の NMF を適用しても，時刻ごとのスペクトルが別の音として[2]

ばらばらに分解されてしまうことになる。この問題意識のも[3]と，観測スペクトログラムを個々の持続音ごとに分解することを目的として，基底スペクトルが時間変化するよう改良された[4]

NMF の拡張モデルが提案されている [10]〜[12]。

[5]6. 4 その他の拡張モデル[6]

NMF の楽音分離や自動採譜への応用においては，音楽スペ
クトログラムを各楽音のノートごとのスペクトログラムに分解することが目的となるが，その成功の鍵は，楽音がもつ構造，
性質，傾向をいかに見つけ，NMF モデルにどう組み込めるか，

[7]
[8]
[9]

にある。このような動機からさまざまな拡張モデルが提案され[10]ている。例えば，各基底スペクトルのゲインが時間方向に連続的となるよう制約が組み込まれたもの [13]，各基底スペクトル[11]が調波構造をなすよう制約されたもの [14]，各基底スペクトルに自己回帰系によるソースフィルタモデルの拘束が組み込まれ[12]たもの [15]，基底スペクトルが音色特徴量空間においていくつ[13]かのクラスタを形成するよう基底のクラスタリングと NMF に[14]

よるスペクトログラム分解が同時に行われるもの [16]，などが検討されている。

[15]

NMF において基底数をいかに決定するかは重要課題の一つである。Cemgil らは，5. で述べたような NMF の生成モデル[16]としての解釈を通して，NMF の基底数決定問題を周辺尤度 (Hと U に関して周辺化したデータ行列の生成確率 p(Y )) に基づ[17]

くモデル選択の問題として定式化している [17]。モデル選択に[18]

より NMF の基底数を決定するためには，さまざまな基底数の[19]

もとで NMF を実行し，情報量規準や周辺尤度を算出して比較する手続きが必要となるが，ノンパラトリックベイズアプロー[20]

チによりパラメータ (基底行列と係数行列) とともにモデルの[21]

複雑度 (基底数) を一挙に推論しようという方法も近年提案されている [12], [18]。

[22]

7. 音響信号処理への応用[23]

NMF は，自動採譜の他，音声強調・分離 [19]，帯域拡張 [20]，[24]

音楽信号からのボーカル成分抽出 [21]・ドラム音抽出 [22]，音
声認識のための音素特徴量抽出 [23]，ホルマントトラッキン
グ [24]，エコーキャンセラ [25] などの音響信号処理問題に応用展開されており，筆者らも NMF のアルゴリズムをヒントにし

[25]
[26]

[27]

たブラインド残響除去法を提案している [26]。また，NMF は
モノラル信号の音源分離のための有効なアプローチとして認識[28]

されるようになって以来，そのモデル化の考え方は多チャネル[29]

信号の音源分離にも効果的であろうという期待から，NMF の多チャネル拡張に関する検討も積極的に進めている [27]〜[30]。

8. ま と め

[30]

[31]

本稿では，NMF の基本性質，アルゴリズムの導出方法，生成モデルとしての解釈，音響信号処理への応用とそのための拡

[32]
[33]

張モデルについて解説した。より理解を深めたい読者は他の著[34]

書や解説記事 (例えば [31]〜[34]) も是非とも参照して頂きたい。

献

D. D. Lee and H. S. Seung, “Algorithms for nonnegative matrix factorization,” in Adv. NIPS, pp. 556–562, 2000.
P. Paatero and U. Tapper, “Positive matrix factorization: A nonnegative factor model with optimal utilization of error estimates ofdata values,” Environmetrics, vol. 5, pp. 111–126, 1994.
C. Févotte, N. Bertin and J.-L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music
analysis,” Neural Computation, vol. 21, no. 3, pp. 793–830, 2009.
亀岡, 後藤, 嵯峨山, “スペクトル制御エンベロープによる混合音中の周期および非周期成分の選択的イコライザ,” 情処研報, 2006-MUS-66-13, pp. 77–84, 2006.
小野, “補助関数法による最適化アルゴリズムとその音響信号処理への応用,” 日本音響学会誌, Vol. 68, No. 11, pp. 566–571, 2012.
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,” in Proc.
MLSP, pp. 283–288, 2010.
T. Hofmann, “Probabilistic latent semantic analysis,” In Proc. UAI,
pp. 289–296, 1999.
P. Smaragdis and J. C. Brown, “Non-negative matrix factorization
for music transcription,” in Proc. WASPAA, pp. 177–180, 2003.
亀岡, 小野, 柏野, 嵯峨山, “複素 NMF: 新しいスパース信号分解表現と基底系学習
アルゴリズム,” 音講論 (秋), 2-8-13, pp. 657–660, 2008.
P. Smaragdis, “Non-negative matrix factor deconvolution; extractionof multiple sound sources from monophonic inputs,” in Proc. ICA
2004, pp. 494–499, 2004.
A. Ozerov, C. Févotte and M. Charbit, “Factorial scaled hidden
Markov model for polyphonic audio representation and source separation,” in Proc. WASPAA 2009, pp. 121–124, 2009.
中野, ルルー, 亀岡, 中村, 小野, 嵯峨山, “スペクトログラムのベイジアンノンパラメ
トリックモデリングに基づく音楽信号の解析,” 情処研報, 2011-MUS-91-6, 2011.
T. Virtanen, “Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria,”
IEEE Trans. ASLP, vol. 15, no. 3, pp. 1066–1074, 2007.
S. A. Raczynski, N. Ono and S. Sagayama, “Multipitch analisys withharmonic nonnegative matrix approximation,” in Proc. of ISMIR
2007, pp. 381–386, 2007.
H. Kameoka and K. Kashino, “Composite autoregressive system forsparse source-ﬁlter representation of speech,” in Proc. ISCAS2009,
pp. 2477–2480, 2009.
H. Kameoka, M. Nakano, K. Ochiai, Y. Imoto, K. Kashino and S.
Sagayama, “Constrained and regularized variants of non-negative
matrix factorization incorporating music-speciﬁc constraints,” inProc. ICASSP, pp. 5365–5368, 2012.
A. T. Cemgil, “Bayesian inference for nonnegative matrix factorization models,” Tech. Rep. CUED/F-INFENG/TR.609, University ofCambridge, 2008.
M. Hoﬀman, D. Blei and P. Cook, “Bayesian nonparametric matrix
factorization for recorded music,” in Proc. ICML, pp.439–446, 2010.
P. Smaragdis, B. Raj and M.V. Shashanka, “Supervised and semisupervised separation of sounds from single-channel mixtures,” inProc. ICA 2007, pp. 414–421, 2007.
P. Smaragdis and B. Raj, “Example-driven bandwidth expansion,”in Proc. WASPAA, pp. 135–138, 2007.
J.-L. Durrieu, G. Richard, B. David and C. Févotte, “Source/ﬁlter
model for unsupervised main melody extraction from polyphonic audio signals,” IEEE Trans. ASLP, vol. 18, no. 3, pp. 564–575, 2010.
M. Helén and T. Virtanen, “Separation of drums from polyphonicmusic using non-negative matrix factorization and support vector
machine,” in Proc. EUSIPCO, 2005.
A. Hurmalainen, J. Gemmeke and T. Virtanen, “Non-negative matrix
deconvolution in noise robust speech recognition,” in Proc. ICASSP,
pp. 4588–4591, 2011.
J.-L. Durrieu, J.-P. Thiran, “Sparse non-negative decompositionof speech power spectra for formant tracking,” in Proc. ICASSP,
pp. 5260–5263, 2011.
戸上, 川口, “準ブラインド非負行列分解を用いたマルチチャンネル非線形エコーキャ
ンセラ,” 音講論 (秋), 3-Q-9, pp. 757–758, 2009.
H. Kameoka, T. Nakatani and T. Yoshioka, “Robust speech dereverberation based on non-negativity and sparse nature of speech spectrograms,” in Proc. ICASSP, pp. 45–48, 2009.
A. Ozerov and C. Févotte, “Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation,” IEEE
Trans. ASLP, vol. 18, no. 3, pp. 550–563, 2010.
Y. Kitano, H. Kameoka, Y. Izumi, N. Ono and S. Sagayama, “Asparse component model of source sinals and its application to blindsource separation,” in Proc. ICASSP, pp. 4122–4125, 2010.
H. Sawada, H. Kameoka, S. Araki and N. Ueda, “New formulations
and eﬃcient algorithms for multichannel NMF,” in Proc. WASPAA,
pp. 153–156, 2011.
H. Sawada, H. Kameoka, S. Araki and N. Ueda, “Eﬃcient algorithms
for multichannel extensions of Itakura-Saito nonnegative matrix factorization,” in Proc. ICASSP, pp. 261–264, 2012.
A. Cichocki, R. Zdunek, A.H. Phan and S. Amari, “Nonnegative matrix and tensor factorizations: Applications to exploratory multi-way
data analysis and blind source separation,” Wiley, 2009.
亀岡, “非負値行列因子分解,” 計測と制御, vol. 51, no. 9, pp. 835–844, 2012.
澤田, “非負値行列因子分解 NMF の基礎とデータ／信号解析への応用,” 信学誌，vol. 95, no. 9，2012.
亀岡, “非負値行列因子分解の音響信号処理への応用,” 日本音響学会誌, Vol. 68,
No. 11, pp. 559–565, 2012.

—6—