SHIROBAKO大好き人間のブログ

SHIROBAKOが好きなエンジニアによる技術ブログ

Adaptive Restricted Boltzmann Machine (ARBM) の実装

前回、音源分離の記事を書いた後、音声処理に関する色々な記事を読んでいたら面白そうな論文がありました。

話者適応型Restricted Boltzmann Machineを用いた声質変換の検討(中鹿ら, 2014)

声質変換とは、入力した話者の音声を別の話者が喋っているような音声に変換することです。

このタイプの研究は以前からあったようですが、今までの研究ではパラレルデータが必要でした。 パラレルデータとは話者Aと話者Bが同じ内容を喋っているような音声データのことです。 そういうデータが用意できれば問題ありませんが、現実的には用意するのは難しいです。

要するに、手法は色々あるけど、使えるデータが限られていました。

この論文では、Adaptive Restricted Boltzmann Machine (ARBM) というモデルを用いて、パラレルデータがなくても声質変換ができるような手法を提案しています。 パラレルデータがなくてもいいなら自分の好きなデータで声質変換が試せそうなので、試してみることにしました。 そのために、まずは論文中で使われているARBMを実装してみました。 実際に声質変換をするのは次回の記事でやります。

そもそもRBMって?

ARBMについて説明する前に、その基本であるRBMについてざっくり説明します。

RBM (Restricted Boltzmann Machine) とは確率モデルの一種で、以下のようなモデルです。

f:id:phoro3:20161119213739p:plain

(図はこちらのブログより引用)


可視層と隠れ層の2つの層があり、可視層と隠れ層の間にはリンクがありますが、可視層同士、および隠れ層同士の間にはリンクがありません。

また、可視層と隠れ層は{0,1}のどちらかの値しかとらないとします。

イメージとしては、可視層に何らかのデータが入力されたらその値が重み付けされて隠れ層に伝わるという感じです。


このモデルにおいて、次のような確率分布p({\bf v}, {\bf h})を考えます。

p({\bf v}, {\bf h})=\frac{1}{Z} \exp(-\Phi({\bf v}, {\bf h}))

\displaystyle \Phi({\bf v}, {\bf h}) = -{\bf b}^{T}{\bf v} - {\bf c}^T {\bf h} - {\bf v}^T {\bf W} {\bf h}

\displaystyle Z = \sum_{{\bf v}, {\bf h}} \exp(-\Phi({\bf v}, {\bf h}))

{\bf v}=(v_1,v_2,...,v_n)^{T}は可視層の値を表すベクトル、{\bf h}=(h_1,h_2,...,h_m)^{T}は隠れ層の値を表すベクトルです。

また、 {\bf b}, {\bf c}, {\bf W}はパラメータで、学習を行いながらこの値を変化させます。

{\bf b}のサイズはn \times 1{\bf c}m \times 1{\bf W}n \times mです。

観測データ{\bf v}に対してp({\bf v}, {\bf h})の尤度が最もが高くなるようなパラメータを求めるのがRBMの学習の目的です。

観測データの尤度が高くなるように学習するので、上手く学習できれば、RBMの隠れ層は入力の隠れた特徴を表した層になっていることが期待できます。




RBMの説明はこれくらいにしておきます。

今は、「可視層と隠れ層がある確率モデル」くらいのことがわかってもらえれば大丈夫です。

RBMに関する詳しいことは、ネットや書籍の方がよくまとまっているのでそちらを参考にしてください。

以下のブログは、理論についても実装についてもうまくまとまっているのでオススメです。

ARBMの概要

では、本題であるARBMを見ていきましょう。

必要なデータ

ARBMを使った声質変換では以下のデータを用います。

  • 参照話者の音声データ

モデルを訓練する為に変換したい話者とは別の話者の音声データを用います。

これを参照話者と呼び、R_1,R_2,...,R_Sとします。

  • 変換元・変換先話者の音声データ

変換したい話者の特徴をモデルに学習させる為に、変換元と変換先の話者の音声データも必要です。

ここでは、変換元話者を R_{source}、変換先話者を R_{target}とします。

ARBMのモデル

ARBMのモデルを表したのが以下の図です。

f:id:phoro3:20161120162711p:plain

(図はこちらの論文の3ページより引用)

論文からの引用なので、先ほどのRBMの図と微妙に見た目は違いますが、下の黒いノードが可視層、上のグレーのノードが隠れ層です。

RBMとの違いは可視層と隠れ層の間にsという識別層が加わっていることです。

