2-8-13複素 NMF： 新しいスパース信号分解表現と基底系学習アルゴリズム∗○亀岡弘和 † , 小野順貴 ‡ , 柏野邦夫 † , 嵯峨山茂樹 ‡( NTT コミュニケーション科学基礎研究所, ‡ 東京大学大学院情報理工学系研究科)†1序論

基底のスパースな線形結合 (結合係数のほとんどがほぼ 0 値) ですべての観測データを良く近似できるようなコンパクトな基底系 (このような基底系を「スパースな分解表現」と呼ぶ。) を統計的，情報論的な規準に基づいて構成するアプローチを総称してスパー
ス信号分解と呼び，スパースコーディング (SC) [1] や非負値行列分解 (NMF) [2] などがその代表例である。
スパース信号分解では，主成分分析や独立成分分析などの多変量解析手法と同様，個々の基底が観測データを構成する基本要素成分となるように学習することが目的となるが，スパース信号分解の特質は，結合係数のスパース性の制約により，観測データの中に混在する特徴的なパターンが個々の基底となって表出される点にある。この性質より，SC も NMF も，モノラル音響信号中に混在する構成音を事前知識なしに分離抽出できる可能性をもった方法論として近年期待され [4, 5, 6]，雑音抑圧などの音声強調 [7] や狭帯域信号からの音声復元 [8] といった通信システムへの用途にも応用され始めている。
NMF が音響特徴抽出法として有効な理由は，パワーないし振幅スペクトルを観測データとして基底分解を行っているためであると考えられる [4]。すな
わち，NMF では，似た形の振幅スペクトルはすべて同一基底で表現可能となる点がポイントであり，このために多様で複雑な現象である実世界音響信号をコンパクトな基底系で表現することができるのである。
振幅スペクトルは，信号がどのような周波数成分で構成されているかを表した情報であるため，同一な振幅スペクトルをもつ音は同一属性と捉える NMF の考え方は，基底の冗長性の観点で多分に合理的である。
また，基底系の学習の際に設ける基底および結合係数に対する非負制約が，副次的に基底系の解を自然にスパースな分解表現へと誘導する効果がある点や，乗法更新アルゴリズム [2] と呼ぶ非常に効率的な反復学習アルゴリズムが存在する点も特徴的である。しかし，NMF は混合音中の構成音の特徴を大まかに捉えるという点においては有効とはいえ，振幅スペクトルは非加法的であるため，振幅スペクトル基底の線形結合による観測データのモデル化をベースにするNMF では混在音の振幅スペクトルをどの程度正確に表出できるかは必ずしも明らかではない。
一方で，時間領域信号基底の線形結合により観測信号をモデル化する SC に基づくアプローチ [5] は，加法性が成り立つ量を基底と考える点で NMF より精∗緻なモデリングに基づいている反面，基底系の冗長性に難点がある。すなわち，波形が異なっても振幅スペクトルが同一ならば同じ属性の音と見なせる NMFとは異なり，波形の微妙な差異に応じて別々の基底が必要になる。シフト不変 SC と呼ぶアプローチ [9] は，基底系に拘束を与えて，時間シフトの関係にある波形群を同一属性の音と見なせるように工夫したものであるが，ある波形に対して時間シフトした波形の集合は，その波形と振幅スペクトルが同形の全波形の中のごく一部に過ぎないため，この分解表現もまた NMF のそれに比べて冗長性は依然大きい。
以上の NMF と SC に対する問題意識を動機として，本稿では，複素スペクトル (加法性の成り立つ) 領域でのモデル化に準拠しながら NMF のように観測信号の中に混在する個々の特徴的な振幅スペクトルパターンを表出できるような性質をもった，新しいタイプのスパース信号分解法「複素 NMF」を提案する。

2

複素 NMF の原理

複素 NMF は，スパース信号分解としての効果を本質的に保持しながら，ベースとなるモデリングが精緻である点で SC の利点を，個々の特徴的な振幅スペクトルパターンを表出できる点と効率的な反復学習アルゴリズムが導ける点で NMF の利点を併せ持つ。以下，音響信号の分解表現をまず定式化する。

2.1

スパース信号分解表現

