DEIM2020 F8-3(day2 p54)

変分ベイズにおける最適解探索効率の検証
岡 威久馬†

若林啓††† 筑波大学情報学群知識情報・図書館学類 〒 305–8550 茨城県つくば市春日 1-2†† 筑波大学図書館情報メディア系 〒 305–8550 茨城県つくば市春日 1-2
E-mail: †s1611490@u.tsukuba.ac.jp, ††kwakaba@slis.tsukuba.ac.jp
あらまし

データの背後に隠れたパターンを発見し，解釈可能な知識を抽出することは重要な課題である．このような課題に対して，確率モデルによる統計的データ分析に基づく手法が有効と考えられており，クラスタリングをおこなう際に用いられる混合ガウスモデル（GMM）や系列データを扱う際に用いられる隠れマルコフモデル（HMM）のようなモデルが提案されている．これらのモデルを最適化する手法として変分ベイズ法が広く用いられているが，変分ベイズ法は決定論的なアルゴリズムであり，一度局所最適解に陥ると抜け出すことができない．本研究では，この問題を解決するための手法として摂動を加えるアプローチに着目し，期待値のモンテカルロ近似を用いた手法や，確率的勾配降下法を用いた手法について，最適解の探索効率を検証する．いくつかの条件の下で行なった実験の結果を示し，それぞれの手法の特性について得られた示唆を報告する．
キーワード 統計的学習，機械学習，クラスタリングアルゴリズム, 変分ベイズ, GMM, HMMなる場合もある. GMM や HMM のようなモデルでは, 多数のパ

1 はじめに近年, 確率モデルの学習アルゴリズムを用いたデータ分析が注目されている. ビッグデータの背後に隠れたパターンを発見
し, 解釈可能な知識を抽出することは重要な課題である. この課
題は, 対象とするデータに応じて, 類似したデータのグループを発見するクラスタリングや, 系列データからパターンを見つける系列データマイニングなどと呼ばれている. このような課題に対して, 確率モデルによる統計的データ分析に基づく手法が有効と考えられており, クラスタリングをおこなう際に用いられる混合ガウスモデル (GMM) や系列データを扱う際に用いられる隠れマルコフモデル (HMM) のようなモデルが提案されている.

GMM や HMM といったモデルは, モデルの関数の形を決めるパラメータが存在する. 確率モデルを用いたデータ分析では,モデルのパラメータを最適化することによって, 未知のデータに対して正確な予測をおこなうことが可能となる. パラメータを最適化する手法は数多く存在し, その中でも反復的数値最適化手法である変分ベイズは, パラメータ最適化手法の中でも広く知られている手法である. そして, 変分ベイズは GMM や HMMのようなモデルに対して用いることができる [1] [2].

GMM や HMM のようなモデルがデータ分析において有効であるということは先行研究によって明らかになっている [3] [4]が, これらの手法は局所最適解に陥りやすいという問題を抱えている. 局所最適解とは, 最適化したい関数のすべての変域のうち, 一部の変域を取り出して考えたとき, その範囲の中で最適解となる場所である. この方法によって求めた最適解は, 一部の変域のみを調べているので, そこで最適化されたパラメータが全体の最適解とは限らない. 局所最適解において定まったパラ
メータが良い精度をもたらす場合もあるが, 精度の悪い結果にラメータを用いるので, 悪い局所最適解に陥りやすいという問
題がある.
局所最適解を抜け出し, それよりも良い最適解を探索する手法はいくつか提案されている. そのうちの一つに, 乱数による摂動を加えるというアプローチがあり, 確率的勾配降下法 (SGD)
や期待値のモンテカルロ近似などの手法が先行研究によって提案されている. しかし, これらの手法は局所最適解を抜け出し,より良い最適解を見つける手法として有効であるということは明らかになっているが, どれだけ効率的に最適解を探索できているかということは明らかになっていない. 探索効率を検証することは, データの学習を高速化することや, データに合わせた最適な学習アルゴリズムの選択することに対して有益である.
また, 機械学習アルゴリズムを用いてデータ分析をする際, 多層
ニューラルネットワークなどの手法において, SGD を合わせた手法がさまざまなデータに対して用いられるが, 実際にはそれが最適な選択肢であるのかはまだあまり明らかにされていない.
データサイズやデータの種類によっては, SGD よりも効率よく最適化がおこなえる手法が存在する可能性がある. よって, 探索
効率を明らかにし, それぞれの手法を比較することでデータに合わせたアルゴリズムを提案することができる. また, それぞれの手法の良いところを取り入れ, 新たなアルゴリズムを提案することも可能である.
本研究では, 変分ベイズ手法を用いて, 摂動を加えない学習ア
ルゴリズム, SGD, 期待値のモンテカルロ近似の 3 手法を GMMと HMM の 2 つのモデルに対して適用し, 最適解の探索効率を検証し明らかにする.

2 関連研究近年, 深層学習や確率モデルを用いた学習アルゴリズムに関する研究は盛んに研究行われており, 多くの手法が提案されている. 多層ニューラルネットワークを利用した学習アルゴリズ
ムや, 潜在変数を持った確率モデルに対して用いられる EM ア