また、先ほどのRBMの例では可視層も隠れ層も0,1の値しか取りませんでしたが、ここでは可視層は実数の値を取ることとします。



このモデルにおいて次のような確率分布を考えます。

p({\bf v}, {\bf h}, s)=\frac{1}{Z} \exp(-\Phi({\bf v}, {\bf h}, s))

\displaystyle \Phi({\bf v}, {\bf h}, s) = \left| \left| \frac{{\bf v} - {\bf b}} {2{\boldsymbol \sigma}} \right| \right|^2 - {\bf c}^T {\bf h} - \left( \frac{{\bf v}}{{\boldsymbol \sigma}^2} \right)^{T} {\bf W}_s {\bf h}

 {\bf W}_s = {\bf A}_s \bar{\bf W} + {\bf B}_s

\displaystyle Z = \sum_{{\bf v},{\bf h},s} \exp(-\Phi({\bf v}, {\bf h}, s))

各行列およびベクトルのサイズはこうなってます。

{\bf v},{\bf b},{\boldsymbol \sigma}:n \times 1

{\bf h},{\bf c}:m \times 1

{\bf W}_s,\bar{{\bf W}},{\bf B}_s:n \times m

{\bf A}_s:n \times n

先ほどと同様、{\bf v}が可視層のベクトル、{\bf h}が隠れ層のベクトルです。

\Phiの中身がRBMとはだいぶ変わっていますね。

\left| \left| \frac{{\bf v} - {\bf b}} {2{\boldsymbol \sigma}} \right| \right|^2の部分は、RBMにおける{\bf b}^T {\bf v}の部分を連続値に対応できるように変えたものです。

また、RBMでは{\bf W}だった重みが{\bf W}_sになっています。

ここがポイントで、ARBMでは話者ごとに可視層と隠れ層の間の重みを変化させます。

{\bf W}_sの添え字sは話者を表しており、話者R_sのデータが可視層に入力された時は重み{\bf W}_sを用います。

重み{\bf W}_sを詳しく見てみると、{\bf A}_s,\bar{{\bf W}},{\bf B}_sが登場しています。

このうち、\bar{{\bf W}}は話者に依存しない重み、{\bf A}_s,{\bf B}_s\bar{{\bf W}}に対する係数およびバイアスを表します。

\bar{{\bf W}}は話者に依存しないので、人が声を発する際に共通する特徴を表すと考えられます。

R_1,R_2,...,R_Sの複数人のデータで訓練を行いながら、それらのデータに共通する\bar{{\bf W}}を求められるのがこのARBMの大きな特徴です。



ARBMのパラメータは{\bf b},{\bf c},{\boldsymbol \sigma},{\bf A}_s,\bar{{\bf W}},{\bf B}_sです。

これらのパラメータを、N個の観測データ{\bf v}に対して、対数尤度

\displaystyle L=\log \prod_{i=1}^{N} p({\bf v}_i, s)=\log \prod_{i=1}^{N} \sum_{{\bf h}} p({\bf v}, {\bf h}, s)

が最大となるように更新していきます。

論文によると、この対数尤度の微分は以下のようになります。

\displaystyle \frac{\partial L}{\partial b_i}= \langle \frac{v_i}{\sigma_i^2} \rangle_{data} - \langle \frac{v_i}{\sigma_i^2} \rangle_{model}

\displaystyle \frac{\partial L}{\partial c_j}= \langle h_j \rangle_{data} - \langle h_j \rangle_{model}

