解説：特集
計測・センシングのアルゴリズム

非負値行列因子分解
亀

岡 弘

＊＊
和＊，

＊東京大学 大学院情報理工学系研究科 東京都文京区本郷 7–3–1＊＊日本電信電話株式会社 NTT コミュニケーション科学基礎研究所
神奈川県厚木市森の里若宮 3–1＊Graduate School of Information Science and Technology, The Universityof Tokyo, 7–3–1 Hongo, Bunkyo-ku, Tokyo, Japan＊＊NTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation, 3–1 Morinosato Wakamiya, Atsugi, Kanagawa,
Japan

1.

＊E-mail: kameoka@hil.t.u-tokyo.ac.jp, kameoka.hirokazu@lab.ntt.co.jp
キーワード：多変量解析 (multivariate analysis)，次元圧縮 (dimensionality reduction)，信号分解 (signal decomposition)，補助関数 (auxiliary
function)．
C 2012 SICE
JL 0009/12/5109–0835 

はじめに

実世界には，パワースペクトル，画素値，頻度，個数な
ど，非負値で表わされるデータが多い．主成分分析や独立
成分分析などの多変量解析では，所与のデータを複数の加法的な成分に分解することを目的とするが，これと同様に
上述のような非負値のデータから構成成分を抽出することが役立つ場面が多い．たとえば，複数の音源の音響信号が混在する多重音のパワースペクトルから個々の音源のパワースペクトルをうまく取り出すことができれば，雑音除去や音源分離などに役立てることができるし，顔画像データを
目や鼻などの顔のパーツに該当する画像データにうまく分解することができれば，顔認証や顔画像合成などに役立てることができる．また，文書データ (文書に出現する各単語の個数のデータ) から「時事」や「スポーツ」や「経済」といった潜在的なトピックに該当するような単語ヒストグラ
ムの構成成分をうまく取り出すことができれば，各文書に対し自動インデキシングを行うことが可能となり，文書検索に大いに役立てることができる．以上のように，非負値のデータを加法的な構成成分に分解することを目的とした多変量解析手法を非負値行列因子分解 (Non-negative Matrix
Factorization; NMF)1) といい，さまざまな分野で近年注目を集めている．本稿では，NMF の定式化，基本的な性
質，アルゴリズムの導出方法，さまざまな視点からの解釈，
おもに音響信号処理の問題に焦点を当てた改良・拡張のア
イディアについて解説する．

2.

NMF とは以後，データをベクトルで表記することとする．たとえば，画像データであれば各ピクセルの画素値がベクトル要素となり，パワースペクトルであれば各周波数におけるパワー値がベクトル要素となる．今，N 個の非負値データベ
クトル y 1 , . . . , y N ∈ R≥0,K が与えられたとしよう．これ
らを観測ベクトルと呼ぶ．ただし，R≥0,K は K 次元の非負値ベクトル全体の集合を表わす．NMF では，これら観測ベクトルが，いずれも M 個の基底ベクトルの適当な重みつ
き和によって表わされたものと見なし，すべての観測ベク
トルを最も良く説明するような M 個の基底ベクトルおよ

計測と制御

第 51 巻 第 9 号

2012 年 9 月号

図 1 スペクトログラムを行列と見なして NMF を行った
場合の例．図のスペクトログラム Y を HU で表わそうとすることは，各時刻の観測スペクトルを，2 種類の基底
スペクトルの適当な重みつき和で表わそうとしていることに相当する．
び重み係数を推定することが目的である．よって，NMF では加法性が成り立つ量のみが対象となる．画素値やパワースペクトルはいずれも厳密には加法性が成り立つ量ではな
いが，NMF を適用する場面では近似的に加法性が成り立つと仮定されているということに注意が必要である．パワー
スペクトルに対して NMF を適用する上でのパワースペク
トルの非加法性についての議論は 6.1 節 で詳しく述べる．
以上の加法性に関する仮定の他に，基底ベクトルと重み係数がいずれも非負値であるという仮定が置かれる点も NMF
における重要なポイントである．すなわち，NMF は，観測ベクトル y n を基底ベクトル h1 , . . . , hM ∈ R≥0,K の非負結合 (結合係数 u1,n , . . . , uM,n が非負値の線形結合)

yn 

M


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) 式は835

できる．基底数 M は観測データが属する凸錐の次元に対応し，本来は未知であるため，観測データから適切な基底数を決定する問題は NMF において重要課題の一つである．
NMF における基底数の推定方法ついては 6.4 節で紹介する．

3.2 共起成分のグルーピング効果
NMF では，観測データの中で共起する成分をひとまと図 2 非負値基底の非負結合

Y  HU(2)と同じ意味である．このように NMF は，観測ベクトルを並べた行列を 2 つの非負値行列の積に分解する問題と捉えることができる．これが，非負値行列因子分解と呼ばれる所以である．NMF による分解表現のイメージを掴んでもらうために，（注1）スペクトログラムを行列と見なして NMFを適用した例を図 1 に示す．
基底ベクトルの非負性の仮定は，観測データを構成する成分も観測データと同じ属性の物理量であり，それゆえ同様に非負値であるべきという考えに基づいている．たとえば，パワースペクトルを構成する成分もパワースペクトルであるべきであり，負値をとるパワースペクトルは物理的に意味をなさない，という考え方である．
重み係数の非負性の仮定は，次章で詳しく述べるが，観測データ行列 Y から基底行列 H と係数行列 U を得る
際に，行列 U の要素がスパース (疎) になるよう誘導する効果をもたらす．係数行列がスパースになるということは，それぞれの観測ベクトルを少数の基底ベクトルだけで表わ
そうとしていることを意味し，どのような基底ベクトルが得られるかということに大きく関係する．

3.

NMF の基本性質

3.1 観測行列の低ランク近似
NMF では通常，基底数 M は観測ベクトルの次元 K や
データの個数 N より小さく設定される．たとえば M = Kの場合，H = I (ただし，I は単位行列) であるような分
解表現 Y = HU を得ることができるが，明らかにこの分解は意味をなさない．また，M = N の場合，U = I であるような分解表現 Y = HU を得ることができるが，この分解からも意味を見出せない．M < min(K, N ) のとき，
NMF は観測行列 Y を低いランクの行列で近似しようとしていることに相当し，その場合の基底行列と係数行列を求めることが重要である．幾何学的には，主成分分析 (特異値分解) では観測データが属する部分空間を見つけようとしているのに対し，NMF では観測データに良くフィットす
るような凸錐 (図 2 参照) を見つけようとしていると解釈（注1）

