２２．NMFで音声を分解してみよう - keep learning blog（キープラーニングブログ）

――難問は、それを解くために適切かつ必要な単位まで分割せよ。――

NMFで音声を分解してみよう
- 非負値行列因子分解（NMF）

NMFで音声を分解してみよう

2月になりました。もう3月は目の前、春はすぐそこ・・・と思っていたところ、全くそんな気配はなく普通に1月より寒くなりました。しかし、今年の北欧は暖冬のようです。年によってはマイナス15℃で海が凍ることもあるらしいです。

研究の方は一つ論文を書き終わったので、別のトピックを検討中です。次はディープラーニングじゃない機械学習をやりたいと思っています。今回は、その検討中にいくつか勉強した機械学習アルゴリズムを記事にまとめたいと思います。

f:id:yuki0718:20200216011207j:plain

非負値行列因子分解（NMF）

これまでディープニューラルネットワーク（DNN）を中心に研究を進めてきましたが、ご存知のとおりDNNは大量の教師データと高いコンピュータ処理能力が必要なので、ノートPCのCPUで気軽に実行したりリアルタイムで実行したりする用途には不向きです。

また、TensorflowやPytorchのようなライブラリさえあれば、ネットワークを構築して終わってしまうので、数学的に解釈性の高いアルゴリズムが好きな人はどこか物足りないかもしれません。

そこで今回は、音声処理やコンピュータビジョンの世界で発展を遂げ、最も美しく定式化されたアルゴリズムの一つである、非負値行列因子分解（NMF: Non negative matrix factorization）をご紹介します。

NMFの有用性と定式化

現代社会では、画像の画素値、文書中の単語の出現頻度、（音波の）振幅を二乗したパワースペクトルetc...、のように、必ず正の値しか取らないデータを取り扱う場面が多々あります。実はこの非負値という性質をうまく利用すると、一つのデータ群を異なる性質を持った複数の群に分解できる可能性があります。

顔の画像を目や耳、鼻などのパーツに分解できれば顔認証やモンタージュに役立ちますし、文書中の単語群から「経済」や「政治」といったトピックに分解できればニュースの分類にも役立ちます。音声の世界であればボーカル音声とBGMの分離や通話中の雑音除去などに応用できるかもしれません。

それでは、肝心の分解はどうやって行うのでしょうか。それこそが今回のメインテーマである、データの非負値性を利用した分解手法、非負値行列因子分解（Non-negative matrix factorization）です。

NMFのアイディア自体は古くから知られていましたが、それが現在のように広く用いられるようになったきっかけは、2000年にLeeらが発表した論文で提案された更新アルゴリズムでした。ここではそれを導くまでの流れに沿ってご紹介します。

まず、データの例として音声のスペクトログラムを考えます。スペクトログラムは、STFT（短時間フーリエ変換）の振幅絶対値を並べた画像のようなもので、行が周波数（行数K）、列が時間フレーム（列数N）を表す（K×N）の大きさを持った行列Sでした。

f:id:yuki0718:20190728035947p:plain — ギター演奏のスペクトログラム

NMFの定式化は、この行列Sを、行数K×列数Mの行列Hと行数M×列数Nの行列Uとの行列積で「近似」できると仮定するところから始めます。すなわち、 $S \approx H \times U$ です。

ただし、任意の基底数Mは、Sの行列サイズよりも小さい（M＜min(K, N)）ものとします。でないと単位行列の自明解に収束してしまうからです。この分解された行列Hを基底行列（Basements matrix）と呼び、行列Uを活性化行列（Activation matrix）と呼びます。

NMFの数学的解釈

音声のスペクトログラムの場合、基底行列Hには「音声に含まれる様々な音色の寄せ集め」という物理的意味があります。すなわち、行数K×列数Mの行列Hは、各周波数に分布を持ったM本の基底ベクトルを横に並べたもの、H=(音色1, 音色2, …, 音色M)だと解釈できます。他方、活性化行列Uは、各時刻における音量（基底ベクトルの重み係数）だと解釈できます。

数学的には、近似値H×Uの計算は正の方向に伸びる基底ベクトルの正の線形結合なので、H×Uの各列ベクトルはK次元空間の原点からHを構成する各基底ベクトルの方向に伸びた三角錐内の点に対応します。

データが三角錐の外にある場合はそのデータから最も近い三角錐内の点（往々にして三角錐の表面）が最もよい近似値になります。

三角錐の表面は、K次元空間中でいずれかの基底ベクトルの重みが0の領域なので、一般的に活性化行列Uはほとんどの成分が0の疎行列（Sparce matrix）になりがちです。反転して考えればHも疎行列になりやすい性質があります。