\displaystyle \frac{\partial L}{\partial (A_s)_{ii'}}= \langle \sum_m \frac{(W_s)_{i' m} v_i h_m}{\sigma_i^2} \rangle_{data} - \langle \sum_m \frac{(W_s)_{i' m} v_i h_m}{\sigma_i^2} \rangle_{model}

\displaystyle \frac{\partial L}{\partial (\bar{W}_s)_{ij}} = \langle \sum_{l, k} \frac{(A_k)_{li} v_l h_j}{\sigma_i^2} \rangle_{data} - \langle \sum_{l, k} \frac{(A_k)_{li} v_l h_j}{\sigma_i^2} \rangle_{model}

\displaystyle \frac{\partial L}{\partial (B_s)_{ij}}= \langle v_i h_j \rangle_{data} - \langle v_i h_j \rangle_{model}

式中の(A_s)_{ij}という表記は行列{\bf A}_s(i,j)成分を表します。



{\boldsymbol \sigma}に関しては学習を安定させるため、\sigma_i^2 = e^{z_i}とし、z_iを以下のように更新します。

\displaystyle \frac{\partial L}{\partial z_i}= e^{-z_i} \langle \frac{(v_i - b_i)^2}{2}- v_i (W_s)_{i:} {\bf h} \rangle_{data}

\hspace{30pt} - e^{-z_i} \langle \frac{(v_i - b_i)^2}{2} - v_i (W_s)_{i:} {\bf h} \rangle_{model}

ただし、(W_s)_{i:}は行列{\bf W}_si行目の成分を表します。



更新式で出てくる\langle \cdot \rangle_{data}は観測データの期待値を表します。

つまり、\langle v_i \rangle_{data}は入力したデータのi番目の成分の期待値を取ればOKです。

では、\langle h_j \rangle_{data}はどうすればいいでしょうか?

実は、p({\bf v}, {\bf h}, s)の定義から、入力{\bf v}sが与えられた時に、h_j=1となる確率p(h_j = 1 | {\bf v}, s)が計算できます。

p(h_j = 1 | {\bf v}, s)= sigmoid ( c_j + (\frac{{\bf v}}{{\boldsymbol \sigma}^2})^{T} (W_s)_{:j})

sigmoid(x) = \frac{1}{1 + \exp(-x)}であり、,(W_s)_{:j}は行列{\bf W_s}j列目の成分を表します。


これを用いてh_j=1となる確率を計算し、それをもとにh_jの値をサンプリングします。

その期待値が\langle h_j \rangle_{data}です。



また、\langle \cdot \rangle_{model}はモデル中での期待値を表します。

例えば、\langle v_i \rangle_{model}はARBMモデル中でのv_iの期待値を意味します。

これを求めるには、\int_{-\infty}^{\infty} v p(v_i = v, {\bf h}, s) dvの値を計算する必要があります。

しかし、この値を計算するためにはp(v_i = v, {\bf h}, s)に出てくるZの計算が必要で、これにはあらゆる{\bf v}{\bf h}の組み合わせに対して\Phi({\bf v}, {\bf h}, s)の値を計算しなくてはいけません。

これだと計算に膨大な時間がかかるので、ここではContrastive Divergence法(CD法)という方法を用いて近似を行います。



CD法は\langle \cdot \rangle_{model}k回のサンプリングの期待値で近似します。

実用的にはk=1が用いられることが多いようです。

CD法の例として、k=1\langle v_i \rangle_{model}を近似する手順を示します。

  1. 入力された{\bf v}をもとにp(h_j = 1 | {\bf v}, s)を計算
  2. p(h_j = 1 | {\bf v}, s)をもとに{\bf h}をサンプリング
  3. サンプリングした {\bf h}を用いてp(v_i = v | {\bf h}, s)を計算
  4. p(v_i = v | {\bf h}, s)をもとにサンプリングしたv_i\langle v_i \rangle_{model}の近似値


p(v_i = v | {\bf h}, s)p(h_j = 1 | {\bf v}, s)と同様に、p({\bf v}, {\bf h}, s)から導出可能で、以下のような分布になります。

p(v_i = v |{\bf h}, s)=N(v|b_i + (W_s)_{i:}{\bf h}, \sigma_i^2)

ただし、N(\cdot|\mu, \sigma^2)は平均\mu、分散\sigma^2正規分布を表します。



これで各パラメータの微分が分かったので、ARBMを学習させることができます。

実装

ここではコードを示しながら実際にARBMをどのように訓練するか見ていきます。

ソースコード全体を見たい方はこちらをどうぞ。

モデルの初期化

class ARBM():
    def initialize(self, visible_size, hidden_size, num_s):
        self.visible_size = visible_size
        self.hidden_size = hidden_size
        self.num_s = num_s
        self.W = np.random.normal(0, 0.1, [visible_size, hidden_size])
        self.b = np.zeros([visible_size, 1], dtype = "float")
        self.c = np.zeros([hidden_size, 1], dtype = "float")
        self.z = np.zeros([visible_size, 1], dtype = "float")

        self.A = np.random.normal(0, 0.1, [visible_size, visible_size, num_s])
        self.B = np.random.normal(0, 0.1, [visible_size, hidden_size, num_s])
        for i in range(num_s):
            vec = np.random.rand(visible_size)
            self.A[:,:,i] = np.diag(vec)

num_sは参照話者の数Svisible_sizeは可視層のサイズnhidden_sizeは隠れ層のサイズmを表しています。

また、記号についてはWは先ほどの数式で出てきた\bar{{\bf W}}を表し、それ以外は数式で出てきたのと同じ文字を使っています。

行列A_sはコード中でA[:,:,s]と表しています。つまりA[i,j,s]とあればこれはA_s(i,j)成分を表します。

bcz0で初期化し、それ以外は平均0、分散0.1正規分布から取った値で初期化します。

また、A[:,:,s]は、過学習を防ぐために対角成分以外は0としました。

隠れ層のサンプリング

def my_exp(self, array):
    array = np.nan_to_num(array)
    array[np.where(array > 50)] = 50
    array[np.where(array < -50)] = -50
    return np.exp(array)
                                                                   
                                                                   
def get_hidden_prob(self, visible, s):
    weight = np.dot(self.A[:,:,s], self.W) + self.B[:,:,s]
    lam = self.c + np.dot(weight.T, visible / self.my_exp(self.z))
    prob = 1 / (1 + self.my_exp(-lam))
                                                                   
    return prob
                                                                   
                                                                   
def get_hidden_sample(self, prob):
    size = prob.shape[0]
    sample = np.zeros([size, 1])
                                                                   
    for i in range(size):
        if np.random.rand() <= prob[i]:
            sample[i] = 1
        else:
            sample[i] = 0
                                                                   
    return sample

ここではp(h_j = 1 | {\bf v}, s)の計算と、それをもとにした隠れ層の値のサンプリングを行っています。

コード中のmy_expは、e^xの計算結果が大きくなりすぎたり小さくなりすぎるのを防ぐための関数です。

get_hidden_probp(h_j = 1 | {\bf v}, s)の計算です。

weightが数式中の{\bf W}_sを表します。

計算の高速化の為にベクトル化をしているのでちょっと分かり辛いコードになっていますが、頑張って数式と照らし合わせながら見ると何をやってるか分かると思います。

get_hidden_sampleでは、get_hidden_probで計算した確率をもとに{\bf h}の値をサンプリングします。

可視層のサンプリング

def get_visible_sample(self, hidden, s):
    weight = np.dot(self.A[:,:,s], self.W) + self.B[:,:,s]
    mean = self.b + np.dot(weight, hidden)
    sample = np.random.normal(loc = mean, scale = self.my_exp(self.z) + np.spacing(1), size = [self.visible_size, 1])
                                                                                                                      
    return sample

p(v_i = v | {\bf h}, s)の計算とサンプリングを行うコードです。

ほぼ数式通りなので、特にこれ以上説明することは無いです。

学習

学習を行う部分のコードは長いので、分割しながら説明します。

def train(self, rate, v_input, s):
    batch_size = v_input.shape[2]
    visible_size = self.visible_size
    hidden_size = self.hidden_size
    update_W = np.zeros([visible_size, hidden_size])
    update_A = np.zeros([visible_size, visible_size])
    update_B = np.zeros([visible_size, hidden_size])
    update_b = np.zeros([visible_size, 1])
    update_c = np.zeros([hidden_size, 1])
    update_z = np.zeros([visible_size, 1])

v_inputは入力された可視層の値を表していて、[visible_size, 1, batch_size]の配列になっています。

update_*は各パラメータの更新量を表す変数で、最初に0に初期化しておきます。

次は、更新を行う為の値のサンプリングと、各パラメータの微分の値の計算です。

    weight = np.dot(self.A[:,:,s], self.W) + self.B[:,:,s]
    for i in range(batch_size):
        visible0 = v_input[:,:,i]
        prob0 = np.nan_to_num(self.get_hidden_prob(visible0, s))
        hidden0 = np.nan_to_num(self.get_hidden_sample(prob0))
        visible1 = np.nan_to_num(self.get_visible_sample(hidden0, s))
        prob1 = np.nan_to_num(self.get_hidden_prob(visible1, s))
                                                                                                                          
                                                                                                                          
        for j in range(self.num_s):
            update_W += np.dot(np.dot(self.A[:,:,s].T, visible0), hidden0.T) \
                - np.dot(np.dot(self.A[:,:,s].T, visible1), prob1.T)
                                                                                                                          
        update_A += np.dot(np.dot(weight, hidden0), visible0.T) \
            - np.dot(np.dot(weight, prob1), visible1.T)
                                                                                                                          
        update_B += np.dot(visible0, hidden0.T) \
            - np.dot(visible1, prob1.T)

                                                                                                                          
        sub = np.nan_to_num(visible0 - self.b)
        sub[np.where(sub > 1e+5)] = 1e+5
        exp_data = np.square(sub) / 2 - visible0 * np.dot(weight, hidden0)
                                                                                                                          
        sub = np.nan_to_num(visible1 - self.b)
        sub[np.where(sub > 1e+5)] = 1e+5
        exp_model = np.square(sub) / 2 - visible1 * np.dot(weight, prob1)

        update_z += exp_data - exp_model
        update_b += visible0 - visible1
        update_c += hidden0 - prob1
                                                                                                                          

visible0が入力された可視層の値、hidden0がそれをもとにサンプリングした隠れ層の値

visible1hidden0をもとにサンプリングした可視層の値、prob1visible1をもとにした隠れ層の値

という風になっています。

このうち、visible0hidden0を数式で出てきた\langle \cdot \rangle_{data}の値として使い、visible1prob1\langle \cdot \rangle_{model}の近似値として使います。

あとは数式通り各パラメータの微分の値を計算して、それをバッチサイズの分だけ繰り返します。



最後は各パラメータの更新です。

    exp = self.my_exp(self.z)
    self.W += rate * update_W / (batch_size * exp)
    self.A[:,:,s] += rate * np.diag(np.diag(update_A)) / (batch_size * exp)
    self.B[:,:,s] += rate * update_B / batch_size
                                                                                                                          
    self.b += rate * update_b / (batch_size * exp)
    self.c += rate * update_c / batch_size
    self.z += rate * update_z / (batch_size * exp)

先ほど計算した微分の値をbatch_sizeで割って、学習率rateをかけて、各パラメータの値を更新します。

A[:,:,s]だけは、過学習を防ぐために対角成分のみを更新します。

実際に実験してみる

このコードでちゃんと訓練できるか試す為に、適当に作ったデータで実験してみました。

if __name__ == "__main__":
    import matplotlib.pyplot as plt
    visible_size = 32
    hidden_size = 256
    num_s = 2
    model = ARBM()
    model.initialize(visible_size, hidden_size, num_s)

    #prepare data
    v1 = np.empty((visible_size, 0), "float")
    v2 = np.empty((visible_size, 0), "float")
    size = 50
    for i in range(size):
        v1 = np.hstack([v1, np.random.normal(0, 0.5, [visible_size, 1])])
        v2 = np.hstack([v2, np.random.normal(0, 1, [visible_size, 1])])

    
    #training
    error_list = []
    batch_size = 10
    repeat_num = 100
    rate = 0.0003
    for i in range(repeat_num):
        for j in range(int(np.ceil(size / float(batch_size)))):
            start = j * batch_size
            end = (j + 1) * batch_size
            model.train(rate, v1[:, start:end].reshape(visible_size, 1, batch_size), 0)
            model.train(rate, v2[:, start:end].reshape(visible_size, 1, batch_size), 1)

        #calc error
        err = 0.0
        for k in range(size):
            input_val = v1[:,k].reshape(visible_size, 1)
            err += np.sum(np.square(input_val - model.construct(input_val, 0, 0)))
            input_val = v2[:,k].reshape(visible_size, 1)
            err += np.sum(np.square(input_val - model.construct(input_val, 1, 1)))
        err /= (size * 2)
        error_list.append(err)
        print "i = " + str(i) + ": " + str(err)

    #plot error
    plt.plot(error_list)
    plt.xlabel("Training Step")
    plt.ylabel("Error")
    plt.show()

可視層のサイズを32、隠れ層のサイズを256とし、32次元のベクトルを50個ずつ用意してv1v2としました。

v1を訓練する際はs=0v2を訓練する際はs=1を用いています。

訓練がちゃんと行われているかの基準は、本当はp({\bf v}, {\bf h}, s)の尤度を求めるべきなのですが、それは計算が大変です。

なので、ここでは入力に使った値をモデルに代入して再構築された値と、入力との2乗誤差で判断しています。

値を再構築するために使うのがconstructメソッドで、中身はこうなっています。

def construct(self, visible, input_s, output_s):
    hidden_prob = self.get_hidden_prob(visible, input_s)
    return self.get_visible_value(hidden_prob, output_s)

実際に実行してみると、誤差はこんな感じになりました。

f:id:phoro3:20161121113957p:plain

ちゃんと誤差が小さくなって、訓練が上手くいってることが分かります。

まとめ

今回はRBMを拡張したARBMを紹介しました。

実際にARBMを使って声質変換をするのは次回の記事で書きます。

参考文献・サイト

話者適応型Restricted Boltzmann Machineを用いた声質変換の検討(中鹿ら, 2014)