各時刻周辺における信号のスペクトルを時間順に並べて表示したものをスペクトログラムと呼ぶ．

836

めにしたものが基底ベクトルの推定結果になる傾向がある．
たとえば，図 3 のようにデータが 3 種類の成分 a，b，c から成っているとき，成分 a，b，c を基底と置けば当然デー
タを完全に表現することは可能である．ここでもし，成分a や成分 b が生起するときいつも成分 c も揃って生起する
ようなら，成分 a と成分 c，および，成分 b と成分 c をひとまとめにしたものを基底と置いても当該データを同様に表現可能である．この方が，少ない基底でデータを表現可
能な分，より「節約的」である．このように，共起する成分をひとまとめにして基底と置いた方が，基底を節約でき，節約した分を他の成分にフィットするのに充てられるようになるわけである．3.1 節で述べたように，NMF では少ない基底で観測データをできるだけ良く表わす必要があるため，上記のような基底が自然に得られる傾向にある．

3.3 係数行列をスパースにする効果
NMF では，係数行列の非負性の制約により，係数行列の要素がスパース (疎) になる傾向がある．これはベクトルの合成を考えるとイメージしやすい．観測ベクトルを近似するには，観測ベクトルと近い方向を向いた基底ベクトル以外の係数はできるだけ小さくした方が良い，ということである．スパース性は優ガウス性と近い概念であり，それゆえ統計的独立性とも概念的に関連が深い．3.2 節で述べた，共起する成分をひとまとめにしたものを基底にする働
きは，なるべく各基底の係数が互いにコヒーレントにならないように (相関をもたないように) 基底を決定しようとする働きと見なすことができる．独立成分分析では，係数の互いの独立性を規準として基底を推定するのが目的であるが，NMF では，非負制約による副次的な効果として係数が互いに独立になるように基底が求まる傾向にある．

3.4 分解の一意性について
一般に，X = HU となる非負値行列 H と U は一意に決まらない．HU を，正則行列 D を用いて H  = HD−1と U  = DU の積の形で表わすことができるが，H  = H,
U  = U ，かつ，H  と U  がいずれも非負値となるようなD が存在すれば，X = HU のような分解は一意に決まらない，ということになる．これを満たす D の中で自明なものは，置換行列，および，要素が正値の対角行列である．前者は，H と U に対し適当な列置換と行置換を行ったものが H  と U  である場合に相当し，後者は，H の各列に定数を乗じ，U の各行にその逆数を乗じたものが H と U  である場合に相当する．NMF において，後者のようなスケールの任意性を除く目的で，H の各列ベクトルのノルムが等しくなるような制約が置かれる場合が多い．

計測と制御

第 51 巻

第9号

2012 年 9 月号

図3

4.

NMF において共起成分がグルーピングされる効果

NMF アルゴリズム

4.1 正値行列因子分解と NMF
観測行列を 2 つの非負値行列の積で表わそうという NMFの基本概念自体は，Paatero らによって最初に提案されている2) ．Paatero らはこの中で，誤差行列 Y − HUの（注2）Frobenius ノルム (行列の要素の二乗和)
DEU (H, U ) = Y − HU 2F
2
 


=
hk,m um,n yk,n −

(3)

m

k,n

で HU の Y からの乖離度を定義し，さらに，H と Uの各行列要素が非負であることを保証する目的でB(H, U ) = −

log hk,m −k,m

log um,n(4)m,nなる対数障壁関数 (各行列要素が 0 になろうとすると無限大のペナルティを課す関数) を定義し，これに適当な係数を掛けて乖離度 DEU (H, U ) に足したものを目的関数とした最適化アルゴリズムを提案している．Paatero らの方法は，最終的に得られる H と U が上記の障壁関数の性質より正値行列に限られることから，正値行列因子分解 (Positive
Matrix Factorization; PMF) と呼ばれる．一方，上述のような障壁関数を用いずとも行列要素の非負性を保証しな
がら Y  HU となる H と U を効率的に得る反復アル
ゴリズムが Lee らによって考案され1), 3) ，これをきっかけに NMF が広く知れ渡るようになった．本章では，Lee らの方法がベースにする補助関数法と呼ぶ最適化手法の原理をまず紹介し，これに基づいて NMF アルゴリズムがどのようにして導かれるかを解説する．

4.2 行列の乖離度
NMF は，HU の Y からの乖離度を表わす規準の定義に応じて異なる最適化問題に帰着する．Lee らは，前述の二乗誤差規準に加え，（注3）I ダイバージェンス4) 規準の
NMF アルゴリズムを導出している3)．当然ながら異なる
乖離度規準を用いた場合では異なる解が最適解として得られる．このため，最適解が所望の解であるためには，背後（注2）

正確には重みつき二乗誤差規準が採用されている2) ．

（注3）一般化 Kullback-Leibler ダイバージェンスとも呼ばれる．

計測と制御

第 51 巻 第 9 号

2012 年 9 月号

図4

二乗誤差
DEU (y|x)，I ダ イ バ ー ジェン ス
DKL (y|x)，板倉斎藤擬距離 DIS (y|x) を x の関数と見たときのグラフにあるデータの生成プロセスに合った適切な規準を設定することが重要である．たとえば，パワースペクトルを対象とした NMF ではしばしば板倉斎藤擬距離5) が乖離度規準として用いられる6) が，これは信号波形の生成プロセスの仮定から導かれる．この理論的裏付けは 6.2 節 で詳しく述
べる．
y, x ∈ R とすると，y の x からの二乗誤差，I ダイバー
ジェンス，板倉斎藤擬距離はそれぞれDEU (y|x) = (y − x)2y
DKL (y|x) = y log − y + xxyyDIS (y|x) = − log − 1
x
x

(5)
(6)
(7)