このNMFの数学的性質に起因する疎性（Sparce性）については、NTT研究所の亀岡先生が日本統計学会誌に寄稿された特集記事の説明が素晴らしく分かりやすいです。

NMFの更新アルゴリズム

話を戻して、とにかくスペクトログラムSのような二次元の行列データをHとUの行列積で近似する、というのがNMFの出発点でした。近似というからには、元のスぺクトログラムSと近似値H×Uとの間に誤差（Error）が生じます。

NMFの場合はこの誤差のことを乖離度D（Deviation）と呼びます（D=0なら近似値が正解と完全一致）。乖離度Dの候補としてすぐに思いつくのは、以下のような各行列要素の差を二乗して和を取った量、いわゆる二乗和誤差（Squared Euclidean distance）でしょう。

$\displaystyle D(S \| H \times U)=\sum_{k=1}^{K} \sum_{n=1}^{N} \{ S_{k,n}-(H \times U)_{k,n} \} ^{2}$

この乖離度を非負値の制約条件下で最小化するような解を求めれば、基底行列Hと活性化行列Uが得られ、音声データSを（近似的に）音色ごとの基底ベクトルに分解することができます。しかし、残念ながらこの最小化問題には解析解が存在しません。

この問題をコンピュータで数値的に解くため、巧妙かつ強力な数学的テクニックとして、補助関数法（Auxiliary function method）と呼ばれる手法が広く知られています。

補助関数法では、最小化したい乖離度 $D(\theta)$ （ $\theta$ はHとUのような求めたい変数）に対して、あえて新たな補助変数 $\lambda$ を導入し、 $D(\theta) \leqq D_{+}(\theta, \lambda)$ を常に満たすような補助関数 $D_{+}(\theta, \lambda)$ を設計します。補助関数が常にDの上限を押さえるので、補助関数を最小化すれば頭を押さえつけられて自ずとDも小さくなります。

さて、二乗和誤差の乖離度Dの場合、Jensenの不等式（1906年にデンマーク人Johan Jensenが提唱）という有名な数学公式を利用して、以下のような補助関数を設計できます。

$\displaystyle D(H, U) \leqq D_{+}(H, U, \lambda)=\sum_{k=1}^{K} \sum_{n=1}^{N} \left( S_{k,n}^{2}-2S_{k,n}\sum_{m=1}^{M} H_{k,m}U_{m,n}+\sum_{m=1}^{M}\frac{H_{k,m}^{2}U_{m,n}^{2}}{\lambda_{k,m,n}} \right)$