ルゴリズムのような反復的数値最適化手法などがその例である.
しかし, これらの手法は, データの学習中に局所最適解に陥り学習が進まなくなるという問題を抱えている. この問題を解決するための手法はいくつか提案されており, 現在でも盛んに研究Nされている対象である.
例えば, 確率的勾配降下法 (SGD) は局所最適解を抜け出すための手法として広く用いられている. SGD は勾配法にランダム

図1

GMM のグラフィカルモデル

性を取り入れた手法のことで, 勾配法はある関数において最小値を求めるための最適化手法である. ある関数 f (x) を最小化しに置き換わり, それ以外の要素は 0 に置き換えられる. このため,たいとき, 初期値を x, f (x) の微分を ∂f, ステップサイズを ρ と∂x

hardEM のように 1 に割り当てられる要素は E ステップを実行するとき, 勾配法の更新式は以下のようになる.

した時点ではわからない. このランダム性があることにより, ノx′ = x − ρ

∂f
∂x

(1)

イズを加えることができ, 局所最適解を抜け出すことができる.

SAEM は, SEM を改良したアルゴリズムで, 負担率の要素を 0と表すことができる. ρ はハイパーパラメータであり, 任意の値と 1 に再割り当てした後, 確率的勾配降下法のようにステップを用いることができる. 式 (2.1) を反復的に用いることで最小サイズを導入し, 以下のように更新する.

値へとたどり着くことが可能となる. これを応用させた手法がSGD である. SGD はランダムに抽出した学習データのみを用いて求めた勾配に基づいてパラメータを更新する．このことから，ノイズを加える効果があり, より良い最適解へたどり着くことが可能となる. Smith [5] らは, 多層ニューラルネットワークにおいて, SGD を用いることで局所解よりも大域解に近づくことができることを示している. 本研究では, SGD の考え方を変分ベイズに適用し, ランダムにサンプリングしたデータのみで求めた勾配を用いる Stochastic Variational Inference(SVI) を最適解探
索の検証の対象とする.
また, Valentin ら [6] は, EM アルゴリズムにおいて, softEM とhardEM という 2 つの手法を比較している. softEM は, EM アル
ゴリズムと同様の手法である. 一方, hardEM は期待値を計算した時, 一番高い確率を持つ要素に 1 を割り当て, それ以外の要素には 0 を割り当てるという方法をとり, 負担率を再割り当てした後に M ステップを実行する. [6] の実験によると, hardEM の方が精度が良いことが示されている. この事実は, 本研究において重要なことである. 本研究では期待値のモンテカルロ法を用いた近似を利用するが, hardEM と同様に負担率を 2 値変数に近似する. 再割り当てを実行する時にランダムな摂動を加えるという点のみが hardEM と異なる点である. このため, 変分ベイズに SEM を導入した手法も hardEM と同様に精度が向上することが見込める.

ii
θi = (1 − ρi )θEM
+ ρi θSEM(2)

i
は EM アルゴリズムの M ス
ρi はステップサイズであり, θEMiは, M ステップ後に負担率
テップによって更新された値, θSEMの再割り当てが行われた値である. 最後に MCEM は, E ステッ
プをモンテカルロ法によって近似する手法である. E ステップの計算が計算量的に困難であるときに用いられる. この手法は,複数のサンプルをモンテカルロ法によって抽出することによって E ステップを実行する. 複数のサンプル抽出をおこなわず, 1
つのみ抽出する場合, MCEM は SEM へ帰着される.
この研究では, SEM が他の手法よりも良い推定をおこなっていることが報告された. [8] でも述べられているように, 通常のEM アルゴリズムは必ず局所最適解に陥るという問題があるが, 確率的に近似することで局所解から抜け出すことができる.

SEM が他の手法よりも良い結果を得られた理由は, この EM ア
ルゴリズムの問題をランダム性を取り入れることによって局所解を抜け出すことができたからである.
本研究では, 上で述べた SVI, SEM を変分ベイズに適用し, 最適解の探索効率を検証する. 変分ベイズとは, 未知パラメータに分布を仮定することで EM アルゴリズムを拡張する手法である.
複雑な確率モデルのパラメータを厳密な計算によって推定することが困難なとき, 変分ベイズを用いて分布を近似することにより, 現実的な時間でパラメータ推定をおこなうことができる.

David ら [7] は, 混合ガウスモデルにおいて, EM アルゴリズ
ム, EM アルゴリズムに期待値の近似を取り入れた手法 (SEM),

3 混合ガウスモデルにおける最適解探索効率の検証

SAEM, MCEM の 4 手法の学習を比較している. SEM は, E ス
テップを行なった後, 1 か 0 の再割り当てをおこなう S ステッ

3. 1 混合ガウスモデル (GMM)

プとして実行される. hardEM の場合では, E ステップで負担率

混合ガウスモデル (GMM) とは, 複数の正規分布が重なり合っを計算した後, それぞれのデータの負担率の 1 番高い確率を 1,たモデルのことで, とても幅広い分野で応用されているモデル

それ以外を 0 に再割り当てする. しかし, SEM では負担率の 1である. GMM のグラフィカルモデルを図 3.1 に示した. GMM