で与えられる．いずれも，x = y のときにのみ 0 となり，xが y から離れるほど増大する関数である．図 4 にそれぞれを x の関数と見たときのグラフを示す．二乗誤差は y を中心に対称であるのに対し，I ダイバージェンスと板倉斎
藤擬距離は非対称であり，x が y を下回る場合に，より過大なペナルティを課す関数であることがわかる．また，板倉斎藤擬距離は y と x の比のみで表わされる関数である
ため，y と x のスケールに非依存である (両者を α 倍しても乖離度は等しいと解釈される)．

837

図 6 Jensen の不等式

θi ← argmin G(θ, θ̄) (i = 1, . . . , I)
θiを繰り返すと，目的関数 D(θ) の値は単調減少する．

図 5 補助関数法の原理
それぞれの規準を採用することがどのようなデータの生成プロセスを背後に仮定したことになっているかは 5 章 で詳しく述べることにして，ひとまずこれらを用いれば，HUの Y からの乖離度を，

D· (H, U ) =







D· yk,n 
hk,m um,n

k,n

m

(8)

で測ることができる．ただし，· は EU，KL，IS を表わす．

4.3 補助関数法
これらを最小化する非負値の H と U を求めるのが目的であるが，いずれも非負制約つき非線形最適化問題であり，解析的に解くことはできない．もし，目的関数が行列
要素 hk,m , um,n ごとの関数の和に分離した形をしていれば，各行列要素ごとに非負値制約の中で解を探索することができて好都合であるが，残念ながらこれらの目的関数ではそうはなっていない．Lee らは，目的関数の上限となる
補助関数を反復的に降下させることで目的関数を間接的に降下していく方法をベースにし，その補助関数として行列
要素 hk,m , um,n ごとの関数の和に分離した形をとるもの
をうまく設計することで，当該制約付き非線形最適化問題の解を見通し良く探索することができると考えた．補助関数の反復降下による目的関数の降下方法を補助関数法と呼び，音響信号処理分野で近年さまざまな最適化問題に適用
7)〜13)

されている．補助関数の定義と補助関数法の原理は以下のとおりである．
定義 1 θ = {θi }1≤i≤I を変量とする目的関数 D(θ) に対し，

D(θ) = min G(θ, θ̄)θ̄(9)が成り立つとき，G(θ, θ̄) を D(θ) の補助関数 (Auxiliary
function)，θ̄ を補助変数と定義する．
定理 1 (補助関数法) 補助関数 G(θ, θ̄) を，θ̄ に関して最小化するステップと，θ1 , . . . , θI に関して最小化するステップθ̄ ← argmin G(θ, θ̄)
θ̄

838

(11)

(10)

4.4 DEU/KL (H, U ) を規準とした NMF
上述の補助関数法の原理に基づき，まず Lee らによる
DEU (H, U ) を規準とした NMF アルゴリズム3) の導出を以下に示す．DEU (H, U ) を展開し，H と U に依らない
項を省略すると，H,U 

DEU (H, U ) =(−2yk,n xk,n + x2k,n )

(12)

k,n

と書ける．ただし，

xk,n =



hk,m um,n

(13)

m
z

である．また，= は z に関係する項のみに関する等号を表わす．(12) 式の中で，x2k,n が H と U の行列要素
hk,1 , . . . , hk,M , u1,n , . . . , uM,n を含んだ非線形関数項であることに気付く．この項に対し，行列要素 hk,m , um,n ごとの関数の和に分離した形をした上限関数を設けたい． 2次関数は凸関数であるため，以下の Jensen の不等式を用いることで，上述のような上限関数が設計可能である．
定理 2 (Jensen の不等式 (図 6)) 任意の凸関数 g ，I 個
の実数 x1 , . . . , xI ， i λi = 1 を満たす I 個の正値の重み係数 λ1 , . . . , λ1 のもとで，

g



 
λi zi ≤λi g(zi )

i

i

(14)

が成り立ち，z1 = · · · = zI のとき等号が成立する．
これを用いると，x2k,n に対し，

2
hk,m um,nλk,m,nm2
hk,m um,n≤λk,m,nλk,m,nm

x2k,n =

λk,m,n(15)のような不等式を立てることができる．ただし，λk,m,n > 0,m λk,m,n = 1 である．ここで重要なのは，右辺が行列要

計測と制御

第 51 巻

第9号

2012 年 9 月号

素 hk,m , um,n ごとの関数の和に分離した形をしているということである．また，もう一つ重要なポイントは，右辺をλk,m,n に関して最小化したとき等号が成立するということ
h u1,nhuM,n= · · · = k,Mであり，等号が成立するのは λk,1
λk,M,n ，k,1,nすなわち，λk,m,n がλk,m,n =hk,m um,n
xk,n

(16)

のときである．
以上より，DEU (H, U ) における x2k,n の項を (15) 式の右辺に置き換えることにより得られる関数

GEU (H, U , λ) =(17)
 h2k,m u2m,n 2− 2yk,nhk,m um,n +yk,nλk,m,nmmmax{ûm,n , 0} で与えられる．ところで，(21)，(22) 式は，H, U , λ の要素がすべて非負であればいずれも非負値となるため，H と U の初期値が非負であることが前提となっている場合は hk,m = ĥk,m , um,n = ûm,n である．
以上より，(18) 式を (19), (20) 式 に代入すると，以下のNMF アルゴリズムが導かれる．
二乗誤差規準の NMF アルゴリズム

1. H, U を非負値行列に初期設定する．
2. 以下の更新を収束するまで繰り返す．
yk,n um,nnhk,m ← hk,m 

k,n

は そ れ ぞ れ DEU (H, U ) の 補 助 関 数 の 要 件 を 満たす．
ただし，λ = {λk,m,n }K×M×N を表わすものとする．
GEU (H, U , λ) が得られれば，あとは定理 1 に従い，

λ ← argmin GEU (H, U , λ)(18)

H ← argmin GEU (H, U , λ)

(19)

λ

H

U ← argmin GEU (H, U , λ)U(20)を反復的に行えば DEU (H, U ) を単調に降下させていくことができる．(18) 式についてはすでに述べたとおり (16) 式で与えられる．(19) 式と (20) 式は hk,m ≥ 0, um,n ≥ 0の制約のもとで解かれなければならない点に注意が必要だ
が，ありがたいことに GEU (H, U , λ) は行列要素ごとの関数の和に分かれた形をしているので，これらはいずれも行
列要素ごとの一変数関数の非負制約つき最小化問題に帰着
し，容易に解くことができる．たとえば，GEU (H, U , λ)は hk,m ごとに見れば二次関数であり，これを最小化する
hk,m は