任意の音響信号の時間周波数成分を Fx,t とする。
ただし，x は周波数ないし対数周波数に，t は時刻に対応するインデックスである。ここで，Fx,t を，Fx,t =K∑
k=1

ak,x,t =∑|ak,x,t |ejφk,x,t(1)kのように要素 ak,x,t の和に分解し，さらに |ak,x,t | を|ak,x,t | = Hk,x Uk,t , Hk,x ≥ 0, Uk,t ≥ 0(2)のように Hk,x と Uk,t の積に分解する。分解のスケー
ルの任意性を除くため，ここで∑Hk,x = 1 (k = 1, · · · , K)
(3)xを仮定しておく。Hk,x は，時刻に依らずグローバルに決定される振幅スペクトルパターンであり，一方，k 番目の振幅スペクトルのアクティビティ(係数)Uk,tComplex NMF: A new sparse signal representation and basis set optimization. by KAMEOKA Hirokazu(NTT Corporation), ONO Nobutaka (University of Tokyo), KASHINO Kunio (NTT Corporation) andSAGAYAMA Shigeki (University of Tokyo)

日本音響学会講演論文集- 657 -

2008年9月

と位相スペクトル φk,x,t はそれぞれ時刻ごとに自由度をもつ。すなわち，実世界の音響信号は，高々K 種類の振幅スペクトルの信号だけで構成されており，それぞれの係数と位相スペクトルが自由に時変して実現したものと見なしているわけである。このモデルのもとで，SC の考え方と同様に，できるだけ係数 Uk,tをスパースに，かつ，観測信号とのモデル化誤差を小さくするように振幅スペクトル Hk,x と係数 Uk,t を学習することで，NMF のように似た振幅スペクトル構造が一つの基底となって表出されるような性質をもった複素スペクトル領域のスパース信号分解が可能になると考えられる。以上より，∑Fx,t =Hk,x Uk,t ejφk,x,t(4)kのような形のスパース信号分解表現を提案する。
ところで，NMF では，時刻に依らずグローバルに決定される振幅スペクトルパターン Hk,x が K種類あり，それぞれの係数 Uk,t が非負制約のもとで自由に時変していると見なしてスペクトログラムをモデル化していることに相当する。すなわち，∑Fx,t = k Hk,x Uk,t なので，F ≡ (Fx,t )X×T , HT ≡(Hk,x )K×X , U ≡ (Uk,t )K×T とすると，F = HU のように行列の形で書けるわけである。一方，位相スペクトルが φk,x,t = ϕk,x のような時不変の拘束をもつ
とき，Vk,x ≡ Uk,x ejϕk,x と置けて，V ≡ (Vk,x )K×Xとすると，F = VU と書けるので，SC のモデルと本質的に一致する。これに対し，式 (4) のモデルは，3 階のテンソル ejφk,x,t の存在により NMF や SC と同じような行列表現ができないことから，従来のスパース信号分解モデルとは性質が根本的に異なる。

2.2基底系学習アルゴリズム

今，観測音響信号の時間周波数成分 Yx,t が与えられたとき，Hk,x , Uk,t , φk,x,t を最適推定する方法について述べる。Y ≡ (Yx,t )X×T , F ≡ (Fx,t )X×T , H ≡(Hk,x )K×X , U ≡ (Uk,t )K×T , φ ≡ (φk,x,t )K×X×T とする。さて，ここで，観測とモデルとの間にYx,t = Fx,t + ²x,t(5)が成り立ち，上記のモデル化誤差 ² を，分散が σ 2 の複素ガウス分布に従う白色雑音と仮定すると，H, U, φの尤度 P (Y |H, U, φ, σ 2 ) はP (Y |H, U, φ, σ 2 ) =∏ 1 − |Yx,t −Fx,t |2σ2e2πσx,t(6)とすると，事後確率 P (H, U, φ, σ 2 |Y ) は，P (H, U, φ, σ 2 |Y ) ∝ P (Y |H, U, φ, σ 2 )P (U )と書ける。ただし，p, b は一般化正規分布の形状を規定する定数であり，0 < p < 2 のとき P (U ) は優ガウス的となり，スパースネスを測るための適切な尺度になる。式 (9) より，事後確率 P (H, U, φ, σ 2 |Y ) をH, U, φ に関して最大化する問題は，σ 2 , b をあらかじめ定めておく定数とすると，

