17.混合ガウスモデルでクラスタリングしよう
――私が述べるすべての文章は、断定ではなく、質問であると理解されるべきである。――
ニールス・ヘンリク・ダヴィド・ボーア
混合ガウスモデルでクラスタリングしよう
冒頭の言葉は、量子力学の祖にしてアインシュタインの永遠のライバル、ニールス・ボーアの名言です。ボーア・アインシュタイン論争という有名な科学史の逸話があります。
実は、ボーアはコペンハーゲン大学出身のデンマーク人です。彼の残した言葉は多くが論文という形でしか残っていないため、数少ない名言の一つではありますが、この言葉から彼の研究者としての謙虚さが伺えます。というわけで、このブログの内容も断定ではなく質問として受け取っていただければ幸いです。
さて、前回の記事では、確率論で定式化される古典的な機械学習モデルとして、混合ガウスモデル(GMM:Gaussian Mixture Model)をご紹介し、最尤推定法によりパラメータを求めるための公式を導きました。今回はそれを使ってプログラムを実装し、実際にクラスタリングを行ってみます。
混合ガウスモデル(GMM)まとめ
GMMは生成モデルの一種であり、その形状を規定するパラメータとしては重み係数と平均ベクトルと共分散行列の3種類が存在します。これらを計算するための公式は最尤推定法により、
ただし、負担率
とまとめられるのでした。
自己無撞着(Self-Consistent)方程式と逐次法
この公式は左辺と右辺の両方にパラメータを含む、いわゆる自己無撞着(Self-Consistent)方程式です。物理でいうと電磁場の平均場近似や多体系量子力学のハートリー=フォック方程式などが自己無撞着方程式として有名です(ので私もよく知っています)。
自己無撞着方程式には多くの場合に解析解が存在しないものの、コンピュータで解く場合には逐次法と呼ばれる手法を用います。
逐次法といっても特別難しいことをするわけではありません。まず、パラメータの初期値を何らかの方法で適当に決めます。そして、その初期値を上述した自己無撞着方程式の右辺に代入して新たなパラメータを得ます。新たなパラメータを初期値として再び右辺に代入して新たなパラメータを得ます。これを何回か繰り返して正しい解に近づけていくという手法です*1。
実は、GMMのパラメータを逐次法で求める計算過程は、機械学習分野で広く用いられているEMアルゴリズム(Estimation-Maximization Algorithm)と呼ばれる手法と全く同じであることが知られています。
EMアルゴリズムは、Estimation StepとMaximization Stepと呼ばれる更新を繰り返すたびに、常に解に向かって収束していくことが数学的に証明されています。今は脇道に逸れるので書きませんが、そのうちEMアルゴリズムのことも紹介できたらと思っています。
Pythonで実装
実装したのでソースコード載せます。まずは多変量正規分布をnumpyの行列演算で計算するためのサブ関数を紹介します。
#■Function for calculating the gaussian■ def get_gaussprob(data, w, m, cov): #data: matrix(nframes x ndim) #w: scalar, m: vector(ndim), v: matrix(ndim x ndim) #Get the number of frames and features nframes = data.shape[0] ndim = data.shape[1] #Duplicate mean m = m[np.newaxis, :] m = np.tile(m, (nframes, 1)) #Avoid the underflow of determinant cov_det = abs(np.linalg.det(cov)) uf = False if cov_det < 1e-150: cov = cov + 1e-6 * np.identity(ndim) cov_det = abs(np.linalg.det(cov)) uf = True #Calculate the logarithm of gaussian with weight C = ndim * np.log(2*np.pi) + np.log(cov_det) EXP = (data - m) @ (np.linalg.inv(cov)) EXP = EXP @ (data - m).T EXP = np.diag(EXP) gaussprob = -0.5*(EXP + C) + np.log(w) #Convert into exponential gaussprob = np.exp(gaussprob) #Return the wk*Norm(x|m,cov) as time-frames vector return gaussprob, uf
引数としてdataと3種類のパラメータw, m, covを入力すると、多変量正規分布を計算して値を返すだけのシンプルな関数です。ufは学習が進んで行列式が64bit浮動小数点で表現できなくなりそうというフラグ(アンダーフロー防止)で、通常はFalseが返ってきます。
ただし、ピッチ推定でもそうだったように、音声データというのは時々刻々と特徴量が変化するデータです(ワインの特徴量のようにある瞬間のデータではない)。そのため、教科書的なGMMがdataとして特徴量次元数のベクトルを入力するのに対し、このプログラムではdataとして時間フレーム数と特徴量次元数の行列を入力しています*2。
続いて、GMMの逐次法によるパラメータ計算を行うサブ関数です。
#■Function for training the GMM■ def GMM_train(train_x, nmix, max_iter, threshold, ini_mode): #Check the parameter if train_x[0].ndim != 2: print("Each data should be 2-dimension tensor(time,feature).") sys.exit() #Get the number of audio files and features nfiles = len(train_x) ndim = train_x[0].shape[1] #Define the covariance delta to avoid underflow delta = np.zeros((nmix, ndim, ndim)) #Call my function for initializing GMM parameters [w, m, cov] = ini_param(train_x, nmix, ini_mode) #Repeat every iteration score = 0 i = 0 while i < max_iter: #Get the start time start = time.time() #Initialize variable total_frames = 0 logL = 0 N = np.zeros(nmix) F = np.zeros((nmix, ndim)) S = np.zeros((nmix, ndim, ndim)) #Repeat every audio file for d in range(nfiles): #Get one data (matrix: nframes x ndim) data = train_x[d] #Total frames is summation of nframe in all dataset nframes = data.shape[0] #nframes is diffrent from audio to audio total_frames = total_frames + nframes #Initialize the variable G = np.zeros((nmix, nframes)) L = np.zeros(nframes) gamma = np.zeros((nmix, nframes)) #gannma: matrix (nmix x nframes) #Repeat every Gaussian for mix in range(nmix): #Call my function for calculating the Gaussian G[mix], uf = get_gaussprob(data, w[mix], m[mix, :], cov[mix, :, :]) #Avoid underflow of covariance if uf == True: cov[mix] = cov[mix] + 1e-6 * np.identity(ndim) delta[mix] = 1e-6 * np.identity(ndim) #Get the summation of Gaussian (Likelihood for one data) L = L + G[mix] #Avoid log(0) error if np.any(L == 0): #print("Data No." + str(d) + ": add 1e-6 to likelihood to avoid log(0) error.") L = np.where(L == 0, 1e-6, L) #Get the summation of log-scale likelihood logL = logL + np.sum(np.log(L)) #logL: scalar #Calcurate the responsibility(γ) L = L[np.newaxis, :] L = np.tile(L, (nmix, 1)) gamma = G / L ### Expectation step ### #data: matrix (nframes x ndim), gamma: matrix (nmix x nframes) N = N + np.sum(gamma, axis=1) #Get the summation along with time axis F = F + gamma @ data #Calculate the multiplication of gamma by data S = S + sum_outer_product(gamma, data) #Calculate the summation of outer product #Avoid the zero division error if np.any(N == 0): #print("Add 1e-6 to N to avoid zero division error.") N = np.where(N == 0, 1e-6, N) ### Maximization Step ### w = N / total_frames #w: vector (nmix) N2 = N[:, np.newaxis] N2 = np.tile(N2, (1, ndim)) #Duplicate in order to Hadamard operation m = F / N2 #m: matrix (nmix x ndim) N3 = N2[:, :, np.newaxis] N3 = np.tile(N3, (1, 1, ndim)) #Duplicate in order to Hadamard operation Imix = np.identity(nmix) #The function "sum_outer_product" results in outer product by using Identity matrix cov = S / N3 - sum_outer_product(Imix, m) + delta #Unless the underflow happens, the delta = 0. #Calculate the total likelihood (score) Likelihood = logL / total_frames #Get summation and divide by N,T if i == 0: diff = 0 else: diff = Likelihood - score #diffrence from the former iteration score = Likelihood #Output the result and process time finish = time.time() - start print("N_mix={}, Iter{}, Likelihood={:.1f}, Gain={:.4f}, Process_time={:.1f}sec".format(nmix, i+1, score, diff, finish)) #Condition of convergence if i != 0 and diff < threshold: break i = i + 1 #Return the Gaussian Parameter for GMM return w, m, cov
「while」のところで、最大繰り返し回数(Max iteration)まで逐次法によるパラメータ更新を繰り返しています。「for d in range(nfiles)」はデータを1つずつ処理しており、「for mix in range(nmix)」は混合ガウス分布の潜在変数について繰り返し、上述の関数「get_gaussprob」を呼び出して多変量正規分布を計算しています。
時間フレームが加わっている分テンソルの次元が1つ増えるので、変数に入っているのがどういう形のテンソルなのか常にコメントに書いています。こうしないと混乱するので、正しく行列演算をプログラミングできません。
なお、この関数中では「sum_outer_product」と「ini_param」という自己定義関数を呼び出しています。これらはΣ和の計算を効率的に行ったり、逐次法の最初でパラメータの初期値を決定したりするためのものです。完全オリジナルで組んでおり、間違っていたら恥ずかしいので、ここには掲載しません(一応ちゃんと動くので正しいと信じています)。
最後に、データをサブ関数に受け渡して結果を表示するためのメイン関数を紹介します。
本来は音声データを加工して特徴量データに変換する処理やデータクレンジングと呼ばれる前処理を書きますが、音声データの加工方法を説明し始めるとGMMから話が脱線してしまうので、今回は乱数を使ってクラスタリング用の試験データを人工的に作ってみました。
import sys import math import time import matplotlib.pyplot as plt import numpy as np from sklearn.cluster import KMeans from scipy.stats import multivariate_normal #■Main script■ if __name__ == "__main__": #Set up GMM_ini = "kmeans" #How to initialize the Gaussian parameters(kmeans++, kmeans, random) niter = 30 #The number of iterations nmix = 4 #The number of Gaussian mixtures threshold = 1e-3 #The likelihood threshold for stopping iteration #Generate demo data N = 500 data = np.concatenate([np.random.multivariate_normal([-2, 0], np.eye(2), round(N/3)), np.random.multivariate_normal([0, 5], np.eye(2), round(N/3)), np.random.multivariate_normal([2, -2], np.eye(2), round(N/3)), np.random.multivariate_normal([4, 3], np.eye(2), round(N/3))]) data = data[:, np.newaxis, :] data = np.tile(data, (1, 100, 1)) #Call my function for training the GMM [w, m, cov] = GMM_train(data, nmix, niter, threshold, GMM_ini) #Plot the Gaussian distribution plt.rcParams["font.size"] = 16 plt.figure(figsize=(12, 8)) colors = np.repeat(['r', 'g', 'b', 'm'], np.ceil(nmix/4)) for k in range(nmix): plt.scatter(m[k, 0], m[k, 1], c=colors[k], marker='o', zorder=3) #Plot the contour graph nfiles = len(data) ndim = data[0].shape[1] ave_x = np.zeros((nfiles, ndim)) for d in range(nfiles): ave_x[d, :] = np.average(data[d], axis=0) x_min, x_max = np.amin(ave_x[:, 0])-0.5, np.amax(ave_x[:, 0])+0.5 y_min, y_max = np.amin(ave_x[:, 1])-0.5, np.amax(ave_x[:, 1])+0.5 xlist = np.linspace(x_min, x_max, 50) ylist = np.linspace(y_min, y_max, 50) x, y = np.meshgrid(xlist, ylist) pos = np.dstack((x,y)) for k in range(nmix): z = multivariate_normal(m[k, 0:2], cov[k, 0:2, 0:2]).pdf(pos) cs = plt.contour(x, y, z, 3, colors=colors[k], linewidths=2, zorder=2) #Plot the time-averaged training data plt.plot(ave_x[:, 0], ave_x[:, 1], 'cx', zorder=1) plt.xlim(x_min, x_max) plt.ylim(y_min, y_max) plt.title('GMM distribution') plt.xlabel('Feature1') plt.ylabel('Feature2') plt.show()
「#Generate demo data」のところで2次元の特徴量からなる試験データを人工的に作っています。2次元にしたのは、グラフに図示できるのが2次元までだからです。音声の特徴量データは40次元くらい余裕で超えたりします*3。試験データは以下のような特徴量データの分布になっています。
今回は4種類くらいの塊に分かれるように、乱数で500個の特徴量データを作ってみました。なお、上述したとおり、このプログラムは時間変化するデータを受け入れ可能ですが、試験データは時間変化なし、つまりずっとこの分布のまま一定にしました(面倒だったので)。
なんとなく私が見ても4種類くらいにグループ分けできそうです。これを混合数4「nmix = 4」の混合ガウス分布でクラスタリングした場合に、4つの多変量正規分布がうまい具合に4種類の塊を認識してくれれば成功です。実際に上述のプログラムを動かしてみた結果、以下のようになりました。
N_mix=4, Iter1, Likelihood=-5.0, Gain=0.0000, Process_time=0.6sec
N_mix=4, Iter2, Likelihood=-4.7, Gain=0.2808, Process_time=0.6sec
N_mix=4, Iter3, Likelihood=-4.5, Gain=0.1635, Process_time=0.6sec
N_mix=4, Iter4, Likelihood=-4.3, Gain=0.2332, Process_time=0.6sec
N_mix=4, Iter5, Likelihood=-4.2, Gain=0.1073, Process_time=0.6sec
N_mix=4, Iter6, Likelihood=-4.2, Gain=0.0216, Process_time=0.6sec
N_mix=4, Iter7, Likelihood=-4.2, Gain=0.0048, Process_time=0.6sec
N_mix=4, Iter8, Likelihood=-4.2, Gain=0.0022, Process_time=0.6sec
N_mix=4, Iter9, Likelihood=-4.2, Gain=0.0011, Process_time=0.6sec
N_mix=4, Iter10, Likelihood=-4.2, Gain=0.0006, Process_time=0.6sec
正規分布の平均ベクトルの位置を「●」で表し、それを囲むように共分散行列に基づいた等高線を表示しています。正規分布の平均ベクトルに近い位置ほど(概ね)負担率が高いことを表します*4。
なんとなくこの辺りに塊があるな、と人間が見て理解できる位置にきちんと正規分布が重なっており、きちんとクラスタリングできていることがグラフィカルに理解できると思います。
以上、今回は前回の記事に書ききれなかったプログラムの実装部分でした。機械学習という手法は様々なデータ分析に役立つ道具だと少しでも伝わっていたら幸いです。
特徴量データとして私のように音声を用いてもいいですし、ワインの特徴量を用いても、画像データを用いてもいいです。入力するデータから何を得たいのか、目的次第で適切な機械学習モデルの選択が可能です。
稚文をお読みいただきありがとうございました。
*1:第8回の記事でHarmonic/Percussive Sound Separation(HPSS)という音声分離アルゴリズムを実装したときにも、逐次法によって解を収束させていました。
*2:「現代はscikit-learnのように便利なライブラリがあるんだから古典のGMMなんて勉強しなくていいじゃないか」と思うかもしれません。しかし、今回のようにちょっとしたアレンジ(時間変化する音声データを受け入れる)を加えたいときに、基礎理論が分かっていないと手も足も出ません。
*3:音声データの分析では音声データをMel-Frequency Cepstrum Coefficient(MFCC)という係数に変換して特徴量とする手法が有名です。
*4:負担率の分子は多変量正規分布に重み係数(混合係数)wを掛けた値なので、本当は重み係数wの分だけグラフとズレます。つまり、厳密にはデータ毎に4種類の正規分布の負担率を計算し、最も負担率が大きい正規分布のグループに属すると判定するのが通常の手続きです。あくまでイメージを掴んでもらうための概略図だとご理解ください。