yk,n um,nĥk,m = nu2m,n /λk,m,n

(21)

n

のように解析的に求まる．この値がもし負であれば hk,m ≥0 の制約下で GEU (H, U , λ) が最小となるのは明らかに
hk,m = 0 のときである．また，um,n については
yk,n hk,mûm,n = k
h2k,m /λk,m,n(22)

xk,n um,nnyk,n hk,mkum,n ← um,n 

xk,n hk,mk

DKL (H, U ) を規準とした NMF アルゴリズムも以上と同様にして導かれる．DKL (H, U ) は，H,U 

DKL (H, U ) =(−yk,n log xk,n + xk,n )

(23)

k,n

と書けるが，この式の中で −yk,n log xk,n が hk,1 ,. . ., hk,M ,u1,n , . . . , uM,n を含んだ非線形関数項である．この項に関し，負の対数関数の凸性を利用して Jensen の不等式

− log xk,n ≤ −
λk,m,n log

mhk,m um,nλk,m,n(24)が立てられ，これにより行列要素 hk,m , um,n ごとの関数の和に分離した形の補助関数を設けることができる．これを用いて以下の NMF アルゴリズムを導くことができる．

I ダイバージェンス規準の NMF アルゴリズム
1. H, U を非負値行列に初期設定する．
2. 以下の更新を収束するまで繰り返す．
yk,n um,n xk,nhk,m ← hk,m n
um,n ← um,n

k



um,n

n

yk,n hk,m xk,n


hk,m

k

k

のとき GEU (H, U , λ) は最小となる．よって，(19) 式と(20) 式 は，それぞれ hk,m = max{ĥk,m , 0}, um,n =

計測と制御

第 51 巻 第 9 号

2012 年 9 月号

4.5 DIS (H, U ) を規準とした NMF
つぎに，著者らによる，板倉斎藤擬距離 DIS (H, U ) を839

規準とした NMF アルゴリズムの導出7) を以下に示す．H,
U に依らない項を省略すると，DIS (H, U ) は
H,U 



DIS (H, U ) =

k,n

yk,n
+ log xk,n
xk,n

1
xk,n



=1
≤



λk,m,nm

λk,m,n 1m
hk,m um,nλk,m,n
hk,m um,nλk,m,ng(x) ≤ g(α) + (x − α)g (α)(27)(等号成立は x = α) が成り立つことを利用すると，第二項の log xk,n に対し，
log xk,n ≤ log αk,n + (xk,n − αk,n )/αk,n

(29)
(30)

のとき等号成立するということである．
以上より，DIS (H, U ) における 1/xk,n の項を (26) 式の右辺に，log xk,n の項を (28) 式の右辺に置き換えることにより得られる関数

   yk,n λ2k,m,n

 hk,m um,nhk,m um,n
αk,n
m
m
k,n

− log yk,n + log αk,n − 2 (31)+は DIS (H, U ) の補助関数の要件を満たす．補助関数がひとたび得られれば，あとは定理 1 に従い，λ, α ← argmin GIS (H, U , λ, α)

840

⎛

yk,n um,n x2k,n(32)

um,n /xk,nn⎛

yk,n hk,m x2k,n⎜ kum,n ← um,n ⎜⎝ hk,m /xk,n

⎞1/2
⎟
⎟
⎠
⎞1/2
⎟
⎟
⎠

k

4.6 Dβ (H, U ) を規準とした NMF(5)〜(7) 式 に 3 種類の乖離度規準を挙げたが，これらの乖離度規準をパラメータ β で統一的に記述した β ダイ
バージェンスと呼ぶ規準14)

Dβ (y|x) = y(28)なる不等式が立てられる．ここでやはり重要になるのは，いずれの不等式も右辺が行列要素 hk,m , um,n ごとの関数の和に分離しているということ，また，λk,m,n と αk,n がλ,α1. H, U を非負値行列に初期設定する．
2. 以下の更新を収束するまで繰り返す．

(26)



GIS (H, U , λ, α) =板倉斎藤擬距離規準の NMF アルゴリズム⎜ nhk,m ← hk,m ⎜⎝ が成り立つ．ただし，先と同様，λk,m,n > 0, m λk,m,n =1 である．つぎに，第二項の log xk,n の上限関数を設計することが目標であるが，第一項と異なり，正の対数関数は凹関数なので，Jensen の不等式では上限関数が作れない．
しかしながら，任意の微分可能な凹関数 g に対し，αk,n = xk,n(34)

を反復的に行えば DIS (H, U ) を単調に減少させることができる．(32) 式は (29), (30) 式のとおりであり，あとは先と同様に H と U の更新式を求めれば，以下の板倉斎藤
擬距離規準の NMF アルゴリズムが導かれる．

hk,m um,nλk,m,n =xk,nU ← argmin GIS (H, U , λ, α)U(25)数の和に分離した形をした上限関数を設計できるかどうかがやはり鍵となる．まず，逆数関数が正領域において凸関数であることに注意すると，第一項の 1/xk,n に関しては，
Jensen の不等式

(33)

H



と書ける．先と同様の方針で補助関数法により DIS (H, U )を目的関数とした非負制約つき最適化問題を解決するためには，DIS (H, U ) に対し，行列要素 hk,m , um,n ごとの関H ← argmin GIS (H, U , λ, α)y β − xβy β−1 − xβ−1
−
β−1
β

(35)