番大きい要素を必ず 1 に置き換えるのではなく, 負担率にしたのパラメータは, 正規分布の平均パラメータ µ, 精度パラメータ

がってランダムに要素を選択する. そして選択された要素が 1

Λ, そして観測することのできない潜在変数 Z が用いられる. 観測データベクトルを X = {x1 , x2 , · · · , xn }, 対応する潜在変数

3. 1. 2 近似分布 q(Z) の計算

を Z = {z1 , z2 , · · · , zn } として, これらのパラメータが与えられ

式 (3.9) のような分解仮定をすることで計算可能な解を得るた時, 観測データベクトルの条件付き分布は以下のように表さことができる. 各因子に対応する逐次更新式は, カルバックライ

れる. ただし, Λ は分散の逆数である.

ブラーダイバージェンスの最小化によって得られる一般的な形

p(x|Z, µ, Λ) =

N ∏
K
∏

を利用することで得ることができる. 具体的には, 最適解を得た
znk
N (xn |µk , Λ−1
k )

(3)

n=1 k=1

い因子以外の因子で期待値をとることで最適解を得ることができる. まず, 因子 q(Z) についての更新式を導出する. 最適な因

ここで,子分布の対数は以下のようになる.

N (xn |µk , Λ−1
k ) =

1
1
1
exp{− (x−µ)T Λ(x−µ)}2(2π)D/2 (|Λ−1 |)1/2(4)

ln q ∗ (Z) ∝ Eπ,µ,Λ [ln p(Z, π, µ, Λ)](11)= Eπ [ln p(Z|π)] + Eµ,Λ [ln p(X|Z, µ, Λ)] (12)である. さらに, 混合比 π が与えられた時の Z の条件付き分布
式 (3.10) から式 (3.11) の式変形は, Z に依存しない項を定数には次のように表すことができる.
N ∏K∏

p(Z|π) =

吸収させている. 式 (3.11) に式 (3.1), 式 (3.3) を代入して計算す
z
πknk

(5)

ると, 以下のようになる.

n=1 k=1

ln q ∗ (Z) ∝混合比 π はクラスタ数 K と同じ要素数のベクトルであり, デー

N ∑
K
∑

znk ln ωnk(13)

n=1 k=1タに対して, それぞれの混合分布がどの程度負担しているかの割合となっており, クラスタ数を K とすると以下の条件を満たす.
K∑

ただし, ln ωnk は次のような式である. また, 式中の D はデータの次元である.

πk = 1(6)E[ln πk ] + 12 E[ln |Λk |] − Dln(2π)
2(14)− 12 Eµk ,Λk [(xn − µk )T Λk (xn − µk )](15)

ln ωnk =k=1

さらに潜在変数 Z は, K 次元の二値確率変数である.
次に，各パラメータに対して事前分布を導入する．内容は [9]式 (3.12) の両辺の指数をとって正規化すると,を参照している.

q ∗ (Z) =平均パラメータ µ, 精度パラメータ Λ, 混合比 π それぞれに事前分布を仮定する. π の事前分布には α0 = {α1 , α2 , · · · , αK }のディリクレ分布を仮定する.

p(π) = Dir(π|α0 )ト分布を導入する.

ωnk
rnk = ∑K
j=1 ωnj

式 (3.15) より, rnk は負担率を計算していることがわかる. これは, EM アルゴリズムでいうところの, E ステップをおこなっている.

3. 1. 3 近似分布 q(π, µ, Λ) の分解に計算できる 3 つの統計量を定義する. これは, 後の計算式を簡

(8)

k=1

3. 1. 1 事後分布の分解仮定

単にするためである.

Nk =GMM で変分ベイズを適用するために全ての確率変数の同時分布を書き下す.

この同時分布は, 図 3.1 のグラフィカルモデルをもとに分解されることがわかる. 次に, X が観測されたときの潜在変数とパラメータに分解した変分近似を考える. 分解は以下のように仮定
される.

(10)

N
∑

rnk

(18)

n=1

x̄k =

N
1 ∑
rnk xnNk n=1

(19)

Sk =

N
1 ∑
rnk (xn − x̄k )(xn − x̄k )T
Nk n=1(20)p(Z, π, µ, Λ) = p(X|Z, µ, Λ)p(Z|π)p(π)p(µ|Λ)p(Λ) (9)q(Z, π, µ, Λ) = q(Z)q(π, µ, Λ)(17)q(π, µ, Λ) を求める前に, 前の節で求めた負担率 rnk をもと

p(µ, Λ) = p(µ|Λ)p(Λ)
N (µk |m0 , (β0 Λk )−1 )W (Λk |W0 , ν0 )(16)

n=1 k=1(7)

度パラメータ ν0 ，D × D の尺度行列 W0 のガウス-ウィシャーK∏

rnk znkとなる. ただし,同様に, パラメータ µ, Λ をもつ正規分布については, 分布の自由

=

N ∏
K
∏

次に, 近似分布 q(π, µ, Λ) について考える. これは q(Z) と同様で, 考えている因子以外の因子で期待値をとることで最適化することができる. 最適化な分布 q ∗ (π, µ, Λ) の対数は以下のようになる.