ただし等号が成立するのは、 $\displaystyle \lambda_{m,k,n}=\frac{H_{k,m}U_{m,n}}{\sum_{m'=1}^{M}H_{k,m'}U_{m',n}}$ のとき

詳細な計算は上述した特集記事にまとまっているのでここでは割愛しますが、これをH,Uのそれぞれについて偏微分して0と置くことで、以下のようなHとUの逐次更新式を得ることができます。

$\displaystyle H_{k,m}^{new}=H_{k,m}\frac{\sum_{n=1}^{N} S_{k,n}U_{m,n}^{T}}{\sum_{n=1}^{N} X_{k,n}U_{m,n}^{T}}$

$\displaystyle U_{m,n}^{new}=U_{m,n}\frac{\sum_{k=1}^{K} H_{k,m}^{T}S_{k,n}}{\sum_{k=1}^{K} H_{k,m}^{T}X_{k,n}}$

ただし、 $\displaystyle X_{k,n} \stackrel{\mathrm{def}}{\equiv} \sum_{m=1}^{M} H_{k,m}U_{m,n}$

この更新式のよいところは、HおよびUの初期値としてすべて非負値を代入しておけば、何回更新しても非負値を保つという点にあります。これにより、非負値という制約を気にせずに最小化問題を解くことができます。

Hを更新→Uを更新→再びHを更新→・・・という逐次更新を数10～数100回程度繰り返す（Iteration）ことで、乖離度Dはおおむね極小値に収束します。

こうして乖離度Dを最小化するように得られたHとUは、音声スペクトログラムSを最も良く近似する行列因子であり、Hの各列ベクトルが音声中に存在する音色の周波数分布を表し、Uの各列ベクトルが各時刻における音量の大きさを表します。

なお、更新式はシグマ記号がたくさん並んでいるので難しく見えるかもしれませんが、実際の和計算は単なる行列の掛け算なので、プログラムに起こすとほんの数行で書けてしまいます。

NMFのバリエーション

上述したとおり、補助関数法を用いて上限となる関数を設計し、解析的に解けない乖離度Dの最小化問題をコンピュータで数値的に解くことができました。

2000年にLeeらがこの二乗誤差基準の乖離度Dに関する更新式を導いてから、NMFには用途に応じて数多くの改良が重ねられてきました。その方向性の一例が、乖離度Dの多様化です。具体的に以下のような乖離度Dのバリエーションが考案されています。

【一般化カルバック・ライブラー情報量（Kullback–Leibler divergence）】

$\displaystyle D_{KL}(S \| H \times U)=\sum_{k=1}^{K} \sum_{n=1}^{N} \left\{ S_{k,n} \log{\frac{S_{k,n}}{(H \times U)_{k,n}}}-S_{k,n}+(H \times U)_{k,n} \right\}$

【板倉・斎藤擬距離（Itakura–Saito divergence）】

$\displaystyle D_{IS}(S \| H \times U)=\sum_{k=1}^{K} \sum_{n=1}^{N} \left\{ \frac{S_{k,n}}{(H \times U)_{k,n}}-\log{\frac{S_{k,n}}{(H \times U)_{k,n}}}-1 \right\}$

それぞれの乖離度Dの定義に対して、二乗誤差のときと同様に偏微分を計算すれば、更新式を得ることができます。更新式自体はシンプルなので、プログラムに起こすことは簡単です。

また、３つの乖離度Dをパラメータβによって統一的に扱えるβダイバージェンス、それをさらに一般化したBregmanダイバージェンスと呼ばれる擬距離も定義可能ですが、式が複雑なのでここでは取り上げません。

なお、音楽の楽器パート分解や自動採譜（音楽からドレミ音階を自動抽出する）などの音声処理用途では、振幅スペクトル（二乗しない）Sに対して上述したKLダイバージェンスを最小化するのが最もよいパフォーマンスを与えると経験的に知られているようです*1。

Pythonで実装

更新式が数学的に導けているので、プログラムに起こすことは簡単です。Numpyライブラリによる行列計算が中心となります。ソースコードは以下のとおりです。

import sys
import soundfile as sf
import IPython.display as ipd
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
import numpy as np
from scipy import signal as sg
from scipy import fftpack as fp

#Function for converting images into animation
def update_anim(i, img):
    if i != 0:
        plt.cla()
    plt.pcolormesh(times, freqs, img[i], cmap = 'rainbow')
    plt.title("Iteration " + str(i+1) + ": the spectrogram approximated by NMF")
    plt.xlabel('Time [sec]')
    plt.ylabel('Frequency [Hz]')

#Function for removing components closing to zero
def get_nonzero(mat1, mat2):
    new_mat1 = np.where(mat1 < 1e-10, 1e-10, mat1)
    new_mat2 = np.where(mat2 < 1e-10, 1e-10, mat2)
    product = new_mat1 @ new_mat2
    return new_mat1, new_mat2, product

#Function for getting basements and activation matrix by NMF
def get_NMF(Y, num_iter, num_base, loss_func, save_anime):
    
    #Prepare the drawing area for animation
    plt.rcParams["font.size"] = 16
    if save_anime == True:
        fig = plt.figure(figsize=(12, 8))
    
    #Remove components closing to zero in Y
    Y = np.where(Y < 1e-10, 1e-10, Y)
    
    #Initialize basements and activation based on the Y size(k, n)
    K, N = Y.shape[0], Y.shape[1]
    if num_base >= K or num_base >= N:
        print("The number of basements should be lower than input size.")
        sys.exit()
    H = np.random.rand(K, num_base) #basements
    H = H / H.sum(axis=0, keepdims=True) #Normalization
    U = np.random.rand(num_base, N) #activation
    
    #Initialize valuables
    loss = np.zeros(num_iter)
    HU = []
    
    #Repeat num_iter times
    for i in range(num_iter):
        
        #In the case of squared Euclidean distance
        if loss_func == "EU":
            #Update the basements
            H, U, X = get_nonzero(H, U)
            H = H * (Y @ U.T) / (X @ U.T)
            H = H / H.sum(axis=0, keepdims=True) #Normalization
            
            #Update the activation
            H, U, X = get_nonzero(H, U)
            U = U * (H.T @ Y) / (H.T @ X)
            
            #Compute the loss function
            H, U, X = get_nonzero(H, U)
            loss[i] = np.sum((Y - X)**2)
        
        #In the case of Kullback–Leibler divergence
        elif loss_func == "KL":
            #Update the basements
            H, U, X = get_nonzero(H, U)
            denom_H = U.T.sum(axis=0, keepdims=True) #(1xM) matrix
            H = H * ((Y / X) @ U.T) / denom_H
            H = H / H.sum(axis=0, keepdims=True) #Normalization
            
            #Update the activation
            H, U, X = get_nonzero(H, U)
            denom_U = H.T.sum(axis=1, keepdims=True) #(Mx1) matrix
            U = U * (H.T @ (Y / X)) / denom_U
            
            #Compute the loss function
            H, U, X = get_nonzero(H, U)
            loss[i] = np.sum(Y * np.log(Y / X) - Y + X)
        
        #In the case of Itakura–Saito divergence
        elif loss_func == "IS":
            #Update the basements
            H, U, X = get_nonzero(H, U)
            denom_H = np.sqrt(X**-1 @ U.T)
            H = H * np.sqrt((Y / X**2) @ U.T) / denom_H
            H = H / H.sum(axis=0, keepdims=True) #Normalization
            
            #Update the activation
            H, U, X = get_nonzero(H, U)
            denom_U = np.sqrt(H.T @ X**-1)
            U = U * np.sqrt(H.T @ (Y / X**2)) / denom_U
            
            #Compute the loss function
            H, U, X = get_nonzero(H, U)
            loss[i] = np.sum(Y / X - np.log(Y / X) - 1)
        
        else:
            print("The deviation shold be either 'EU', 'KL', or 'IS'.")
            sys.exit()
        
        #Collecting reconstructed image for each iteration
        if save_anime == True:
            HU.append(H @ U)
    
    #Save the images as animation
    if save_anime == True:
        #Pass the H @ U list into "update_anim" and repeat "frames" times
        anime = FuncAnimation(fig, update_anim, fargs=[HU], interval=200, frames=num_iter)
        anime.save('NMF_spectrogram.gif', writer=PillowWriter(fps=5))
        plt.close()
    
    return H, U, loss

#■Main■
if __name__ == "__main__":
    
    #Setup
    audiolen = None        #Cropping time (second) [Default]None(=without cropping)
    frame_length = 0.04    #STFT window width (second) [Default]0.04
    frame_shift = 0.02     #STFT window shift (second) [Default]0.02
    num_iter = 100         #The number of iteration [Default]100
    num_base = 3           #The number of basements [Default]3
    spec_type = "amp"      #The type of spectrum (amp: amplitude, pow: power) [Default]amp
    loss_func = "KL"       #Select either EU, KL, or IS divergence [Default]KL
    save_anime = False     #Save the reconstructed spectrogram for each iteration as anime
    
    #Compute the Mel-scale spectrogram
    data, Fs = sf.read("./piano.wav")
    
    #Transform stereo into monoral
    if data.ndim == 2:
        wavdata = 0.5*data[:, 0] + 0.5*data[:, 1]
    else:
        wavdata = data
    
    #Down sample and normalize the wave
    wavdata = sg.resample_poly(wavdata, 8000, Fs)
    Fs = 8000
    
    #Original sound
    print("Original sound")
    ipd.display(ipd.Audio(data=wavdata, rate=Fs))
    
    #Calculate the index of window size and overlap
    FL = round(frame_length * Fs)
    FS = round(frame_shift * Fs)
    OL = FL - FS
    
    #Execute STFT
    freqs, times, dft = sg.stft(wavdata, fs=Fs, window='hamm', nperseg=FL, noverlap=OL)
    arg = np.angle(dft) #Preserve the phase
    Adft = np.abs(dft) #Preserve the absolute amplitude
    if spec_type == "amp":
        Y = Adft
    elif spec_type == "pow":
        Y = Adft**2
    else:
        print("The spec_type must be either 'amp' or 'pow'.")
        sys.exit()
    
    #Display the size of input
    print("Spectrogram size (freq, time) = " + str(Y.shape))
    print("Basements for NMF = " + str(num_base) + "\n")
    
    #Call my function for updating NMF basements and activation
    H, U, loss = get_NMF(Y, num_iter, num_base, loss_func, save_anime)
    
    #Pick up each basement vector
    for j in range(num_base):
        zero_mat = np.zeros((H.shape[1], H.shape[1]))
        zero_mat[j, j] = 1
        Base = (H @ zero_mat) @ U #Extract an only dictionary
        Base = Base * np.exp(1j*arg) #Restrive the phase from original wave
        _, Base_wav = sg.istft(Base, fs=Fs, window='hamm', nperseg=FL, noverlap=OL)
        print("Basement " + str(j+1))
        ipd.display(ipd.Audio(data=Base_wav, rate=Fs))
        if j ==0:
            Sum_wav = Base_wav
        else:
            Sum_wav = Sum_wav + Base_wav
    print("Sum all basements (approximation by NMF)")
    ipd.display(ipd.Audio(data=Sum_wav, rate=Fs))
    
    #Plot the loss curve
    plt.figure(figsize=(18, 6))
    plt.subplot(1, 2, 1)
    plt.title('An original spectrogram')
    plt.xlabel('Time [sec]')
    plt.ylabel('Frequency [Hz]')
    plt.pcolormesh(times, freqs, Y, cmap = 'jet')
    plt.colorbar(orientation='horizontal').set_label('Power')
    
    plt.subplot(1, 2, 2)
    plt.title('The spectrogram approximated by NMF')
    plt.xlabel('Time [sec]')
    plt.ylabel('Frequency [Hz]')
    plt.pcolormesh(times, freqs, (H @ U), cmap = 'jet')
    plt.colorbar(orientation='horizontal').set_label('Power')
    
    #Plot the loss curve
    plt.figure(figsize=(10, 5))
    plt.plot(np.arange(num_iter)[:], loss[:], marker='.')
    plt.title(loss_func + '_loss curve')
    plt.xlabel('Iteration')
    plt.ylabel('Loss value')

まず注意点として、ところどころにある「ipd.display(ipd.Audio()」というのは、IPython.displayというライブラリを使ってjupyter notebook上に音楽の再生バーを表示する処理なので、それ以外のエディタで実行する場合には動きません。そういう方は「soundfile.write()」などを使って音楽データ自体をwav形式で保存してください。

メインのNMFの行列計算は「get_NMF」で行われています。それぞれの乖離度Dに対して場合分けをしているのでちょっと長くなっていますが計算はシンプルです。ただし、更新式の分母がゼロになると割り算でエラーが出るので、それを防ぐために「get_nonzero」という関数を作って行列要素がゼロにならないように処理しています*2。

また、更新を繰り返すたびにNMFの近似精度が上がっていく様子を可視化してみたかったので、「update_anim」という関数でmatplotlibライブラリのanimationクラスを呼び出してgif形式のアニメーションを作りました。以下がその様子です。

f:id:yuki0718:20200216000424g:plain

KLダイバージェンス基準で更新していますが、かなり収束は速いことが伺えます。目視ではものの20回くらいの繰り返し回数（Iteration）でほぼ収束しているように見えます。

続いて、実際に以下のようなピアノ音源をサンプル入力にしてNMFで分解した結果を示します。このサンプルは、初めにピアノ音が鳴り、しばらくして別の高さのピアノ音が鳴り、最後に2つのピアノ音が同時に鳴る、という10秒強の音源です。

音源がかなり単純なので、基底ベクトルの本数（基底行列Hの列数）Mは3本に設定して分解しました。

このコードではNMFで得られた基底ベクトルが何を表しているのか確認するのが目的なので、「#Pick up each basement vector」のところで、基底行列Hのうち特定列の基底ベクトル以外を0にした行列を順次作り、右からUを掛け、逆STFTでスペクトログラムからwav音楽データに戻しています。

Basement1

Basement2

Basement3

どうやら、一つ目の基底ベクトルは最初の音だけを拾っており、二つ目の基底ベクトルは別の高さの音だけを拾っているようです。二つの音が同時に鳴っているところでも、きちんと分離抽出できています。三つ目の基底ベクトルは鍵盤を叩いた時の打音だと思われます。

Basement1+2+3

上述したBasement1～3の分離音源を足し合わせたものは、ようするに $H \times U$ を計算した結果（NMFによる近似スペクトログラム）です。実際に聴いてみると近似とは思えないくらい完璧に再現できていて、機械学習ってすごいなあと思いました。

今回のところは以上です。これまでの定式化とコードを見ていてお気づきかと思いますが、NMFはいわゆる教師なし機械学習（Unsupervised machine learning）です。ディープラーニングとは違って、一切の前提情報を与えずともいきなり本番一発勝負でデータを分離できます。これはなかなかすごいことです。

次回はこのNMFを教師データ「あり」の場合に拡張した、半教師ありNMF（SSNMF: Semi-supervised non-negative matrix factorization）についてご紹介する予定です。

稚文をお読みいただきありがとうございました。

*1:D. Fitzgerald et al., “On the use of the Beta Divergence for Musical Source Separation”, Proc. Irish Signals Syst. Conf., 2009（https://ieeexplore.ieee.org/document/5524688）

*2:他にもっといいやり方があるんでしょうが、特にこれで精度がガタ落ちしたりすることはないので安直に小さい値1e-10を加えるという解決策を採用しました。