がある．ただし，β は β = 0, β = 1 の実数である．
lim (xβ − y β )/β = log(x/y) であることを利用すれば，
β→0(35) 式は β → 0 のとき板倉斎藤擬距離，β → 1 のとき I ダイバージェンス，β = 2 のとき二乗誤差となることが確認できる．これを (8) 式のように HU の Y からの乖
離度の規準 Dβ (H, U ) とした NMF アルゴリズムが Kompass らや中野らによって提案されている15), 16) ．Kompass
らは 1 ≤ β ≤ 2 においてのみ，中野らはすべての β において収束性が保証される NMF アルゴリズムを導出しており，以下に中野らによる導出16) の概略を示す．
(35) 式の第一項の (y β−1 − xβ−1 )/(β − 1) は，x の関数と見れば β ≤ 2 のとき凸関数，β ≥ 2 のとき凹関数となる．一方，第二項の −(y β − xβ )/β は β ≤ 1 のとき凹関数，β ≥ 1 のとき凸関数となる．前述の著者らのアイディ
ア7) と同様，β の条件に応じて凸関数の項に対しては (14)式を，凹関数の項に対しては (27) 式の不等式を適用することで Dβ (H, U ) の補助関数を設計することができる．詳細は省略するが，あとは前述の板倉斎藤擬距離規準の NMF
アルゴリズムの導出に従えば，以下の β ダイバージェンス
規準の NMF アルゴリズムが導かれる．

計測と制御

第 51 巻

第9号

2012 年 9 月号

log p(yk,n |xk,n ) はいずれも xk,n = yk,n のときに最大となる．すなわち，L(yk,n ) ≥ L(xk,n ) である．よって，対数尤度差 L(yk,n ) − L(xk,n ) は，xk,n = yk,n のときにのみ 0 になる，yk,n と xk,n の近さを表わす非負の尺度と見なすことができる．実は，(36)〜(39) 式の場合の対数尤度
差 L(yk,n ) − L(xk,n ) は，それぞれ (5)〜(7) 式, (35) 式とβ ダイバージェンス規準の NMF アルゴリズム
1. H, U を非負値行列に初期設定する．
2. 以下の更新を収束するまで繰り返す．
⎛yk,n xβ−2k,n um,n⎜ nhk,m ← hk,m ⎜
⎝ 
n

⎛

yk,n xβ−2k,n hk,m⎜ kum,n ← um,n ⎜⎝ kϕ(β) =xβ−1
k,n um,nxβ−1k,n hk,m

⎞ϕ(β)
⎟
⎟
⎠

等しい．

6.

⎞ϕ(β)

（注4）

⎟
⎟
⎠

⎧
⎪
⎪1/(2 − β) (β < 1)⎨1
(1 ≤ β ≤ 2)⎪⎪⎩
1/(β − 1) (β > 2)

以上のアルゴリズムは，β = 0, 1, 2 のとき，前述の板倉
斎藤擬距離規準，I ダイバージェンス規準，二乗誤差規準の NMF アルゴリズムと等価である．

5.

観測行列の生成プロセス

以上のような規準で NMF を行うことは，観測データの背後にどのような生成プロセスを仮定したことになっているだろうか．実は，二乗誤差，I ダイバージェンス，板倉斎
藤擬距離，β ダイバージェンスを規準とした NMF は，観測行列 Y の各行列要素 yk,n が，xk,n を平均とした正規分布，Poisson 分布，指数分布，Tweedie 分布

yk,n ∼ Normal(yk,n ; xk,n , σ 2 )(36)yk,n ∼ Poisson(yk,n ; xk,n )(37)

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

yk,n ∼ Tweedie(yk,n ; xk,n , φ)(39)

に従って独立に生成されたと仮定した場合の H, U の最尤推定問題とおのおの等価である．ただし，
221e−(z−μ) /2σ
Normal(z; μ, σ 2 ) = √2πσμz e−μPoisson(z; μ) =(z ≥ 0)z!
1Exponential(z; μ) = e−z/μ (z ≥ 0)
μ

(40)
(41)
(42)

1

(43)
Tweedie(z; μ, φ) = a(z, φ)e φ (zρ(μ)−κ(μ))
⎧
⎧
β−1
β
⎨ μ −1 (β = 0)⎨ μ −1 (β = 1)β−1β, κ(μ) =ρ(μ) =⎩log μ (β = 0)⎩log μ(β = 1)
である．このことは，以下により確かめられる．(36)〜(39) 式により立てられる xk,n の対数尤度 L(xk,n ) =

計測と制御

第 51 巻 第 9 号

2012 年 9 月号

音響信号処理への応用

モノラル音源分離のアプローチとして多重音信号のスペクトログラム (振幅スペクトログラムまたはパワース
ペクトログラム) を観測行列と見なして NMF を適用する
方法17), 18) が Smaragdis らと Virtanen によって提案され，それ以降 NMF は，自動採譜17) ，音声強調・分離19) ，帯域拡張20) などさまざまな問題への有効なアプローチとして近年応用展開されている．
瞬時瞬時の観測スペクトルが構成音のスペクトルの重みつき和によって表される (スペクトルは加法的である)，という仮定と，各構成音は周波数成分比が時不変でゲインのみが時間変化する，という仮定の下で，観測スペクトログ
ラムから個々の構成音のスペクトル (周波数成分比) と各時
刻におけるゲインを推定する問題は，NMF と同形の問題と見なすことができる．たとえば，音楽は平均律音階に該当する限られた種類の音高の音で構成されるため，各音高の音が時刻に依らずいつも同じ周波数成分比を有すると仮定できる場合，音楽音響信号のスペクトログラムに NMFを適用すると各基底ベクトルが一つの音階の音のスペクトルとなり，結合係数が各時刻におけるそれぞれの音階の音のゲインとなるように得られる．以上のアプローチの有効性は経験的に広く知られているが，上述の，核となる 2 つの仮定は実際には成り立たない．本章では，これらの仮定
をより現実に即したものに改変した改良モデルを紹介する．

6.1 複素 NMF21)
音の信号波形は加法的である．短時間 Fourier 変換やウェーブレット変換などの時間周波数分解は基本的には信号波形の線形変換であるため，時間周波数分解により得られる複素スペクトログラムも加法的である．しかし，複素ス
ペクトログラムから振幅 (またはパワー) スペクトログラムへの変換は非線形であるため，振幅スペクトログラムもパ
ワースペクトログラムも実際には非加法的である．端的に言えば，波形 A と波形 B の和のスペクトルは波形 A のス
ペクトルと波形 B のスペクトルの和とは等しくなると限らない．したがって，振幅 (ないしパワー) スペクトログラム
を加法的な成分に分解したところで，それぞれの成分が物（注4）異なる位置に配置された複数のマイクロホンで取得した観測信号から個々の音源信号を分離することを目的とした多チャンネル音
源分離では，音源の空間的な位置の違いを手がかりにすることができるが，モノラル信号からの音源分離は，これを手がかりにすることができないため，より難しい不良設定問題となる．