ln q ∗ (π, µ, Λ) ∝ ln p(π) +K∑

Wk−1 = W0−1 + Nk Sk +ln p(µk , Λk )k=1

+ EZ [ln p(Z|π)]
+

K ∑
N
∑

β0 Nk
(x̄k − m0 )(x̄k − m0 )T
β 0 + Nk(31)

νk = ν0 + Nk

E[znk ] ln N (xn |µk , Λ−1
k )

(21)

(32)

これらの更新式は, EM アルゴリズムでいうところの, M ステップの計算をしている. 最後に, 式 (3.26) 中の E[znk ] の具体的な

k=1 n=1

この式の右辺について, π と µ, Λ が別々の項に分解されている計算であるが, これは式 (3.13) の ln ωnk を正規化することで得ことがわかる. また, µ, Λ が k 個の和で表されていることより,られる.

指数をとると k 個の積になることがわかる. 以上の議論により,の Algorithm1 のようになる.

q(π, µ, Λ) は以下のように分解される.
q(π, µ, Λ) = q(π)

K∏

以上より, GMM における変分ベイズのアルゴリズムは以下

q(µk , Λk )

(22)

k=1

Algorithm 1 GMM における変分ベイズアルゴリズム
Require: Data {xn }N
n=1 , iteration number T
Initialize distributions q(π), q(Z), q(µ, Λ)

3. 1. 4 近似分布 q(π) の計算

for i = 1 to T do

式 (3.20) の q(π) について, 最適な分布を求める. 式 (3.19) から π に関係する項を取り出すと ln q ∗ (π) は,

ln q ∗ (π) ∝ ln p(π) + EZ [ln p(Z|π)]Update q(Z) by Eq.(16)
Update q(π), q(µ, Λ) by Eq.(25), (28)
end for(23)

右辺第 1 項は式 (3.5), 第 2 項は式 (3.3) を代入して期待値を計算すると,

3. 2 期待値のモンテカルロ近似の導入 (GMM)
期待値のモンテカルロ近似の手法を変分ベイズに導入する.

ln q ∗ (π) ∝ (α0 − 1)K∑

ln πk +K ∑N∑

rnk ln πk(24)k=1 n=1k=1q (π) = Dir(π|α)

手続きをおこなって近似分布を最適化するので, E ステップにおいて負担率を確率的 EM アルゴリズムと同様の考え方でハー両辺の指数をとると, ディリクレ分布になる.
∗

上でも述べたように, 変分ベイズでも EM アルゴリズムと同じドに再割り当てをすることは可能である. 実際に, David ら [10](25)は, トピックモデルにおいてモンテカルロ近似を用いている. 変分ベイズに負担率の近似を導入した手法のアルゴリズムは以下

ここで, α の k 番目の要素 αk は αk = α0 + Nk である.

の Algorithm2 のようになる.

3. 1. 5 近似分布 q(µ, Λ) の計算
q(π) の計算と同様, 式 (3.20) から µ, Λ に関係する項を取り出して計算する. そうすると, ln q ∗ (µk , Λk ) は,

ln q ∗ (µk , Λk ) ∝K∑

ln p(µk , Λk )+

K ∑
N
∑

手法
Require: Data {xn }N
n=1 , iteration number T

Initialize distributions q(π), q(Z), q(µ, Λ)E[znk ] ln N (xn |µk , Λ−1k )k=1 n=1

k=1

for i = 1 to T do

(26)
ここで, 式 (3.25) 第 1 項について, 式 (3.6) と同じ分解を考えて,p(µk , Λk ) = p(µk |Λk )p(Λk ) より,
ln q ∗ (µk , Λk ) ∝K∑for all n do
Choose random element rnk
rnj ← δj=k for all j = 1 to K
end for{ln p(µk |Λk ) + ln p(Λk )}Update q(Z)Update q(π), q(µ, Λ)

k=1+

Algorithm 2 変分ベイズに期待値のモンテカルロ法を導入した

K ∑
N
∑

E[znk ] ln N (xn |µk , Λ−1k )end for(27)

k=1 n=1

式 (3.26) に式 (3.6) を代入して両辺の指数をとり整理すると, この近似分布はガウス-ウィシャート分布となる.

3. 3 確率的勾配降下法 (SGD) の導入 (GMM)
変分ベイズにおいて SGD を導入する. Hoffman ら [11] は,q ∗ (µk , Λk ) = N (µk |mk , (βk Λk )−1 )W (Λk |Wk , νk ) (28)
各パラメータの更新式は以下のように定義される.

LDA というモデルにおいて, 変分ベイズに確率的勾配降下法を導入した手法 (SVI) を提案した. また, Foulds ら [12] は, LDAにおいて, SVI を発展させた手法を提案している. このように,

βk = β0 + Nk

(29)

1
mk =
(β0 m0 + Nk x̄k )
βk

(30)

SVI は活発に研究され有用な手法であることがわかる. また,
Chris [13] は, GMM において SVI を用いた計算手法を示している. 以下の式の展開は, [13] をもとにしている.

SVI では, まずデータをランダムに抽出する. 抽出するデータ

