keep learning blog(キープラーニングブログ)

自分が興味を持ったことを備忘録として残すブログです。

17.混合ガウスモデルでクラスタリングしよう

――私が述べるすべての文章は、断定ではなく、質問であると理解されるべきである。――

ニールス・ヘンリク・ダヴィド・ボーア

 

 

混合ガウスモデルでクラスタリングしよう

冒頭の言葉は、量子力学の祖にしてアインシュタインの永遠のライバル、ニールス・ボーアの名言です。ボーア・アインシュタイン論争という有名な科学史の逸話があります。

実は、ボーアはコペンハーゲン大学出身のデンマーク人です。彼の残した言葉は多くが論文という形でしか残っていないため、数少ない名言の一つではありますが、この言葉から彼の研究者としての謙虚さが伺えます。というわけで、このブログの内容も断定ではなく質問として受け取っていただければ幸いです。

f:id:yuki0718:20191119034942j:plain

さて、前回の記事では、確率論で定式化される古典的な機械学習モデルとして、混合ガウスモデル(GMM:Gaussian Mixture Model)をご紹介し、最尤推定法によりパラメータ \theta を求めるための公式を導きました。今回はそれを使ってプログラムを実装し、実際にクラスタリングを行ってみます。

keep-learning.hatenablog.jp

 

混合ガウスモデル(GMM)まとめ

GMMは生成モデルの一種であり、その形状を規定するパラメータ \theta としては重み係数 w_{k} と平均ベクトル \boldsymbol{\mu}_{k} と共分散行列 \boldsymbol{\Sigma}_{k} の3種類が存在します。これらを計算するための公式は最尤推定法により、

 \displaystyle w_{k}=\frac{\sum_{i=1}^{N} \gamma_{i,k}}{N}

 \displaystyle \boldsymbol{\mu}_{k}=\frac{\sum_{i=1}^{N} \gamma_{i,k}\boldsymbol{x}_{i} }{\sum_{i=1}^{N} \gamma_{i,k}}

 \displaystyle \boldsymbol{\Sigma}_{k}=\frac{\sum_{i=1}^{N} \gamma_{i,k}(\boldsymbol{x}_{i}-\boldsymbol{\mu}_{k})(\boldsymbol{x}_{i}-\boldsymbol{\mu}_{k})^{T} }{\sum_{i=1}^{N} \gamma_{i,k}}

ただし、負担率 \displaystyle \gamma_{i,k} \stackrel{\mathrm{def}}{\equiv} \frac{w_{k}Norm(\boldsymbol{x}_{i} \| \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})}{\sum_{k=1}^{C} w_{k}Norm(\boldsymbol{x}_{i} \| \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})}

とまとめられるのでした。

 

自己無撞着(Self-Consistent)方程式と逐次法

この公式は左辺と右辺の両方にパラメータを含む、いわゆる自己無撞着(Self-Consistent)方程式です。物理でいうと電磁場の平均場近似や多体系量子力学のハートリー=フォック方程式などが自己無撞着方程式として有名です(ので私もよく知っています)。

自己無撞着方程式には多くの場合に解析解が存在しないものの、コンピュータで解く場合には逐次法と呼ばれる手法を用います。

逐次法といっても特別難しいことをするわけではありません。まず、パラメータ w_{k}, \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k} の初期値を何らかの方法で適当に決めます。そして、その初期値を上述した自己無撞着方程式の右辺に代入して新たなパラメータを得ます。新たなパラメータを初期値として再び右辺に代入して新たなパラメータを得ます。これを何回か繰り返して正しい解に近づけていくという手法です*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)」は混合ガウス分布の潜在変数 k について繰り返し、上述の関数「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。試験データは以下のような特徴量データの分布になっています。

f:id:yuki0718:20191104033003p:plain

今回は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

f:id:yuki0718:20191104033024p:plain

正規分布の平均ベクトルの位置を「●」で表し、それを囲むように共分散行列に基づいた等高線を表示しています。正規分布の平均ベクトルに近い位置ほど(概ね)負担率が高いことを表します*4

なんとなくこの辺りに塊があるな、と人間が見て理解できる位置にきちんと正規分布が重なっており、きちんとクラスタリングできていることがグラフィカルに理解できると思います。

 

以上、今回は前回の記事に書ききれなかったプログラムの実装部分でした。機械学習という手法は様々なデータ分析に役立つ道具だと少しでも伝わっていたら幸いです。

特徴量データとして私のように音声を用いてもいいですし、ワインの特徴量を用いても、画像データを用いてもいいです。入力するデータから何を得たいのか、目的次第で適切な機械学習モデルの選択が可能です。

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

*1:第8回の記事でHarmonic/Percussive Sound Separation(HPSS)という音声分離アルゴリズムを実装したときにも、逐次法によって解を収束させていました。

*2:「現代はscikit-learnのように便利なライブラリがあるんだから古典のGMMなんて勉強しなくていいじゃないか」と思うかもしれません。しかし、今回のようにちょっとしたアレンジ(時間変化する音声データを受け入れる)を加えたいときに、基礎理論が分かっていないと手も足も出ません。

*3:音声データの分析では音声データをMel-Frequency Cepstrum Coefficient(MFCC)という係数に変換して特徴量とする手法が有名です。

*4:負担率の分子は多変量正規分布に重み係数(混合係数)wを掛けた値なので、本当は重み係数wの分だけグラフとズレます。つまり、厳密にはデータ毎に4種類の正規分布の負担率を計算し、最も負担率が大きい正規分布のグループに属すると判定するのが通常の手続きです。あくまでイメージを掴んでもらうための概略図だとご理解ください。