841

理的に何に対応するかは必ずしも定かではない．
NMF をスペクトログラムに適用するアプローチではス
ペクトログラムが行列の積でモデル化されるが，このモデルでは，各構成音の周波数成分比は時不変という仮定が置かれていることになる．つまり，周波数成分比さえ同じであれば波形 (ゲインや位相スペクトル) が異なっていたとしても同一音と見なそう，という考え方がベースとなっている．著者らはこの考え方を基に複素スペクトログラムをモ
デル化し，観測複素スペクトログラムを加法的な成分に分解する「複素 NMF」と呼ぶ方法を考案した21)．以下，その原理について概説する．
複素スペクトログラムの加法性より，M 個の音から成る
多重音の複素スペクトログラムは

xk,n =

M


(m)

m=1

sk,n

(44)
(m)

のように各構成音の複素スペクトログラム sk,n の和で表わすことができる．各構成音の周波数成分比が時不変と仮(m)定することは，構成音の振幅スペクトログラム |sk,n | を(m)|sk,n | = hk,m um,n(45)と表わすことと等しく，これを (44) 式に代入すると，多重音の複素スペクトログラムのモデル

xk,n =

(m)

hk,m um,n ejφk,n

(46)

m

という目的関数を最小化する問題を考え，これを局所的に最小化する反復アルゴリズムを補助関数法に基づき導出した．ただし，0 < γ < 2， k hk,m = 1 とする．導出の詳
細や実験結果についてはこちらの文献21) を参照されたい．

6.2 板倉斎藤擬距離規準の NMF6)
先述のとおりパワースペクトル自体は加法的ではないが，各構成音の信号が互いに統計的に独立のときそれぞれのパ
ワースペクトルの期待値は加法的となる．これは，NMF に
おいて置かれるスペクトルの加法性の仮定が，期待値の意
味では正当化される場合があるということを示唆している．
平均が 0 の定常 Gauss 過程に従う信号は，Fourier 変換領域において各周波数成分が独立に 0 平均の複素正規分布に従う．したがって，多重音信号中の各構成音の信号が短時
間区間ごとに平均が 0 の (巡回) 定常 Gauss 過程に従って生成されたと仮定すると，各区間における離散 Fourier 変換の各周波数成分は平均が 0 の複素正規分布に従う．すなわち，
区間 n における m 番目の構成音の周波数 k の成分 (構成(m)(m)
音 m の複素スペクトログラム) を sk,n とすると，sk,n は
(m)

k,n

k,t

(47)

842

(48)
2

(m)

値 E[|sk,n |2 ] を表わすパラメータである．ここで，観測信号の複素スペクトログラム yk,n が(m)(m)
が得られる．ただし，j は虚数単位であり，φk,n = arg(sk,n )yk,n − xk,n 2 + 2ηD(H, U , φ) =|um,n |γ(m)に従う．ただし，NormalC (z; μ, ν) = (1/πν)e−|z−μ| /ν で(m)ある．νk,n は構成音 m のパワースペクトログラムの期待

yk,n =である．以上のモデルを観測複素スペクトログラム yk,n (本節における yk,n は，5 章までの yk,n と異なり複素数であることに注意されたい．) にフィットするように未知パラ(m)
メータ hk,m , um,n , φk,n を推定することが目的であるが，NMF と同様に各基底スペクトル h1,m , . . . , hK,m が独立に生起する個々の音のスペクトルに対応するようになるためには係数 um,n がスパースになるかどうかが鍵となる．し
たがって，以上のモデルのもとで，できるだけ係数 um,nをスパースに，かつ，xk,n が yk,n にフィットするような(m)hk,m , um,n , φk,n を推定することができれば，NMF と同様の基底スペクトル獲得機能を実現できる可能性がある．
ただし，3.3 節で述べたように NMF では非負制約の副次的効果により係数 um,n が自然にスパース化される傾向があったのに対し，上記のモデルでは NMF とは違い実部と虚部における構成成分同士の減算を許すため，そのままでは係数 um,n がスパース化される効果は望めないことに注意が必要である．著者らは，yk,n と xk,n の二乗誤差に係
数 um,n に関するスパース正則化項を付加した(m)sk,n ∼ NormalC (sk,n ; 0, νk,n )

 (m)
sk,n

(49)

m

のように構成音の複素スペクトログラムの和で与えられ，
(m)
(m )
sk,n と sk,n (m = m ) が互いに独立と仮定できるなら，
yk,n は (48), (49) 式より

 (m) yk,n ∼ NormalC yk,n ; 0,
νk,n

(50)

m



(m)

に従う．よって，xk,n = m νk,n と置くと，yk,n が与えられたもとでの xk,n の対数尤度はL(xk,n ) = − log πxk,n −|yk,n |2
xk,n

(51)

と書かれる．この対数尤度は xk,n = |yk,n |2 のとき最大となるので L(|yk,n |2 ) ≥ L(xk,n ) である．実は，対数尤度差
L(|yk,n |2 ) − L(xk,n ) ≥ 0 は，|yk,n |2 と xk,n の板倉齋藤
擬距離 DIS (|yk,n |2 |xk,n ) と等しい．従って，各構成音のパ(m)
ワースペクトログラムの期待値を表わす νk,n に関し，周波数成分比が時不変であるような構造(m)

νk,n = hk,m um,n(52)を仮定すれば，H = (hk,m )K×M と U = (um,n )M×N の最尤推定問題は，観測パワースペクトログラム |yk,n |2 を計測と制御

第 51 巻

第9号

2012 年 9 月号

要素にもつ行列 Y に対し板倉斎藤擬距離規準の NMF を行うことと等価となる6) (4 章における yk,n は本節における |yk,n |2 に対応していることに注意されたい)．以上のよ(m)うな生成モデルの見地に立てば，sk,n を潜在変数と見なして Expectation-Maximization (EM) アルゴリズムの考え方に従って H と U の最尤解を局所探索するアルゴリズムを導くことができる6) が，4 章で示したアルゴリズムの方が収束が早いことが実験的に確認されている16)．