サイズは 1 つ, もしくは複数である. 複数のデータを抽出する場
合, 抽出したデータサイズをバッチサイズという. そして, E ス

表 1 各分散における各手法の平均対数尤度 (GMM)
分散

データ数

テップと M ステップは通常通りおこない, 最後にステップサイ
ズ ρ を用いて M ステップで更新したパラメータを更新する. tρ(t) = (t + τ )−κ-3.9024-4.6508-3.9304-4.6662

300 件 -4.9067-4.8740-5.04163,000 件 -4.8618-4.8618-4.8762

30,000 件 -4.8608-4.8608-4.8650

300 件 -5.0861-5.0840-5.1969

3,000 件 -5.0891-5.0895-5.2184

30,000 件 -5.0927-5.0920-5.2241

τ と κ はハイパーパラメータである.
更新する近似分布のパラメータをまとめて λ, t 回目のイテレー

5.0

ションにおいて M ステップで更新されたパラメータを λ̂ とす-4.4453

3,000 件 -3.9274

3.0

(33)

-4.0246

30,000 件 -3.9064

1.0

回目のイテレーションにおいて, ρ は次のように更新する.

VI VI+SEM VI+SGD

300 件 -4.0256

ると更新式は以下のようになる.
(t)λ
3. 4 実= (1 − ρt )λ(t−1)+ ρt λ̂

π

Π

Z1

Z2

X1

X2

・・・

(34)

・・・

Zn

験
Xn

3.2 節から 3.4 節で述べたそれぞれの手法を用いて, 最適解探索に関する実験をおこなった.

3. 5 実 験 方 法

µ1

Λ1

µ2

µn

図2

HMM のグラフィカルモデル

Λ2Λnそれぞれのモデルのデータへの当てはまりの良さを示す指標
として対数尤度を用いる. 対数尤度は対象としているモデルがどれだけデータに適合できているかを示す指標であり, データ

それらの値を平均した値となっている.
各値について見てみると, 分散が 1.0 かつデータが 3,000 件やがそれぞれ独立に生成されると仮定すると次のように表すこと

分散 3.0 かつデータが 3,00 件の値では SEM+VI の手法が最もができる.

良い値となっている．一方で，VI ＋ SGD の対数尤度はどの分

ln p(X|θ) =N∑

ln p(xi |θ)

(35)

n=1

散の値，データ件数を比較しても最も低い値となった.
次に，各手法において正解パラメータと推定パラメータの差

X = {x1 , x2 , · · · , xn } であり, θ はパラメータベクトルである.

を計算した値を下の表に示す．ここでは，データ数 30,000 件，

GMM におけるパラメータベクトルは, θ = {π, µ, λ} である.

分散 5.0 の場合についての結果となっている．

イテレーション回数は 1,000 回で固定とする. 収束判定を用いる場合, 各手法で収束する速さが異なるため比較する際に不表 2 正解パラメータと推定パラメータのノルムの差 (GMM)µ

都合である. イテレーション回数を固定すれば, 各手法を同じ条VI

Λ0.3506

1.222

件で比較することができる. また，対数尤度の他に，パラメー

VI+SEM 0.0795

0.0191

タ推定をおこないその推定したパラメータ行列のフロベニウス

VI+SGD 2.0551 16.3804

距離を求め，正解パラメータのフロベニウス距離と比較する．

3. 5. 1 実験データ
実験に用いるデータは, 3 つの等方的な正規分布が重なった 2

表の値から，SEM を用いた手法の推定結果が最も正解パラ
メータとの距離が近いという結果となった．

次元混合分布からランダムに生成したデータを使う. 生成した分布の平均パラメータは, [0, 0], [4, 4], [7, 1] とし, 分散パラメー
タは, 対角成分が全て 1.0 の対角行列を用いる. また, 分散パラ

4 隠れマルコフモデル (HMM) における最適解探
索効率の検証メータは, 1.0 以外に 3.0, 5.0 と値を変えて実験をおこなった. また, データサイズは, 300 件, 3,000 件, 30,000 件の 3 種類である.
この実験は, 分散がそれぞれ 1.0, 3.0, 5.0 の場合で異なるデータを用いて 5 回ずつ実験をおこなった.

4. 1 隠れマルコフモデル (HMM)
天気データやセンサーデータなど, 系列データを扱う際に用いられるモデルが隠れマルコフモデル (HMM) である. HMM の3. 5. 2 実 験 結 果

グラフィカルモデルは図 4.1 のようになる. 観測データ xi はそ実験をおこなった結果, 各分散における平均対数尤度の値

れぞれに対応する潜在変数 zi , µi , Λi から生成される. GMM では表 1 のようになった. また, 表について, VI は変分ベイズ,

は潜在変数 Z を導入していたが, HMM でも同様に潜在変数を

VI+SEM は変分ベイズに期待値のモンテカルロ法を合わせた手導入する. しかし, HMM では, 各潜在変数 zn はそれぞれ前の状法, VI+SGD は変分ベイズに SGD を合わせた手法を表してい態の潜在変数に依存するという仮定をおく. そうすると, 遷移確る. 縦軸はデータの件数を表している. 表 1 に記載されている対率 Π が与えられたときの潜在変数 Z の条件付き確率は以下の数尤度の値は, 5 種類の異なるデータに対して対数尤度を求め,ようになる.

