初期値の違いによる Baum-Welch アルゴリズムの収束過程

Baum-Welch アルゴリズムは、初期値によって収束の過程が異なります。出力記号系列が複数の局所解を持つ場合には、初期値によっていずれか一つの解が求まります。また、収束に要する反復回数も初期値によって大きく異なることがあります。今回は、初期値を変えながら Baum-Welch アルゴリズムによるパラメータ推定を繰り返し、この様子をグラフに描画して眺めてみます。

パラメータ推定の対象として、1 と 2 を交互に繰り返す出力記号系列を考えます。系列長は 1,000 とします。MATLAB では次のように定義できます。

>> x = repmat([1 2], 1, 500);

パラメータ A, B, ρ を次のように設定すれば、この系列が得られる確率は P(x) = 1 になります。前向きアルゴリズムで α を計算すれば簡単に確認できます。これらの A, B, ρ が、この出力記号系列に対する最尤推定量になります。

>> A = [
  0 1
  1 0
];

>> B = [
  1 0
  0 1
];

>> rho = [1 0];

>> alpha = Forward(A, B, rho, x);

>> exp(alpha(:, end))
ans =
     0
     1

プログラムの実装変更

収束の様子を詳細に調べられるように、前回の記事で実装した Baum-Welch アルゴリズムのコードを次のように変更します。引数として、初期値 A1, B1, rho1 を与えられるようにするほか、収束条件の maxiter, epsilon も指定できるようにします。関数の戻り値は、推定結果の A, B, rho だけではなく、反復処理中の A, B, rho をすべて戻すように変更します。このほかの関数の実装は、前回の記事で実装したものと同じです。

function [A, B, rho, logLH] = BaumWelchTrace(x, A1, B1, rho1, maxiter, epsilon)
    A(:,:,1) = A1;
    B(:,:,1) = B1;
    rho(:,:,1) = rho1;
    a = Forward(A1, B1, rho1, x);
    b = Backward(A1, B1, rho1, x);
    logLH(1) = calcLikelihood(a);
    for i = 2:maxiter
        [A(:,:,i), B(:,:,i), rho(:,:,i)] = maximize(A(:,:,i-1), B(:,:,i-1), x, a, b);
        a = Forward(A(:,:,i), B(:,:,i), rho(:,:,i), x);
        b = Backward(A(:,:,i), B(:,:,i), rho(:,:,i), x);
        logLH(i) = calcLikelihood(a);
        if logLH(i) - logLH(i-1) < epsilon
            break;
        end
    end
end

遷移確率行列 A の初期値による収束過程の変化

はじめに、A の初期値によって収束過程が変化する様子をグラフに描いてみます。B と ρ については、すべての要素を 0.5 としたものを初期値にします。以下のようなグラフが描かれました。左のグラフは反復終了時の対数尤度です。a11 と a21 は行列 A の 1 行 1 列、2 行 1 列の各要素の値を表します。他の要素は a12 = 1 - a11, a22 = 1 - a21 で計算できるため、グラフには含めていません。a11 が 0 に近く a21 が 1 に近いほど最適解に近く、そのような部分では対数尤度が 0 になっています。これは所定の反復回数で最適解が得られたことを表します。右のグラフは収束までの反復回数です。最大反復回数を 100 回、収束条件を ε = 1e-10 として実行しており、右のグラフで 100 回で平らになっている部分は、所定の反復回数で収束しなかったことを表しています。

このグラフの 0.45 から 0.55 の範囲を拡大したものが次のとおりです。最大反復回数を 1,000 回に増やして計算しました。計算量が増えた分だけメッシュを少し粗くしています。

さらに 0.49 から 0.51 の範囲を拡大すると次のグラフになります。

次は、B の初期値を少しずらして B = [0.4 0.6; 0.55 0.45] としたグラフです。先ほどのグラフで中央部分に筋のように見えていたものが、B をずらすとこのように曲がってきました。

さらに B をずらして、B = [0.3 0.7; 0.6 0.4] にしたものです。

出力確率行列 B の初期値による収束過程の変化

次に、A の初期値を固定した状態で、B の初期値による収束過程の変化を見てみます。先ほどのグラフで a11 + a21 = 1 を満たすところで筋状に特徴的な様子が見られましたので、これに沿って A の初期値を設定してみました。

まず A = [0.2 0.8; 0.8 0.2] としたグラフです。ほとんどの B で最適解に収束していますが、b11 == b21 の場合には計算が進まないようです。

A = [0.5 0.5; 0.5 0.5] にすると次のようなグラフに変化します。b11 と b21 が近い値をとるときに、より収束に時間がかかるようになっています。

さらに A を動かして A = [0.6 0.4; 0.4 0.6] にすると、次のように 100 の反復では収束しないようになりました。

初期状態 ρ による収束過程の変化

ρ を変化させたときのグラフも作成してみました。以下は ρ = [0.6 0.4; 0.4 0.6] として計算したグラフです。B を変えたときと異なり、直線上のまま中心からずれている様子が見えます。

Baum-Welch アルゴリズムによる軌跡のプロット

最後に、特定の初期値から収束するまでの A, B, ρ の軌跡をプロットして、収束の過程を見てみます。以下のように適当に初期値を決めて、プログラムを実行します。

>> A1 = [
  0.8 0.2
  0.3 0.7
];

>> B1 = [
  0.45 0.55
  0.50 0.50
];

>> rho1 = [0.5 0.5];

>> [A, B, rho, logLH] = BaumWelchTrace(x, A1, B1, rho1, 1e8, 1e-10);

得られた結果から (a11, a21) と (b11, b21) をプロットすると、以下のようになりました。反復回数は 72,516 回で、最終的には最適解が得られています。

各反復での ρ の値は次のグラフのとおりでした。グラフの値は初期状態が 1 になる確率を表します。

グラフからも読み取れるように、Baum-Welch アルゴリズムで得られた解は以下のようになります。B と ρ は、この記事の最初に計算した最適解と列が逆になっています。これは、状態の 1 と 2 がそれぞれ出力記号の 2 と 1 に対応付けられたためで、実質的には同一の解を表します。出力記号系列 x が得られる確率も P(x) = 1 と求まり、最適解になっています。

>> A(:,:,end)
ans =
    0.0000    1.0000
    1.0000    0.0000

>> B(:,:,end)
ans =
    0.0000    1.0000
    1.0000    0.0000

>> rho(:,:,end)
ans =
    0.0000    1.0000