く説明する，各 HMM の状態ごとの出力スペクトル，状態
遷移系列，時変ゲインを推定する問題として定式化される．
なお，このモデルは，生成モデルの観点で見れば Factorial
HMM28) の一種と見なすことができる．

6.3 NMFD22) ・可変基底 NMF24)〜27)
NMF をスペクトログラムに適用するアプローチでは基6.4 その他の拡張モデル・他手法との関連
NMF において基底数をいかに決定するかは応用と基礎の両面で重要課題の一つである．Schmidt らや Cemgil らは，5 章で述べたような NMF の生成モデルとしての解釈
を通して，NMF の基底数決定問題を周辺尤度 (H と U に関して周辺化したデータ行列の生成確率 p(Y )) に基づくモ

本的に，分解したい各構成音の周波数成分比は時不変であるという仮定が置かれるが，音楽音響信号を楽音ごとに分解する問題を扱う上では必ずしもこの仮定は成り立たない．
歌声や音声はもちろん，ピアノやヴァイオリンの各音符に対応する音響信号のスペクトルは時々刻々と変化する．このように同一音源の音が持続していたとしてもスペクトルが時間変化する場合，通常の NMF を適用しても，時刻ごとのスペクトルが別の音としてばらばらに分解されてしま

デル選択の問題として定式化している29), 30) ．モデル選択に
より NMF の基底数を決定するためには，さまざまな基底数のもとで NMF を実行し，情報量規準や周辺尤度を算出
して比較する手続きが必要となる．これに対し，パラメー
タの推論と同時にモデルの複雑度の推論を行える枠組であるノンパラトリックベイズアプローチ31) により，NMF の基底数を推論する方法も提案されている27), 32)〜34) ．
NMF をスペクトログラムに適用するアプローチは，モう．そこで，観測スペクトログラムを各楽音の持続音ごとに分解することを目的として，楽音スペクトルの時変構造を取り入れた NMF の拡張モデルが提案されている．
(m)
各楽音の振幅 (またはパワー) スペクトログラム sk,n を時不変なスペクトルに時変のゲインを乗じたものとしてモデル化する代わりに，ある時間幅をもったスペクトログラ
ムに時変なゲインを畳み込んだもの

ノラル信号の音源分離のための有効なアプローチであるため，そのモデル化の考え方は多チャンネル信号の音源分離にも効果的であろうという期待から，NMF の多チャンネ
ル拡張に関する検討が近年進められている35)〜39)．
ところで，文書データのクラスタリングを目的として開発された probabilistic Latent Semantic Analysis (pLSA)40)
は，NMF と類似したアイディアに基づくデータ解析手法である．両者の関連については紙面の都合上割愛せざるをえなかったが，実は pLSA は I ダイバージェンス規準の NMFと H と U を正規化するかどうかの違いを除けばアルゴリ
ズム的に等価である．

(m)sk,n =

N
−1

hk,m,n−n um,n

(53)

n =0

としてモデル化し，観測振幅 (またはパワー) スペクトロ
 (m)グラム yk,n が与えられたもとで yk,n 
m sk,n となるように H = (hk,m,n )K×M×N  と U = (um,n )M×Nを推定する，Non-negative Matrix Factor Deconvolution(NMFD) と呼ぶ方法が提案されている22) ．また，複素版も提案されており，音声と背景音楽を分離する技術に応用されている23) ．
NMFD では，各楽音スペクトルの時間変化の仕方はどの発音時においても同じである，という仮定が暗に置かれるが，この仮定も必ずしも実際には成り立ない．たとえばビブラートのレート (揺れの周期) や深さは発音のたびに変わるかもしれない．そこで，各楽音スペクトルの変化のタ
イミングはいつも同じとは限らないが，変化の順序だけは保たれている，という仮定のもと，基底スペクトルが内部状態の遷移に応じて時間変化するよう改良された NMF の拡張モデルが近年提案されている24)〜27)．このモデルでは，
M 個の楽音スペクトログラムがそれぞれ隠れマルコフモ
デル (Hidden Markov Model; HMM) からの出力系列と見なされ，観測スペクトログラムはそれらの重ね合わせと解釈される．このモデルの下で観測スペクトログラムを楽音ごとに分離する問題は，観測スペクトログラムを最も良

計測と制御

第 51 巻 第 9 号

2012 年 9 月号

7.

まとめ

NMF は，音響信号処理やデータマイニングや画像処理
や機械学習などの分野で市民権を得始めているデータ解析
ツールであり，本稿では，NMF の基本性質，アルゴリズムの導出方法，生成モデルとしての解釈，音響信号処理への応用とそのための拡張モデルについて解説した．NMF に関する解説記事は澤田氏のもの41) もあるのでより理解を深めたい読者はぜひとも参照していただきたい．
（2012 年 7 月 27 日受付）
参

考 文

献

1） D.D. Lee and H.S. Seung: Learning the parts of objects with
nonnegative matrix factorization, Nature, 401, 788/791(1999)
2） P. Paatero and U. Tapper: Positive matrix factorization: A
non-negative factor model with optimal utilization of error
estimates of data values, Environmetrics, 5, 111/126 (1994)
3） D.D. Lee and H.S. Seung: Algorithms for nonnegative matrix factorization, in Adv. NIPS, 556/562 (2000)
4） I. Csiszár: I-divergence geometry of probability distributions and minimization problems, The Annals of Probabil-843