p(Z|Π) = p(z1 |π)N∏p(zn |zn−1 , Πn )ただし, bn,k , ak,j は以下である.

(36)

n=2

遷移確率は, Π2 , ..., ΠK を行にもつ K × K の行列である. π は

ak,j = exp(ln E[ln πk,j ]) = Ψ(αk,j ) − Ψ(

αk,j )(46)

j=1

初期確率である. 観測モデルは扱うデータの種類によって異な

bn,k = exp(ln E[p(xn |µk , Λk )])るが, 本研究では観測データは連続量とする. したがって, 観測(47)となる. 式 (4.13) の最右辺 Ψ はディガンマ関数を表している.

モデルは正規分布とする.

p(x|Z, µ, Λ) =

K
∑

N ∏
K
∏

詳細は [9] などを参照されたい.
znk
N (xn |µk , Λ−1
k )

(37)

4. 2. 4 近似分布 q(µ, Λ) の計算
q(µ, Λ) も q(Z) と同様に潜在変数 Z 以外の変数で期待値をn=1 k=1

計算することにより求めることができる. 計算結果は, GMM で

ここで, 正規分布は式 (3.2) の通りである.

導出したときと同様, ガウス-ウィシャート分布となる.

4. 2 HMM における変分ベイズ
4. 2. 1 事前分布の導入

q ∗ (µk , Λk ) = N (µk |mk , (βk Λk )−1 )W (Λk |Wk , νk ) (48)

HMM も GMM と同様に, パラメータに分布を仮定すること

4. 2. 5 近似分布 q(Π) の計算

で変分ベイズを用いることができる. 遷移確率 Π には, ディリ

変数 Π 以外で期待値を計算すると, ディリクレ分布となる.

クレ分布を仮定する. 式は,K∏

q ∗ (Π) = Dir(Π|α)
Dir(Πk |α0k )(38)

α0k = {α0k,1 , ..., α0k,K }(39)p(Π) =

k=1となる. また, µ と Λ の同時分布 p(µ, Λ) は, GMM と同様の事

4. 3 E ステップの計算
GMM の E ステップでは負担率を計算していた. HMM の E
ステップでも同様に負担率を計算する. 負担率 rn,k は以下のよ
うな式になる.

αn,k βn,krn,k = ∑K
k=1 αn,k βn,k

前分布なので,p(µ, Λ) = p(µ|Λ)p(Λ)
=

K
∏

(49)

(50)

式 (4.15) の αn,k と βn,k は, Baum-Welch アルゴリズムを用いる−1

N (µk |m0 , (β0 Λk ))W (Λk |W0 , ν0 ) (40)ことで表すことができる. 導出方法は, [15] や [16] を利用している. [15] では, 正規分布ではなく, ポアソン分布を用いた導出方

k=1

4. 2. 2 事後分布の分解仮定

法であるが, 負担率の計算方法は同じである. αn,k と βn,k は以

まず, 全ての確率変数の同時分布の分解は以下のようになる.

下のように表される.

p(X, Π, Z, µ, Λ) = p(X|Z, µ, Λ)p(Z|Π)p(Π)p(µ|Λ)p(Λ)αn,k = bn,j

(41)

K
∑

an,k αn−1,k(51)

βn+1,k an,k bn+1,k

(52)

k=1

次に, データ X が観測されたときの変分近似を考える. 分解は

βn,k =

以下のように仮定する.

K
∑
k=1

q(Z, Π, µ, Λ) = q(Z)q(Π)q(µ, Λ)= q(Z)K∏q(Πk )

k=1

K
∏

(42)

式 (4.15), (4.16), (4.17) より, 負担率は再帰計算によって求められる.

q(µk , Λk )

(43)

k=1

4. 3. 1 M ステップの計算
E ステップと同様, GMM と同じように計算する. まず, あと

4. 2. 3 近似分布 q(Z) の計算
GMM において潜在変数の近似分布を最適化したときと同様に, HMM で扱う潜在変数 Z の近似分布は, Z 以外の変数で期待値を計算することにより最適化することができる.

の計算を簡単にするため, 3 つの統計量を定義する. これは, 実
際には GMM と同じ式で表すことができる. よって, 式 (3.16) から式 (3.18) を HMM の M ステップでも用いる.
さらに, HMM における観測データには正規分布を仮定していln q ∗ (Z) ∝ EΠ,µ,Λ [ln p(X, Z, Π, µ, Λ)]るので, HMM で更新する正規分布のパラメータも GMM で示= EΠ,µ,Λ [ln p(X|Z, µ, Λ)p(Z|Π)p(Π)p(µ|Λ)p(Λ)] した式と同様の式になる. よって, ガウス-ウィシャート分布の更新式は GMM と全く同じ式になり, 式 (3.28) から式 (3.31) に= Eµ,Λ [ln p(X|Z, µ, Λ)] + EΠ [ln p(Z|Π)](44)
よって表される.
式 (4.11) の計算を進め, 両辺の指数をとると以下のようになる.
また, 遷移確率 Π のディリクレ分布 q(Π) の更新式は以下のよ
なお, 以下の式の表記は部分的に [14] を用いている.