f (H, U, φ) ≡とし，P (H), P (φ), P (σ 2 ) はいずれも一様分布，スパース性を表す U の事前確率を一般化正規分布

P (U ) =日本音響学会講演論文集∏|U|p1− k,tbp)
e2Γ 1 + p1 bk,t((8)(10)k,tを H, U, φ に関して最小化することと等しいので，以上の最適化問題は以下のようにまとめられる。

minimizesubject tof (H, U, φ)∑Hk,x = 1 (k = 1, · · · , K)(11)x

この最適化問題を解くアルゴリズムを補助関数法により導く。その準備として，まず，補助関数法の原理を示す。目的関数 G(z) に対し，

G(z) = min G+ (z, z̄)

(12)

z̄

が成り立つとき，G+ (z, z̄) を G(z) の補助関数，z̄ を補助変数と定義すると，次の定理が成り立つ．
定理 1 (補助関数法). 補助関数 G+ (z, z̄) を，z̄ に関して最小化するステップと，z に関して最小化するス
テップを繰り返すと，目的関数値は単調収束する。
証明： z̄t+1 = argminz̄ G+ (zt , z̄) とし，zt+1 =argminz G+ (z, z̄t+1 ) とする．ただし，t は反復計
算 の ス テップ 数 を 表 す．z = zt , z̄ = z̄t から z = zt+1 , z̄ = z̄t+1 に更新されたときに，G(z) が増加しないことを示す。明らかに G(zt ) =G+ (zt , z̄t+1 ) であり，zt+1 = argminz G+ (z, z̄t+1 ) より，G+ (zt , z̄t+1 ) ≥ G+ (zt+1 , z̄t+1 ) である．更に補助関数の定義より，G+ (zt+1 , z̄t+1 ) ≥ G(zt+1 ) なのだか
ら，結局，G(zt ) ≥ G(zt+1 ) である．
以上の補助関数法を適用するため，式 (10) で与えられる目的関数 f (H, U, φ) の補助関数を次の 2 つの不等式 (証明略) を用いて導く。
¯¯¯¯2 ∑ ¯Ȳk,x,t − Hk,x Uk,t ejφk,x,t ¯2¯Yx,t − Fx,t ¯ ≤(13)βk,x,t|Uk,t |p ≤(7)∑¯∑¯
¯Yx,t −Fx,t ¯2 +2λ|Uk,t |px,tで与えられる。また，H, U , φ, σ 2 は独立

P (H, U, φ, σ 2 ) = P (H)P (U )P (φ)P (σ 2 )(9)p|Ūk,t |2k
p−2Uk,t 2 + |Ūk,t |p −p|Ūk,t |p2(14)∑ただし，式 (13) は k Ȳk,x,t = Yx,t の条件のもとで，式 (14) は 0 < p < 2 の条件のもとで成り立つ。また，∑βk,x,t は 0 < βk,x,t < 1,k βk,x,t = 1 を満たす任意の定数である。Ȳ ≡ (Ȳk,x,t )K×X×T , Ū ≡ (Ūk,t )K×T- 658 -2008年9月

とし，式 (13), (14) より，f (H, U, φ, Ȳ , Ū ) ≡++λ∑(¯¯∑ ¯Ȳk,x,t − Hk,x Uk,t ejφk,x,t ¯2βk,x,tk,x,tp|Ūk,t |p−2 Uk,t 2 + 2|Ūk,t |p − p|Ūk,t |p)k,tとすれば f (H, U, φ) ≤ f + (H, U, φ, Ȳ , Ū ) が成り立ち，()Ȳk,x,t = Hk,x Uk,t ejφk,x,t + βk,x,t Yx,t − Fx,t (15)Ūk,t = Uk,t(16)のとき f (H, U, φ) = f + (H, U, φ, Ȳ , Ū ) となるため，
f + (H, U, φ, Ȳ , Ū ) は補助関数の定義を満たす。
以上の結果を用いて，H, U, φ の更新式が求められる。Hk,x に関しては，本来，式 (11) の条件で f + を最小化する必要があるが，ここでは簡単のため，まず無制約で f + を最小化する Hk,x を解析的に求めたのちに式 (11) を満たすように規格化する。
∂f +まず， ∂H= 0 を解くと，k,x∑ Uk,tHk,x =[]∗Re Ȳk,x,tejφk,x,tβk,x,t∑ Uk,t 2

ttのときに f + (H, U, φ, Ȳ , Ū ) は φk,x,t に関して最小になる。従って，ejφk,x,t の更新式は[][]Re Ȳk,x,t
Im Ȳk,x,tȲk,x,tjφk,x,t¯¯¯¯¯ (20)e= ¯+j ¯= ¯¯Ȳk,x,t ¯Ȳk,x,t ¯Ȳk,x,t ¯で与えられる。
この更新式より Hk,x と Uk,t の更新値の正値性が容易に示される。式 (17), (18) に式 (20) を代入すると，¯¯∑ Uk,t ¯Ȳk,x,t ¯βk,x,tHk,x = t(21)∑ Uk,t 2βk,x,tt¯¯∑ Hk,x ¯Ȳk,x,t ¯βk,x,tx
Uk,t =(22)∑ Hk,x 2¯¯p−2+ λp¯Ūk,t ¯
βk,x,txとなり，Hk,x , Uk,t の初期値が正であれば Hk,x , Uk,tは常に正値に更新されることが分かる。
以上より、アルゴリズムは次のようにまとめられる。

1. 初期設定: H, U, φ を初期設定
2. Ȳ の更新: 式 (15) により Ȳ を更新
0
3. φ の更新: 式 (20) により ejφk,x,t を更新
4. H の更新: 式 (17) により H を更新後，式 (11)を満たすように規格化
5. U の更新: 式 (18) により U を更新
6. 2. に戻る

(17)

βk,x,tを得る。ここで，∂ 2 f + /∂Hk,x ∂Hk0 ,x0 ≥ 0 なので式(17) は f + を最小にする更新である。
∂f +続いて，Uk,t についても， ∂U= 0 を解くと，k,t∑ Hk,xUk,t =[]∗Re Ȳk,x,tejφk,x,tβk,x,t∑ Hk,x 2+ λp|Ūk,t |p−2βk,x,txx

2.3(18)を得る。ここで，∂ 2 f + /∂Uk,t ∂Uk0 ,t0 ≥ 0 なので式(18) は f + を最小にする更新である。
次に，φk,x,t の更新式を導く。まず，f + のうち φk,x,tに依らない項をまとめて c と置くと，f + は+f (H, U, φ, Ȳ , Ū ) =∑c−2|Ak,x,t | cos(φk,x,t − Ck,x,t )k,x,tNMF と等価になるアルゴリズムの条件以上で導出した複素 NMF の基底系学習アルゴリズ
ムは，特定の条件下で Lee らが導出した NMF の乗法更新アルゴリズム [2] と等価となる。まず，位相ス
ペクトル φ の更新を，式 (20) により反復的に行う代わりに，ここでは，φ を初期設定時にejφk,x,t =(23)と置き，以後その値のまま固定して H, U のみを反復的に更新する場合について考察する。これにより，
Yx,t ∑
Fx,t =Hk,x Uk,t(24)|Yx,t |kと書ける。ただし，Ȳk,x,t Hk,x Uk,tβk,x,t]][[Re Ȳk,x,tIm Ȳ¯ , sin Ck,x,t = ¯ k,x,t¯cos Ck,x,t = ¯¯¯Ȳk,x,t ¯Ȳk,x,t ¯

Ak,x,t =である。cos(φk,x,t − Ck,x,t ) = cos φk,x,t cos Ck,x,t +sin φk,x,t sin Ck,x,t = 1 のとき，すなわち，

cos φk,x,t = cos Ck,x,t , sin φk,x,t = sin Ck,x,t (19)となるから，式 (15) にこれを代入すると Ȳ は∑()][Hk,x Uk,tn Hn,x Un,t+ βk,x,t 1 −Ȳk,x,t = Yx,t|Yx,t ||Yx,t |のような更新がなされることになる。これと式 (23)を式 (17) に代入すると，右辺は[)](∑ Uk,t∑Hn,x Un,tHk,x Uk,t + βk,x,t |Yx,t | −βk,x,ttn∑ Uk,t 2t

日本音響学会講演論文集

Yx,t|Yx,t |- 659 -βk,x,t

2008年9月

Fig. 1

目的関数値の収束の様子

Frequency bin #50040030020010050100150200150200

Frame #

Frequency bin #500400300200

Fig. 310050100

図 2 における H, U の学習結果

Frame #Fig. 2

観測信号 Yx,t (上段) と反復計算 30 回後のモデル
Fx,t (下段) のスペクトログラム∑となる。ところで βk,x,t は，βk,x,t > 0, k βk,x,t = 1を満たす範囲で任意に定めて良い定数なので，H やU とともに値を更新したとしてもアルゴリズムの収束性には何ら影響を与えない。そこで，βk,x,t をHk,x Uk,tβk,x,t = ∑Hn,x Un,t(25)nのように更新することにすると，Hk,x の更新則は，∑Uk,t |Yx,t |Hk,x ← Hk,x ∑ttUk,t∑(26)
Hn,x Un,tnと書け，Lee らが導出した Frobenius ノルム規準のもとでの乗法更新式 [2] と一致する。U に関しても上記と同様のことが言える。以上の考察より，複素 NMFアルゴリズムは，(1) NMF と同等に高効率で安定的な動作をすること，(2) βk,x,t を固定するのではなく H, U , φ とともに式(25) のとおりに更新するのが効果的であること，(3) (2) を行う場合，φ を式 (23) に固定すれば NMFと等価なアルゴリズムになるので，反復計算の初期段階で H, U を NMF により推定しておけば局所解問題が極めて効果的に回避できること，

レーム周期 32ms，Hanning 窓) により計算し，提案法の条件については，K = 10，p = 1.2，λ =∑21−p/2
× 10−5 ，反復計算回数 30 とし
x,t |Yx,t | /Kた。H, U の初期値はいずれも乱数生成により設定し，ejφ の初期値は Yx,t /|Yx,t | とした。2.3 節での議論のとおり，βk,x,t は式 (25) により更新した。
図 1 は，以上の条件のもとで目的関数 f (H, U, φ)が反復計算ごとにどのように推移していくかを示したものである。2.2 節で導出したアルゴリズムは厳密に収束性が保証されているわけではないが，実際にはほぼ収束していく様子が確認できる。
図 2 は，Yx,t と 30 回の反復計算によって推定された Fx,t をそれぞれスペクトログラム表示したものである。このときの H, U の推定結果を併せて図 3 に示す。H のそれぞれの推定結果から，調波構造らしいスペクトルの特徴的構造が事前知識なしに表出されていることが分かり，提案法が NMF に似た性質をもつことが確認できた。

4

本稿では，複素スペクトル (加法性の成り立つ) 領域でのモデル化に準拠しながら NMF のように観測信号の中に混在する個々の特徴的な振幅スペクトルパターンを表出できるような性質をもった，新しいタイプのスパース信号分解法「複素 NMF」を提案し，音声データに対してその基本動作を確認した。

参考文献などが予想される。

3動作実験[1] Olshausen & Field,Nature, 381, 607–609, 1996.
[2] Lee & Seung, Nature, 401, 788–791, 1999.
[3] Plumbley, IEEE Trans. NN, 14(3), 534–543, 2003.
[4] Smaragdis et al., WASPAA2003, 177–180, 2003.
[5] Abdallah et al., IEEE Trans. NN, 17(1), 179–196,2006.

ATR 音声データベースの音声データ (女性話者)を用いて提案法の基本動作の確認を行った。Yx,t はSTFT(標本化周波数 16kHz， フレーム長 64ms，フ日本音響学会講演論文集

まとめ- 660 -[6] Schmidt, Interspeech2006, 2614–2617, 2006.
[7] Wilson et al., ICASSP2008, 4029–4032, 2008.
[8] Bansal et al., Interspeech2005, 1505–1508, 2005.
[9] Blumensath et al., ICASSP2004, 5, 497–500, 2004.

2008年9月