ity, 3–1, 146/158 (1975)
5） 板倉：統計的手法による音声分析合成系に関する研究，博士論
文，名古屋大学大学院工学研究科 (1972)
6） C. Févotte, N. Bertin and J.-L. Durrieu: Nonnegative matrix factorization with the Itakura-Saito divergence. Withapplication to music analysis, Neural Computation, 21–3,
793/830 (2009)7） 亀岡，後藤，嵯峨山：スペクトル制御エンベロープによる混合
音中の周期および非周期成分の選択的イコライザ，情報処理学
会研究報告，2006-MUS-66-13, 77/84 (2006)
8） H. Kameoka: Statistical approach to multipitch analysis,
Ph.D. Thesis, The University of Tokyo, Mar. (2007)9） 亀岡，鎌本，原田，守谷：予測誤差の Golomb-Rice 符号量を最小化する線形予測分析，信学論，J91-A–11, 1017/1025 (2008)
10） H. Kameoka, T. Nakatani and T. Yoshioka: Robust speech
dereverberation based on non-negativity and sparse natureof speech spectrograms, in Proc. ICASSP, 45/48 (2009)11） N. Ono, H. Kohno, N. Ito and S. Sagayama: Blind alignment of asynchronously recorded signals for distributed microphone array, in Proc. WASPAA, 161/164 (2009)12） N. Ono and S. Miyabe: Auxiliary-function-based independent component analysis for super-Gaussian sources, in
Proc. LVA/ICA, 165/172 (2010)13） N. Yasuraoka, H. Kameoka, T. Yoshioka and H. G. Okuno:
I-divergence-based dereverberation method with auxiliary
function approach, in Proc. ICASSP, 369/372 (2011)14） S. Eguchi and Y. Kano: Robustifying maximum likelihood
estimation, Technical report, Institute of Statistical Mathematics, Research Memo. 802 (2001)
15） R. Kompass: A generalized divergence measure for nonnegative matrix factorization, Neural Computation, 19–3,
780/791 (2007)16） M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono
and S. Sagayama: Convergence-guaranteed multiplicativealgorithms for non-negative matrix factorization with betadivergence, in Proc. MLSP, 283/288 (2010)17） P. Smaragdis and J. C. Brown: Non-negative matrix factorization for music transcription, in Proc. WASPAA, 177/180(2003)
18） T. Virtanen: Sound source separation using sparse codingwith temporal continuity objective, in Proc. ICMC, 231/234(2003)
19） P. Smaragdis, B. Raj and M.V. Shashanka: Supervised andsemi-supervised separation of sounds from single-channel
mixtures, in Proc. ICA 2007, 414/421 (2007)20） P. Smaragdis and B. Raj: Example-driven bandwidth expansion, in Proc. WASPAA, 135/138 (2007)21） 亀岡，小野，柏野，嵯峨山：複素 NMF: 新しいスパース信号分解
表現と基底系学習アルゴリズム，音講論 (秋)，2-8-13, 657/660(2008)
22） P. Smaragdis: Non-negative matrix factor deconvolution;
extraction of multiple sound sources from monophonic inputs, in Proc. ICA 2004, 494/499 (2004)
23） 北野，亀岡，柏野，小野，嵯峨山：複素 NMFD による音声抽出マスクの設計と背景音楽抑圧への応用，音講論 (春)，3-9-6,
719/720 (2009)24） A. Ozerov, C. Févotte and M. Charbit: Factorial scaled hidden Markov model for polyphonic audio representation andsource separation, in Proc. WASPAA 2009, 121/124 (2009)25） 中野，北野，ルルー，亀岡，小野，嵯峨山：可変基底 NMF に基づく音楽音響信号の解析，情処研報，2010-MUS-84-10 (2010)

84426） G.J. Mysore, P. Smaragdis and B. Raj: Non-negative hidden Markov modeling of audio with application to source
separation, in Proc. LVA/ICA 2010, 140/148 (2010)27） 中野，ルルー，亀岡，中村，小野，嵯峨山：スペクトログラムのベイジアンノンパラメトリックモデリングに基づく音楽信号の解析，情処研報，2011-MUS-91-6 (2011)28） Z. Ghahramani and M.I. Jordan: Factorial hidden Markov
models, Machine Learning, 29, 245/273 (1997)29） M.N. Schmidt, O. Winther and L.K. Hansen: Bayesian nonnegative matrix factorization, in Proc. ICA, 540/547 (2009)30） A.T. Cemgil: Bayesian inference for nonnegative matrix factorization models, Technical Report CUED/FINFENG/TR.609, University of Cambridge (2008)31） P. Orbanz and Y. W. Teh: Bayesian nonparametric models,in Encyclopedia of Machine Learning. New York: Springer(2010)
32） M.N. Schmidt and M. Mørup: Inﬁnite non-negative matrix
factorization, in Proc. EUSIPCO (2010)33） M. Hoﬀman, D. Blei and P. Cook: Bayesian nonparametric matrix factorization for recorded music, in Proc. ICML,
439/446 (2010)34） K. Yoshii and M. Goto: A nonparametric Bayesian multipitch analyzer based on inﬁnite latent harmonic allocation,
IEEE Trans. ASLP, 20–3, 717/730 (2012)35） A. Ozerov and C. Févotte: Multichannel nonnegative matrix factorization in convolutive mixtures for audio source
separation, IEEE Trans. ASLP, 18–3, 550/563 (2010)36） Y. Kitano, H. Kameoka, Y. Izumi, N. Ono and S. Sagayama:
A sparse component model of source sinals and its application to blind source separation, in Proc. ICASSP, 4122/4125(2010)
37） H. Sawada, H. Kameoka, S. Araki and N. Ueda: New formulations and eﬃcient algorithms for multichannel NMF,
in Proc. WASPAA, 153/156 (2011)38） 武田，亀岡，澤田，荒木，宮部，山田，牧野：混合 DOA モデルに基づく多チャンネル複素 NMF による劣決定ブラインド音
源分離，音講論 (春)，2-1-9, 747/750 (2012)39） H. Sawada, H. Kameoka, S. Araki and N. Ueda: Eﬃcient algorithms for multichannel extensions of Itakura-Saito nonnegative matrix factorization, in Proc. ICASSP, 261/264(2012)
40） T. Hofmann:
Probabilistic latent semantic analysis,
In Proc. Uncertainty in Artiﬁcial Intelligence (UAI’99),
289/296 (1999)41） 澤田：非負値行列因子分解 NMF の基礎とデータ／信号解析への応用，電子情報通信学会誌，95–9 (2012)［著 者 紹 介］
かめ

亀

おか

岡

ひろ

弘

かず

和君

1978 年 6 月 25 日生．2002 年東京大学工学部計数工学科卒業．07
年東京大学大学院情報理工学系研究科博士課程修了．情報理工学博士．
同年日本電信電話株式会社入社．NTT コミュニケーション科学基礎研究所配属．11 年東京大学大学院情報理工学系研究科客員准教授．音声・音楽を対象とした音響信号処理・機械学習の研究に従事．日本音響学会，情報処理学会，電子情報通信学会，IEEE の会員．

計測と制御

第 51 巻

第9号

2012 年 9 月号