q ∗ (Z) =N ∏K∏(bn,k )

n=1 k=1

N ∏
K ∏
K
∏

(ak,j )zn,k ,z(n+1),j

n=1 k=1 j=1

うになる.

(45)

(0)

αk,j = αk,j +

N
−1
∑
n=1

h(zn,k , z(n+1),j )(53)

ただし, h(zn,k , z(n+1),j ) は以下のように定義される.

表 3 各分散における各手法の平均対数尤度 (HMM)
分散α(n−1),k an,k bk,j βn,jh(zn,k , z(n+1),j ) = ∑K ∑K
l=1
m=1 α(n−1),l al,m bn,m βn,m
(54)

1.0

データ数

GMM と同様, HMM でも期待値を近似する. 式 (4.17) の計算によって求めた負担率を, 1 と 0 の 2 値に再割り当てする.

5.0

4. 5 確率的勾配降下法 (SGD) の導入 (HMM)
HMM に対して SVI を導入する. この手法は, Foti ら [17] によって提案されている. まず, バッチサイズ S を定め, 全データ
から S だけデータをサンプリングする. E ステップではサンプリングされたデータを使って負担率を計算し, M ステップではE ステップで求めた負担率から統計量を計算し, 近似分布のパラメータを更新する. また, M ステップにおける各パラメータの更新では, GMM と同様ステップサイズ ρ を用いる.

4. 6 実-3.7750-4.2346

3,000 件 -3.7712-3.7708-3.7781

30,000 件 -3.7665-3.7664-3.7774

300 件 -4.9067-4.8740-5.04163,000 件 -4.8618-4.8618-4.8762

30,000 件 -4.8608-4.8608-4.8650

4. 4 期待値のモンテカルロ近似の導入 (HMM)
3.0

VI VI+SEM VI+SGD

300 件 -3.9238

300 件 -5.3593-5.3363-5.52393,000 件 -5.3503-5.3503-5.3635

30,000 件 -5.3494-5.3496-5.3522

SEM+VI が最も良い値を示し，VI+SGD が最も低い値を示した.
次に，正解パラメータと推定パラメータの距離の差を求めた
値を次の表に示す．また，データ数，分散の値は GMM と同じである．
表 4 正解パラメータと推定パラメータのノルムの差 (HMM)

µ

験

VI

4.2 節から 4.3 節について述べた手法について, GMM と同様

Λ

0.0050 0.0558

VI+SEM 0.0004 0.0249の実験をおこなった.

VI+SGD 0.3099 0.1009

4. 6. 1 実 験 方 法
GMM 同様, データへの当てはまりの良さの指標として対数
尤度を計算する. 式は (3.37) のようになる. また, イテレーショ

表より，SEM を用いた手法によって推定されたパラメータが最も正解パラメータに近いことがわかる．

ン回数は 1,000 回で固定とする. さらに，正解パラメータと推定パラメータのフロベニウス距離を比較する．

5 結論

本実験で用いるパラメータは, 平均パラメータ, 分散パラメー
タ, 遷移確率である. 正解の平均パラメータは以下のような行列
である.



0.0


 0.0
µ=

 9.0
11.0

0.0



た. GMM においては, 部分的に SEM が対数尤度, パラメータ推
11.0 

10.0 

定の両方で最も良い性能を示した. これは, 通常の変分ベイズでは, 対数尤度がある値で収束してしまい局所最適解に陥ってしまうが, SEM では期待値をランダムに近似することで摂動が加わ−1.0

さらに, 正解の分散パラメータは以下のような対角行列である. また, 分散は 1.0 の要素が 3.0, 5.0 の場合の合計 3 種類で実験をおこなった.



0



0.1

 0
Λ=

 0

0

0

0.1

0

0

0.1


0 


0 

0

0

0

0.1

0.7


 0.3
Π=

 0.0
0.2

り, 局所最適解をうまく抜け出すことができたからである. 一方
で SGD は，SEM 同様通常の変分ベイズよりも良い性能を示すかと思われたが, 比較した手法の中で最も性能が低い結果となった. SGD を用いた手法は，各イテレーションで全データから 1つもしくは複数データをランダムにサンプル抽出するが, 実験
結果の平均対数尤度やパラメータ推定がうまくいっていないことから，SGD による摂動の加え方が大きすぎると言える. SEM
は毎回全データを用いて学習するが, SGD は毎回異なるデータ

最後に遷移確率 Π は次のような値を持つ行列である.



本研究では, 変分ベイズにおいて, 摂動を加える各手法の最適解探索効率を GMM と HMM の 2 つのモデルにおいて検証しを抽出して学習するので SEM に比べてうまく学習が進まないのである. 同様の結果が HMM でも示された. HMM では，大部

SGD を用いた手法が最も性能の低い結果となった. GMM の場

0.5


0.0 


0.2 

0.2

0.6

0.2

0.0

0.5

0.2

0.3
0.0

0.1

データサイズは 300 件, 3,000 件, 30,000 件の 3 種類である.

分の実験設定において SEM を用いた手法が最も良い値を示し,合と同様の理由から SGD はうまく学習をおこなうことができなかったと言える.

2 つのモデルにおいて変分ベイズの各手法を検証した結果, 摂

4. 6. 2 実 験 結 果

動が比較的小さい SEM は性能がよく, 摂動を大きく加える SGD

実験の結果, 平均対数尤度は表 3 のようになった.

は加えない手法よりも性能が低いことがわかった. SGD は多層表を見てみると, 大部分の分散の値, データ件数においてニューラルネットワークなど様々な手法で応用されているとて

も広く知られた手法であり, その手法の有効性は数多くの文献によって示されている. しかし, 変分ベイズにおいて局所最適解を抜けだすための手法として考察してみると, 本研究で示した[13]

ように SGD よりも良い性能を示す手法が存在する. SGD は全データを使って学習することが計算量的に困難である際に用いると効果を発揮するが, 全データを使って学習をおこなうこと[14]ができる際には, SEM のような摂動をあまり加えない手法を用
[15]
[16]

いる方が良い.
今後の課題としては, 各手法の実行時間を計測し比較することである．SGD は大量のデータから一部のデータを取り出す操作をおこなうことにより学習を高速化する効果があるため，実行時間において比較すると良い結果となる可能性がある．また，
本研究では，GMM と HMM という比較的扱いやすいモデルを対象に実験をおこなったが，さらに複雑なモデルを対象として実験し，最適解探索効率に関する一般的な結論を得る必要がある．さらに, SEM を拡張した手法である MCEM や SEM を改良した SAEM を変分ベイズに導入し, 摂動の加わり方を調べる必
要もある.

謝

辞

本 研 究 の 一 部 は ，JSPS 科 研 費（ 課 題 番 号 JP16H02904，

JP19K20333）の助成によって行われた．
文

献[1] Zoubin Ghahramani and Matthew J. Beal. Variational inference forbayesian mixtures of factor analysers. In S. A. Solla, T. K. Leen,and K. Müller, editors, Advances in Neural Information Processing
Systems 12, pp. 449–455. MIT Press, 2000.
[2] ZOUBIN GHAHRAMANI. An introduction to hidden markov models and bayesian networks. International Journal of Pattern Recognition and Artificial Intelligence, Vol. 15, No. 01, pp. 9–42, 2001.
[3] Noam Shental, Aharon Bar-Hillel, Tomer Hertz, and Daphna Weinshall. Computing gaussian mixture models with em using equivalence constraints. In Proceedings of the 16th International Conference on Neural Information Processing Systems, NIPS’03, pp. 465–
472, 2003.
[4] Padhraic Smyth. Clustering sequences with hidden markov models.
In Proceedings of the 9th International Conference on Neural Information Processing Systems, NIPS’96, pp. 648–654, 1996.
[5] Sam Smith and Quoc V. Le. A bayesian perspective on generalizationand stochastic gradient descent. 2018.
[6] Valentin I. Spitkovsky, Hiyan Alshawi, Daniel Jurafsky, and Christopher D. Manning. Viterbi training improves unsupervised dependency parsing. In Proceedings of the Fourteenth Conference on Computational Natural Language Learning, pp. 9–17, Uppsala, Sweden,
July 2010. Association for Computational Linguistics.
[7] Didier Chauveau Gilles Celeux and Jean Diebolt. On stochastic versions of the em algorithm. Technical report, INRIA, 1995.
[8] 修功上田, 良平中野. 確定的アニーリング em アルゴリズム. 電
子情報通信学会論文誌. D-2, 情報・システム 2-情報処理, Vol. 80,
No. 1, pp. 267–276, jan 1997.
[9] C.M. Bishop. パターン認識と機械学習 (下). 丸善, 2012.
[10] David M. Mimno, Matthew D. Hoffman, and David M. Blei. Sparse
stochastic inference for latent dirichlet allocation. In Proc. ICML,
2012.
[11] Matthew D. Hoffman, David M. Blei, Chong Wang, and John Paisley. Stochastic variational inference. Journal of Machine Learning
Research, 2013.
[12] James Foulds, Levi Boyles, Christopher DuBois, Padhraic Smyth,and Max Welling. Stochastic collapsed variational bayesian infer-[17]

ence for latent dirichlet allocation. In Proceedings of the 19th ACM
SIGKDD International Conference on Knowledge Discovery andData Mining, KDD ’13, pp. 446–454, 2013.
In-Depth Variational Inference Tutorial. https://chrisdxie.
files.wordpress.com/2016/06/in-depth-variationalinference-tutorial.pdf.
C. A. McGrory and D. M. Titterington. Variational bayesian analysis
for hidden markov models. Australian & New Zealand Journal of
Statistics, Vol. 51, No. 2, pp. 227–244, 2009.
須山敦志. ベイズ推論による機械学習入門. 講談社, 2017.
Christian Gruhl and Bernhard Sick. Variational bayesian inference
for hidden markov models with multivariate gaussian output distributions. 05 2016.
Nick Foti, Jason Xu, Dillon Laird, and Emily Fox. Stochastic variational inference for hidden markov models. In Z. Ghahramani,
M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 27, pp.
3599–3607. Curran Associates, Inc., 2014.