2020-12-26

数理最適化をしっかり学ぶために2段階単体法を実装

連続最適化アルゴリズム

はじめに

数理最適化には興味があってこのブログでもいくつか記事を書いてきたが、つまみ食いばかりして理解があいまいな部分もあるので、もう少しちゃんと基本から学ぶため以下の書籍を読み進めている。

しっかり学ぶ数理最適化モデルからアルゴリズムまで (KS情報科学専門書)

作者:梅谷俊治
発売日: 2020/10/26
メディア: 単行本

ということで今回は、この書籍の2章に書かれている線形計画問題を解くアルゴリズム『2段階単体法』のまとめとPythonによる実装を記す。

線形計画問題

解きたい問題

目的関数が線形関数で、制約条件が全て線形の等式もしくは不等式である最適化問題を線形計画問題という。線形計画問題は全て以下の標準形として定式化できる。

$\displaystyle \begin{aligned} \max_{{\bf x}} &\ \ {\bf c}^T {\bf x} \\ {\rm s.t.} &\ \ {\bf A}{\bf x} \leq {\bf b} \\ &\ \ {\bf x} \geq {\bf 0} \end{aligned} \tag{1}$

最適化問題の変数が $n$ 個、非負制約を除く制約条件が $m$ 個であれば、 ${\bf x} \in \mathbb{R}^n, \ {\bf A} \in \mathbb{R}^{m \times n}, \ {\bf b} \in \mathbb{R}^m$ である。ここで、スラック変数 ${\bf s} \in \mathbb{R}^m$ を導入し不等式制約を等式制約に変換する。

$\displaystyle \begin{aligned} \max_{{\bf x}, \ {\bf s}} &\ \ {\bf c}^T {\bf x} \\ {\rm s.t.} &\ \ {\bf A}{\bf x} + {\bf s} = {\bf b} \\ &\ \ {\bf x} \geq {\bf 0} \\ &\ \ {\bf s} \geq {\bf 0} \end{aligned} \tag{2}$

この問題（2）の解は元の問題と同じ解が得られる。

最適解

線形計画問題（1）の解の存在性は、 ${\bf c}, \ {\bf A}, \ {\bf b}$ の値に応じて以下の3パターンに分かれる。

制約を満たす解が存在しない（実行不能）
目的関数を無限に大きくできる（非有界）
最適解がある

そして最適解があるパターンの場合、その最適解は不等式制約で形成される $n$ 次元凸多面体上の頂点上にある。従って最適化する上では凸多面体上の頂点のみを調べればよい。これを効率よく探索するのが単体法である。

また、凸多面体上の頂点ということは式（1）において、非負制約を含む $n + m$ 個の制約条件のうち少なくとも $n$ 個は等式を満たすことになる。さらに、これを式（2）で考えると等式を満たすということは対応する変数もしくはスラック変数が $0$ になることを意味する*1。つまり、最適解 ${\bf x}^{*}, \ {\bf s}^{*}$ の合わせて $n + m$ 個の要素のうち少なくとも $n$ 個は $0$ になる。

アルゴリズム

単体法の考え方

細かい解説は上記の書籍に譲るとして、ここでは大まかな考え方を述べる。まずは ${\bf x} = {\bf 0}$ が制約条件を満たすとして話を進める。この前提は ${\bf b}$ に負値が含まれていると成り立たないが、その場合の対処法は後述する。とりあえず ${\bf b}$ に負値は含まれないとする。

単体法では最適解候補である凸多面体上の1つの頂点から始め（上記の前提が満たされているなら、原点 ${\bf x} = {\bf 0}$ も頂点の1つなので原点を初期値としてよい）、目的関数が改善する隣接する頂点に移動するという操作を改善できなくなるまで繰り返す。先述の通り、頂点ということは $n$ 個の制約条件が等式を満たすことになるが、隣接する頂点に移動するということは、この等式を満たす制約を1つ入れ替える操作に対応する。

この操作を数式を使って記述する。まず式（2）を少し記号などを書き換える。

$\displaystyle \begin{aligned} \max_{{\bf x}_B, \ {\bf x}_N} &\ \ {\bf c}_N^T {\bf x}_N + {\bf c}_B^T {\bf x}_B \\ {\rm s.t.} &\ \ {\bf N}{\bf x}_N + {\bf B}{\bf x}_B = {\bf b} \\ &\ \ {\bf x}_B \geq {\bf 0} \\ &\ \ {\bf x}_N \geq {\bf 0} \end{aligned} \tag{3}$

${\bf x}_N = {\bf x}, \ {\bf x}_B = {\bf s}, \ {\bf c}_N = {\bf c}, \ {\bf c}_B = {\bf 0}, \ {\bf N} = {\bf A}, \ {\bf B} = {\bf I}$ とすれば式（2）になる。 ${\bf x}_N$ は等式を満たすとした制約に対応する変数を表す。先述の通り等式を満たす制約に対応する変数は $0$ になるので、 ${\bf x}_N$ は ${\bf 0}$ に固定される。この ${\bf x}_N$ に割り当てれた（つまり $0$ に固定する）変数を非基底変数、 ${\bf x}_B$ に割り当てられた変数を基底変数という。また、 ${\bf N}$ を非基底行列、 ${\bf B}$ を基底行列と呼ぶ。単体法は ${\bf x}_N$ を ${\bf 0}$ にするのが本当に最適なのかを調べ、最適でない（まだ大きくできる）なら非基底変数と基底変数を1つ入れ替える。それに伴って対応する ${\bf N}$ と ${\bf B}$ の列を交換し、対応する ${\bf c}_N$ と ${\bf c}_B$ の要素を交換する。これが凸多面体上の頂点を移動する操作に当たる。

例えば、初期値が原点から始まるとすれば、最初は非基底変数には元の最適化問題に含まれる変数 ${\bf x}$ が割り当てられており、基底変数にはスラック変数 ${\bf s}$ が割り当てられている。そして、この状態が最適でなければ基底変数、非基底変数の中から変数を1つずつ適切に選び交換する。また、それらに対応する基底行列、非基底行列の列および係数の要素を交換する。この時の操作のイメージはこんな感じ。

f:id:YamagenSakam:20201229094122p:plain

では、どうやって現在の状態が最適か否かを判断するか？入れ替える基底変数、非基底変数はどう選ぶか？というのが次のお話。

交換する変数の選択と停止、非有界の判定

まず ${\bf B}$ が正則だとすると式（3）における等式制約より、

$\displaystyle {\bf x}_B = {\bf B}^{-1} {\bf b} - {\bf B}^{-1} {\bf N} {\bf x}_N\tag{4}$

が成り立つ。これを式（3）の目的関数に代入すると、

$\displaystyle {\bf c}_B^T {\bf B}^{-1} {\bf b} + {\bar{\bf c}_N} ^T {\bf x}_N \tag{5}$

が得られる。ここで、

$\displaystyle {\bar{\bf c}_N} = {\bf c}_N - {\bf N}^T {\bf B}^{-T} {\bf c}_B \tag{6}$

である。式（5）, （6）より次のことがわかる。

${\bar{\bf c}_N} \leq {\bf 0}$ なら ${\bf x}_N = {\bf 0}$ が最適でこれ以上大きくできない。つまり、この時にアルゴリズムは終了。
${\bar{\bf c}_N}$ に正の要素が含まれるなら、その要素に対応する ${\bf x}_N$ の要素を $0$ より大きな値にすると改善できる。従ってその変数は基底変数にする方がよい。

これで、アルゴリズムの終了条件と交換するべき非基底変数がわかった。では次に交換の対象となる基底変数はどう選ぶべきか？今、選んだ非基底変数を $x_k$ とすると、当然 $x_k$ は大きくした方がよい。一方で、式（4）より、 $x_k$ を大きくしすぎて ${\bf x}_B$ の要素が負になってはいけない（非負制約があるため）。 ${\bf B}^{-1} {\bf N}$ の $x_k$ に対応する列を ${\bar{\bf a}}_k$ とすると、

$\displaystyle {\bf B}^{-1} {\bf b} - x_k {\bar{\bf a}}_k \geq {\bf 0} \tag{7}$

を満たすように $x_k$ を定める必要がある。ここで、 ${\bar{\bf a}}_k \leq {\bf 0}$ ならば $x_k$ をいくらでも大きくできる。つまり、非有界ということになる。

一方、 ${\bar{\bf a}}_k$ の中に正の要素が $m_{+}$ 個含まれる場合、その正の要素の1つを ${\bar a}_{ik}^+$ 、それに対応する ${\bar{\bf b}} = {\bf B}^{-1} {\bf b}$ の要素を ${\bar b}_{i}^+$ とすると、

$\displaystyle 0 \leq x_k \leq \frac{\bar{b}_{i}^+}{\bar{a}_{ik}^+} \tag{8}$

が全ての $i = 1 \cdots m_{+}$ について成り立っている必要がある。実行可能解を移動してる限り式（4）と ${\bf x}_N = {\bf 0}$ より ${\bar{\bf b}} \geq {\bf 0}$ なので式（8）を満たす $x_k$ は存在し、 $x_k$ は以下の値にまで大きくできる。

$\displaystyle x_k = \min_i{\frac{\bar{b}_i^+}{\bar{a}_{ik}^+}} \tag{9}$

そして、この式（9）を最小にする $i$ に対応する基底変数は $0$ になる。つまり非基底変数になる。よって、交換対象の基底変数はこの $i$ に対応する変数を選ぶことになる。

交換する非基底変数、基底変数の選び方をまとめると、

非基底変数からは式（6）の正の要素に対応する変数を選ぶ。なければその時が最適であるためアルゴリズム停止。
基底変数からは式（9）を最小にする $i$ に対応する変数を選ぶ。 ${\bar{\bf a}}_k \leq {\bf 0}$ ならば非有界。

初期化（2段階単体法）

ここまでは、 ${\bf x} = {\bf 0}$ が制約を満たすとしこれを初期値として話を進めてきたが、 ${\bf b}$ に負値が含まれる場合 ${\bf x} = {\bf 0}$ は制約を満たさない。そもそも制約を満たす解などなく実行不能である可能性もある。そういことを考慮すると、まず実行可能解が存在するかを判断し、存在するなら適切な初期値を見つける必要がある。そのために以下の補助問題を考える。*2

$\displaystyle \begin{aligned} \min_{{\bf x}, \ x_a} &\ \ x_a \\ {\rm s.t.} &\ \ {\bf A}{\bf x} - x_a {\bf I} \leq {\bf b} \\ &\ \ {\bf x} \geq {\bf 0} \\ & \ \ x_a \geq 0 \end{aligned} \tag{10}$

追加した変数 $x_a$ は人工変数と呼ばれ、式（10）はこの $x_a$ を大きくすれば制約を満たせるので実行可能解は存在する。補助問題は最適化問題（1）が解を持つためには、どれくらい制約条件を修正するべきか？を求める問題とみなせる。従って最適解が $x_a^{*}=0$ であれば制約条件を修正する必要がない、つまり制約条件を満たす実行可能な解が存在する。そして、この時の ${\bf x}^{*}$ は制約条件を満たすのでこれを初期値として本来解きたい問題を解くようにすればよい。逆に最適解で $x_a^{*}>0$ であれば実行不能である。

実装としては、補助問題を先述の単体法の手続きで解いて、最適解 $x_a^{*}=0$ となった時の ${\bf B}, \ {\bf N}, \ {\bf c}_B, \ {\bf c}_N$ から $x_a$ に対応する列および要素を抜いたものを初期値としてもう一度単体法を行えばよい。このようにまず補助問題を解いて初期値を求め、その後解きたい問題を解く手続きを2段階単体法という。

ただし、 ${\bf b}$ には負値が含まれる場合、この補助問題についても原点 ${\bf x} = {\bf 0}, x_a = 0$ は制約を満たさない。そこで、 ${\bf b}$ の最小要素に対応するスラック変数を非基底変数、人工変数を基底変数とする入れ替えを行う。こうすれば実行可能になるのでこれを補助問題の初期値とすればよい。

処理の流れまとめ

ということで、2段階単体法の処理をまとめると以下のようになる。

2段階単体法の流れ

式（10）の補助問題を単体法で解く。
最適値が0でなければ実行不能として終了。
最適値が0になった場合、その時の非基底変数、基底変数の割り当てに基づき、 ${\bf B}, \ {\bf N}, \ {\bf c}_B, \ {\bf c}_N$ を並び替え。
これを初期値として単体法で式（1）の問題を解く。

そして単体法本体部分はこんな感じ。

単体法の流れ

${\bar{\bf b}} = {\bf B}^{-1} {\bf b}$ と式（6）で ${\bar {\bf c}}_N$ を求める。
${\bar {\bf c}}_N \leq {\bf 0}$ であれば最適解 ${\bf x}_B = {\bar{\bf b}}, {\bf x}_N = {\bf 0}$ 、最適値 ${\bf c}_B^{T} {\bf x}_B$ として終了。
そうでなければ、 ${\bar {\bf c}}_N$ の正の要素 $c_k$ に対応する非基底変数 $x_k$ を選択。
${\bf B}^{-1} {\bf N}$ の $x_k$ に対応する列 ${\bar{\bf a}}_k$ が ${\bar{\bf a}}_k \leq {\bf 0}$ であれば非有界として終了。
そうでなければ、式（9）の $i$ に対応する基底変数を選択。
選択した非基底、基底変数およびそれらに対応する ${\bf N}, {\bf B}$ の列、 ${\bf c}_N, {\bf c}_B$ の要素を交換。
1.に戻る。

実装

以下が実装した2段階単体法のpythonコード。

import numpy as np
import copy

def swapNB(cB, cN, B, N, Bidx, Nidx, swap_pos_B, swap_pos_N):
    # idx swap
    tmp = Nidx[swap_pos_N]
    Nidx[swap_pos_N] = Bidx[swap_pos_B]
    Bidx[swap_pos_B] = tmp
    
    # N, B swap
    tmp = np.copy(N[:, swap_pos_N])
    N[:, swap_pos_N] = B[:, swap_pos_B]
    B[:, swap_pos_B] = tmp

    #cN, cB, swap
    tmp = cN[swap_pos_N]
    cN[swap_pos_N] = cB[swap_pos_B]
    cB[swap_pos_B] = tmp   

# 単体法本体
def Simplex(cB_init, cN_init, B_init, N_init, Bidx_init, Nidx_init, b):
    Bidx = copy.copy(Bidx_init)
    Nidx = copy.copy(Nidx_init)
    cB = copy.copy(cB_init)
    cN = copy.copy(cN_init)
    B = copy.copy(B_init)
    N = copy.copy(N_init)
    while(1):
        invB = np.linalg.inv(B)
        b_ = invB @ b
        y = invB.T @ cB
        cN_ = cN - N.T @ y
        # 全て0以下なら停止
        if (cN_ <= 0).all() == True:
            xB = b_
            return cB @ xB, xB, cB, cN, B, N, Bidx, Nidx
        # cN_の正の要素の中で一番最初の要素に
        # 対応する変数交換する非基底変数として選択
        k = np.where(cN_ > 0)[0][0]
        
        ak_ = invB @ N[:, k]
        # 全て0以下なら非有界
        if((ak_ <= 0).all() == True):
            raise Exception('unbounded error')
        # ak_の中で、b_i/a_k_iが最小になる要素iに
        # 対応する変数を交換する基底変数として選択
        tmp = np.ones(b_.shape) * np.inf
        tmp[ak_ > 0] = b_[ak_ > 0] / ak_[ak_ > 0]
        l = np.argmin(tmp)
        swapNB(cB, cN, B, N, Bidx, Nidx, l, k)

# 初期値を探す（2段階単体法の1段目）
def SearchFeasibleDicstionary(cB_init, cN_init, B_init, N_init, 
                              Bidx_init, Nidx_init, b):
    # bに負が含まれる場合、実行可能な辞書を探す
    bmin = np.min(b)
    if bmin < 0:
        Bidx_aux = copy.copy(Bidx_init)
        Nidx = copy.copy(Nidx_init)
        cB_aux = copy.copy(cB_init)
        cN = copy.copy(cN_init)
        B_aux = copy.copy(B_init)
        N = copy.copy(N_init)
        
        # 人工変数のindex
        artificial_var_idx = cB_init.shape[0] + cN_init.shape[0]

        # 人工変数を非基底変数にする
        Nidx_aux = copy.copy(Nidx)
        Nidx_aux.append(artificial_var_idx)

        # 人工変数に対応する要素と列を追加
        cN_aux = np.zeros(cN.shape[0] + 1)
        # min x_a -> max -x_a　という最適化問題を解くため係数は-1
        cN_aux[cN.shape[0]] = -1
        N_aux = np.concatenate([N, np.ones((N.shape[0], 1)) * -1], axis = 1)

        # bの最小値に対応する基底変数と人工変数を入れ替え
        k = np.argmin(b)
        swapNB(cB_aux, cN_aux, B_aux, N_aux, 
                Bidx_aux, Nidx_aux, k, N_aux.shape[1] - 1)

        # 補助問題を解いて実行可能解を見つける
        z, _, res_cB, res_cN, res_B, res_N, res_Bidx, res_Nidx = \
         Simplex(cB_aux, cN_aux, B_aux, N_aux, Bidx_aux, Nidx_aux, b)
        
        # 最適解が負なら実行不能
        if z < 0:
            raise Exception('infeasible error')  
        # 得られた辞書のidxから人工変数を削除
        if artificial_var_idx in res_Nidx:
            # 人工変数が非基底変数にあるならそのまま削除
            res_Nidx.remove(artificial_var_idx)
        else:
            #退化して最適解の基底変数に人工変数が含まれる場合は
            # 非基底変数1つと交換する
            r = res_Bidx.index(artificial_var_idx)
            for i in range(len(res_Nidx)):
                swapNB(res_cB, res_cN, res_B, res_N, res_Bidx, res_Nidx, r, i)
                # 基底行列が正則か確かめる
                if(np.linalg.matrix_rank(res_B) == res_B.shape[0]):
                    res_Nidx.remove(artificial_var_idx)
                    break
                # 正則でなければ元に戻して、交換するindexを変えてもう一度
                swapNB(res_cB, res_cN, res_B, res_N, res_Bidx, res_Nidx, i, r)

        # 得られた基底、非基底変数の割り当てに基づいて
        # 基底、非基底行列、係数を再構成
        res_cNB = np.concatenate([cN_init, cB_init])[res_Nidx + res_Bidx]
        res_cN = res_cNB[:cN_init.shape[0]]
        res_cB = res_cNB[cN_init.shape[0]:]
        res_NB = np.concatenate([N_init, B_init],
                                        axis = 1)[:, res_Nidx + res_Bidx]
        res_N = res_NB[:, :cN_init.shape[0]]
        res_B = res_NB[:, cN_init.shape[0]:]
        return res_Bidx, res_Nidx, res_cB, res_cN, res_B, res_N
        
    # 負がなければそのままreturn
    return Bidx_init, Nidx_init, cB_init, cN_init, B_init, N_init

# max c.T @ x
# s.t. A @ x <= b
#      x >= 0
def SolveLP(c, A, b):
    cB = np.zeros(A.shape[0])
    cN = np.copy(c)
    N = np.copy(A)
    B = np.eye(A.shape[0])
    # 実際の変数のインデックス(0~n-1)
    actual_variable_idx = list(range(c.shape[0]))
    # スラック変数のインデックス(n~n+m-1)
    slack_variable_idx = list(range(c.shape[0], c.shape[0] + A.shape[0]))
    Nidx = copy.copy(actual_variable_idx)
    Bidx = copy.copy(slack_variable_idx)
   
    # 2段階単体法
    new_Bidx, new_Nidx, new_cB, new_cN, new_B, new_N = \
        SearchFeasibleDicstionary(cB, cN, B, N, Bidx, Nidx, b)
    z, xB_opt, cB_opt, cN_opt, B_opt, N_opt, Bidx_opt, Nidx_opt = \
        Simplex(new_cB, new_cN, new_B, new_N, new_Bidx, new_Nidx, b)
    
    # 得られた解から実際の変数（スラック変数以外の変数）を並び替えて出力
    x_opt = np.concatenate([xB_opt, np.zeros(cN_opt.shape[0])])
    x_actual_opt = x_opt[np.argsort(Bidx_opt + Nidx_opt)][actual_variable_idx]
    return z, x_actual_opt

実装は上述の処理を割と素直に書いた。現在どの変数が非基底変数でどの変数が基底変数にであるかについては、変数の添え字をNidx, Bidxというリストで持つことで管理している。 numpyは行優先にもかかわらず列を取り出す操作が多々あるので、効率的には行列を転置して行を取り出す実装にした方がよかったかもしれない*3。もう少し効率面で改善できそうだったり数値計算的に問題がありそうな箇所はあるが、まぁとりあえず書籍に載ってる演習問題を入力してみて、正しい答えを返すのでOKとする。

上述のソースコードに加えて検証用コード、検証に使ったテストケースを書いたテキストファイルを以下のgithub リポジトリに上げている。

github.com

1つこの実装方法でハマったのは、1段階目の補助問題を解く際に退化*4して人工変数が基底変数の中に含まれる最適値 $0$ の最適解が得られるケースがあったこと。補助問題の最適値が $0$ なのでもちろん実行可能なのだが、こうなると人工変数に対応する列を削除する際、基底行列から削除してしまい基底行列が正方でなくなってしまう。こうなった場合の対処として、人工変数といずれかの非基底変数を入れ替える処置を取っている。補助問題の定義より最適解が $0$ なら人工変数も $0$ になっているはずなので、同じく $0$ となる非基底変数を入れ替えても問題ない。ただし、入れ替えた結果の基底行列 ${\bf B}$ が正則でなくなってしまうこともある。そのため、1つずつ入れ替えを試して正則になったものを採用するようにした（この問題に対するもう少しスマートな解決方法はないものか）。

最後に

今回は単体法の実装を行ったが、それ以外にも内点法と呼ばれる線形計画問題を解くアルゴリズムもあり、こちらは非線形計画問題にも拡張されてるらしい。非線形計画問題の章を読んだらこちらの方も実装してみたいところ。

*1:非負制約の等式を満たす場合は当然それに対応する変数が $0$

*2:実装では目的関数に $-1$ をかけて $\max$ を求める問題に変換して解いている

*3:np.arrayのorderに'F'を指定してやればいいかも

*4:基底変数の値に $0$ が含まれる実行可能解が得られること

2020-04-16

再構成誤差における平均ベクトルとの差はどのように評価されるか

機械学習線形代数

はじめに

異常検知などのタスクでは、あらかじめ取得した学習データにより形成される主部分空間にテストデータを射影し、元のテストデータと射影後のデータの距離を異常度とすることがしばしば行われる。これを再構成誤差という（詳細は下記書籍を参照）。

入門機械学習による異常検知―Rによる実践ガイド

作者:剛, 井手
発売日: 2015/02/19
メディア: 単行本

もう少し具体的に書くと、学習データを ${\bf X} = \left\{{\bf x}_1, \cdots {\bf x}_N | {\bf x}_i \in {\bf R}^d, i = 1, \cdots, N \right\}$ 、テストデータを ${\bf x}_{*} \in {\bf R}^d$ とすると、学習データの平均ベクトルと共分散行列は

$\displaystyle {\bf m} = \frac{1}{N} \sum_{i = 1}^N {\bf x}_i \tag{1}$ $\displaystyle {\bf C} = \frac{1}{N} \sum_{i = 1}^N ({\bf x}_i - {\bf m}) ({\bf x}_i - {\bf m})^T \tag{2}$

として推定される。そして、この ${\bf C}$ の固有値分解 ${\bf C} = {\bf U} {\bf \Lambda} {\bf U}^T$ により、 $d$ 個の固有ベクトル ${\bf U} = \left[{\bf u}_1, \cdots , {\bf u}_d \right]$ を得るが、このうち固有値が大きい方から $p$ 個だけ対応する固有ベクトルを取ってくる。この $p$ 個の固有ベクトルが射影先の主部分空間の基底である。これを列ベクトルに並べた行列を ${\bf U}_p = \left[{\bf u}_1, \cdots , {\bf u}_p \right]$ とすると射影行列は ${\bf U}_{p} {\bf U}_{p}^{T}$ なので、テストデータ ${\bf x}_{*}$ の再構成誤差は、

$\displaystyle a({\bf x}_{*} ) = \parallel ({\bf x}_{*} - {\bf m}) - {\bf U}_{p} {\bf U}_{p}^{T} ({\bf x}_{*} - {\bf m}) \parallel_{2}^2 \tag{3}$

として計算される。さて、式（3）では平均ベクトル ${\bf m}$ で引いているが、このテストデータと平均ベクトルとの差 $({\bf x}_{*}-{\bf m})$ は再構成誤差にどのように効いてくるだろうか？直感的には ${\bf x}_{*}$ が ${\bf m}$ から遠ければ遠いほど再構成誤差は大きくなりそうだが*1実はそうではない、というのが今回のお話。結論を言ってしまうと、

テストデータ ${\bf x}_{*}$ が主部分空間の元である場合、平均ベクトルとは無関係に一定値になる
テストデータ ${\bf x}_{*}$ が主部分空間の元でない場合、平均ベクトルとの距離が大きくても再構成誤差は大きくならないことがある
主部分空間を平行移動させてできるアフィン部分空間内のテストデータは、全て同じ再構成誤差になる

である。以下、これを詳しく見ていく。

テストデータ ${\bf x}_{*}$ が主部分空間の元である場合

このケースはたまたま ${\bf x}_{*}$ が主部分空間の基底の線形結合で表せる、つまり ${\bf x}_{*} = {\bf U}_p {\bf a}$ となる ${\bf a} \in {\bf R}^{p}$ が存在するケースである。式（3）を少し変形すると、

$\displaystyle a({\bf x}_{*} ) =\parallel ({\bf x}_{*} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf x}_{*} ) - ({\bf m} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} ) \parallel_{2}^2 \tag{4}$

となる。式（4）の中で $({\bf x}_{*} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf x}_{*} )$ に着目し ${\bf U}_{p}^{T} {\bf U}_{p} = {\bf I}$ であることに注意すると、

$\displaystyle {\bf x}_{*} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf x}_{*} = {\bf x}_{*} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf U}_{p} {\bf a} = {\bf x}_{*} - {\bf U}_{p} {\bf a} = {\bf x}_{*} - {\bf x}_{*} = {\bf 0}$

なので、再構成誤差は ${\bf x}_{*}$ によらない一定値 $\parallel - ({\bf m} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} ) \parallel_{2}^2$ になる。従って、 ${\bf x}_{*}$ が主部分空間の元の場合 ${\bf m}$ からどれだけ離れていようが近くにあろうが全く関係ない。

このことを以下のpythonコードで確かめてみる。

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt

if __name__ == '__main__':
    # 学習データの生成とプロット
    np.random.seed(1234)
    N = 1000
    mu = [3, 10]
    Sig = [[8, 4],[4, 7]]
    X = np.random.multivariate_normal(mu, Sig, N).T
    plt.scatter(X[0, :], X[1, :], color = 'r', label="Train Data")

    # 学習データのサンプル平均推定とそのプロット
    m = np.average(X, axis = 1)
    plt.scatter(m[0], m[1], color = 'c', 
                s = 70, marker = "*", label = "Train Sample Mean")

    # 主成分分析で主部分空間推定とそのプロット
    Z = X - m[:, np.newaxis]
    C = Z @ Z.T/N
    lam, U = np.linalg.eig(C)
    Up = U[:, 0][:, np.newaxis]
    x = [-100 * U[0, 0], 0, 100 * U[0, 0]]
    y = [-100 * U[1, 0], 0, 100 * U[1, 0]]
    plt.plot(x, y, color = 'b', label="Principal Axes")
    
    # 学習データの平均から離れたテストデータの再構成誤差
    test_x1 = Up.squeeze() * -20
    plt.scatter(test_x1[0], test_x1[1], color = 'g', 
                s = 70, marker = "^", label = "Test Data 1")
    a1 = np.linalg.norm((test_x1 - m).T - Up @ Up.T @ (test_x1 - m)) ** 2
    print("Test Data 1 Reconstruct Error", a1)

    # 学習データの平均に近いテストデータの再構成誤差
    test_x2 = Up @ Up.T @ m 
    plt.scatter(test_x2[0], test_x2[1], color = 'g',
                s = 70, marker = "v", label = "Test Data 2")
    a2 = np.linalg.norm((test_x2 - m).T - Up @ Up.T @ (test_x2 - m)) ** 2
    print("Test Data 2 Reconstruct Error", a2)

    plt.xlim(-25, 25)
    plt.ylim(-25, 25)
    plt.legend()
    plt.show()

このコードにより、以下の図のような2次元のデータを生成し、主部分空間上にある2つのテストデータの再構成誤差を算出している。 f:id:YamagenSakam:20200415224630p:plain

そして、Test Data 1とTest Data 2の再構成誤差の計算結果は以下の通り。

Test Data 1 Reconstruct Error 28.670763684482974
Test Data 2 Reconstruct Error 28.670763684482964

Test Data 1はTest Data 2に比べて明らかに学習データの平均値から遠いにもかかわらず再構成誤差は同じ値になっている。

テストデータ ${\bf x}_{*}$ が主部分空間の元でない場合

一般の場合を考える。式（3）を展開すると以下を得る。

$\displaystyle \begin{aligned} a({\bf x}_{*} ) &= ({\bf x}_{*} - {\bf m})^T ({\bf I}- {\bf U}_{p} {\bf U}_{p}^{T})^T ({\bf I}- {\bf U}_{p} {\bf U}_{p}^{T}) ({\bf x}_{*} - {\bf m}) \\ &= ({\bf x}_{*} - {\bf m})^T ({\bf I}- {\bf U}_{p} {\bf U}_{p}^{T}) ({\bf x}_{*} - {\bf m}) \\ &= ({\bf x}_{*} - {\bf m})^T({\bf x}_{*} - {\bf m}) - ({\bf x}_{*} - {\bf m})^T{\bf U}_{p} {\bf U}_{p}^{T}({\bf x}_{*} - {\bf m}) \\ &= {\bf x}_{*}^T ({\bf I} -{\bf U}_{p} {\bf U}_{p}^T ){\bf x}_{*} - 2 {\bf m}^T ({\bf I} -{\bf U}_{p} {\bf U}_{p}^T ){\bf x}_{*} +{\bf m}^T ({\bf I} -{\bf U}_{p} {\bf U}_{p}^T ){\bf m} \end{aligned} \tag{5}$

第1項はテストデータの主部分空間からのはみだし具合だけで決まるので平均ベクトル ${\bf m}$ は無関係、第3項は定数、従って第2項に着目する。この第2項が大きければ再構成誤差としては小さく、逆に第2項が小さければ再構成誤差は大きくなる。 $({\bf I} -{\bf U}_{p} {\bf U}_{p}^T )$ は主部分空間の直交補空間への射影行列である。すなわち、第2項はテストデータ ${\bf x}_{*}$ を直交補空間へ射影したベクトルと平均ベクトルの内積を計算している。従って、どんなにテストデータが平均ベクトルから離れていても直交補空間へ射影されてから内積を計算するので射影された先が同じなら第2項は同じ値になる。また内積なので、射影後のベクトルと平均ベクトルの方向が重要で、方向が同じであるほど第2項は大きな値となる（再構成誤差は小さくなる）。

このことを確かめるために以下のようなpythonコードで実験。なお、主部分空間推定までは上のコードと共通なので省略。

    # 主部分空間推定までは上のコードと同じなので省略

    # 学習データの平均から離れているが学習データの平均と同じ方向に
    # 主部分空間から飛び出ているテストデータの再構成誤差
    test_x1 = (Up.squeeze() * -20) + 5 * U[:, 1]
    plt.scatter(test_x1[0], test_x1[1], color = 'g', 
                s = 70, marker = "^", label = "Test Data 1")
    a1 = np.linalg.norm((test_x1 - m).T - Up @ Up.T @ (test_x1 - m)) ** 2
    print("Test Data 1 Reconstruct Error", a1)

    # 学習データの平均から離れているが学習データの平均と逆方向に
    # 主部分空間から飛び出ているテストデータの再構成誤差
    test_x2 = (Up. squeeze() * -20) - 5 * U[:, 1]
    plt.scatter(test_x2[0], test_x2[1], color = 'g', 
                s = 70, marker = "v", label = "Test Data 2")
    a2 = np.linalg.norm((test_x2 - m).T - Up @ Up.T @ (test_x2 - m)) ** 2
    print("Test Data 2 Reconstruct Error", a2)

    # 学習データの平均と近く学習データの平均と同じ方向に
    # 主部分空間から飛び出ているテストデータの再構成誤差
    test_x3 = Up @ Up.T @ m + 5 * U[:, 1]
    plt.scatter(test_x3[0], test_x3[1], color = 'g', 
                s = 70, marker = ">", label = "Test Data 3")
    a3 = np.linalg.norm((test_x3 - m).T - Up @ Up.T @ (test_x3 - m)) ** 2
    print("Test Data 3 Reconstruct Error", a3)

    # 学習データの平均と近く学習データの平均と逆方向に
    # 主部分空間から飛び出ているテストデータの再構成誤差
    test_x4 = Up @ Up.T @ m - 5 * U[:, 1]
    plt.scatter(test_x4[0], test_x4[1], color = 'g', 
                s = 70, marker = "<", label = "Test Data 4")
    a4 = np.linalg.norm((test_x4 - m).T - Up @ Up.T @ (test_x4 - m)) ** 2
    print("Test Data 4 Reconstruct Error", a4)

    # 以降のプロットする部分も上のコードと同じなので省略

この実験では以下のように、主部分部分空間からそれと直交する方向に平行移動させた4つのテストデータで実験。 f:id:YamagenSakam:20200415230950p:plain

そして、この再構成誤差の計算結果は以下の通り。

Test Data 1 Reconstruct Error 0.12567643599399253
Test Data 2 Reconstruct Error 107.21585093297189
Test Data 3 Reconstruct Error 0.12567643599399056
Test Data 4 Reconstruct Error 107.21585093297192

このように、Test Data 1とTest Data 3は平均ベクトルと同じ方向に同じ量だけ主部分空間から平行移動させたため、再構成誤差としてかなり小さい同じ値になっている（Test Data 3はほとんど平均ベクトルと一致しているので再構成誤差は小さくなるのは当然だが、平均ベクトルから離れているTest Data 1もそれと全く同じ値になっている）。逆にTest Data 2とTest Data 4については、平均ベクトルと逆方向に主部分空間から平行移動させているため、再構成誤差は大きな値になっている（Test Data 4はTest Data 1よりも平均ベクトルとの距離という意味では小さいのに、再構成誤差はTest Data 4の方が大きい）。

ちなみに平均値 ${\bf m}$ が主部分空間の元である場合、直交補空間の元との内積は ${\bf 0}$ になるので、平均値とは無関係に元のデータの主部分空間からのはみ出具合のみで評価される。これは式（4）からも明らか。

再構成誤差の等高線

ここまでの結果を見ると主部分空間からどれだけ、どの方向に平行移動させたかが再構成誤差に効いてくるように思える。これを確かめるため再構成誤差の等高線を引いてみる。そのコードと実行結果は以下の通り。

    # 主部分空間推定までは上のコードと同じなので省略

    # 再構成誤差の等高線の計算と描画
    x1 = np.arange(-35, 35, 1) 
    x2 = np.arange(-35, 35, 1)    
    X1, X2 = np.meshgrid(x1, x2)
    A = np.zeros(X1.shape)
    for i in range(X1.shape[0]):
        for j in range(X1.shape[1]):
            test_x = np.array([X1[i, j], X2[i, j]])
            # 2乗があると等高線の間隔が広がりすぎてわかりにくいので省く
            A[i, j] = np.linalg.norm((test_x - m).T - Up @ Up.T @ (test_x - m))
    plt.contour(X1, X2, A, levels=50, cmap="viridis_r")
    plt.colorbar()

    # 以降のプロットする部分も上のコードと同じなので省略

f:id:YamagenSakam:20200415233632p:plain

これを見ると確かに主部分空間と平行な線上に同じ再構成誤差の値が乗っているように見える。実際、主部分空間を平行移動させた平面（アフィン部分空間）上のテストデータは同じ再構成誤差になる。以下これを示す。

任意の ${\bf x}_{*}$ は、主部分空間上の元 ${\bf U}_p {\bf a}$ とその直交補空間の元 ${\bf b}$ を使うと、

$\displaystyle {\bf x}_{*} = {\bf U}_p {\bf a} + {\bf b} \tag{6}$

と表せる。 ${\bf b}$ を加えるということは、主部分空間上の元を ${\bf b}$ だけ平行移動させて主部分空間から飛び出させることに等しい。ここで式（6）を式（3）に代入し ${\bf U}_{p}^T {\bf b} = {\bf 0}$ であることに注意して展開すると、

$\displaystyle \begin{aligned} a({\bf x}_{*} ) &= \parallel ( {\bf U}_p {\bf a} + {\bf b} - {\bf m}) - {\bf U}_{p} {\bf U}_{p}^{T} ( {\bf U}_p {\bf a} + {\bf b} - {\bf m}) \parallel_{2}^2 \\ &= \parallel ( {\bf U}_p {\bf a} + {\bf b} - {\bf m}) - ({\bf U}_{p} {\bf U}_{p}^{T} {\bf U}_p {\bf a} + {\bf U}_{p} {\bf U}_{p}^{T}{\bf b} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf m}) \parallel_{2}^2 \\ &= \parallel {\bf U}_p {\bf a} + {\bf b} - {\bf m} - {\bf U}_p {\bf a} + {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} \parallel_{2}^2 \\ &= \parallel {\bf b} - {\bf m}+ {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} \parallel_{2}^2 \\ &= {\bf b}^T {\bf b} - 2 {\bf b}^T{\bf m} + {\bf m}^T {\bf m} + 2 {\bf b}^T {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} - 2{\bf m}^T {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} + {\bf m}^T {\bf U}_{p}{\bf U}_{p}^{T} {\bf m}\\ &= \parallel {\bf b} - {\bf m} \parallel_{2}^2- {\bf m}^T {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} \end{aligned} \tag{7}$

となり、再構成誤差は平行移動成分 ${\bf b}$ と平均ベクトル ${\bf m}$ の距離のみに依存することがわかる。すなわち、平行移動成分 ${\bf b}$ が一定なら、主部分空間上の元をどれを取ってきても再構成誤差は同じである。つまりそれは、主部分空間を ${\bf b}$ だけ平行移動させたアフィン部分空間上の元の再構成誤差は、全て等しい値になることを意味する。

まとめ

ここで示したように平均ベクトルとテストデータの距離が大きくても再構成誤差は小さくなりうる。再構成誤差を使って異常検知などを行う場合、このことを頭の片隅にでも置いておかないと落とし穴にはまるかもしれない。

*1:というか、異常検知の枠組みだとそうなって欲しいケースが多いと思う

2020-03-31

2つの部分空間の共通部分の基底を求める

線形代数

はじめに

しばらくブログを更新してなかったので線形代数のネタを1つ。やりたいのは線形独立な $m$ 個のベクトル ${\bf u}_{1}, \cdots, {\bf u}_{m} \in R^{d}$ が張る部分空間と線形独立な $n$ 個のベクトル ${\bf v}_{1}, \cdots, {\bf v}_{n} \in R^{d}$ が張る部分空間の共通部分 ${\rm span} \left\{ {\bf u}_{1}, \cdots, {\bf u}_{m} \right\} \cap {\rm span} \left\{ {\bf v}_{1}, \cdots, {\bf v}_{n} \right\}$ の基底を求めること。

以降、行列 ${\bf U}$ を ${\bf U} = \left[ {\bf u}_{1}, \cdots, {\bf u}_{m} \right]$ 、行列 ${\bf V}$ を ${\bf V} = \left[ {\bf v}_{1}, \cdots, {\bf v}_{n} \right]$ とする。

方法1：零空間を使う方法

まず、 ${\rm span} \left\{ {\bf u}_{1}, \cdots, {\bf u}_{m} \right\}$ の任意の元は適当なベクトル ${\bf a} \in R^{m}$ を使って ${\bf U} {\bf a}$ と表せる。同様に ${\rm span} \left\{ {\bf v}_{1}, \cdots, {\bf v}_{n} \right\}$ の任意の元は適当なベクトル ${\bf b} \in R^{n}$ を使って ${\bf V} {\bf b}$ と表せる。共通部分の元は、

$\displaystyle {\bf U} {\bf a} = {\bf V} {\bf b}$

を満たす。これを変形すると、

$\displaystyle \begin{aligned} {\bf U} {\bf a} - {\bf V} {\bf b} &= &{\bf 0} \\ \begin{bmatrix} {\bf U} & {\bf V} \end{bmatrix} \begin{bmatrix} {\bf a} \\ -{\bf b} \end{bmatrix} & = & {\bf 0} \\ {\bf W} {\bf c} &=& {\bf 0} \end{aligned}$

最後の式は ${\bf W} = \left[{\bf U} \ \ {\bf V}\right]$ 、 ${\bf c} = \left[ {\bf a}^T \ \ -{\bf b}^T \right]^T$ とおいた。最後の式より、この条件を満たす ${\bf c}$ は ${\bf W}$ の零空間 $\ker ({\bf W})$ の元である。従って、この $\ker ({\bf W})$ の基底を求め、それぞれの基底について上から $m$ 次元分取り出し（これを ${\hat {\bf a}}_1, \cdots, {\hat {\bf a}}_k$ とする）、 ${\bf U}{\hat {\bf a}}_1,, \cdots, {\bf U}{\hat {\bf a}}_k$ とすれば共通部分の基底が求まる。

これをpythonで実装する場合は、scipy.linalg.null_spaceを使えば一発で零空間の基底を求められるのでこれを活用すればよい。ちなみに零空間の基底は、行列を特異値分解してその特異値ゼロに対応する右特異ベクトルとして得られる。実際、scipy.linalg.null_spaceも中では特異値分解を使ってるとのこと。これを試してみたコードは以下の通り（コメントにある通りこのコードの実験だと共通部分は3次元になるはず*1）。

import numpy as np
from scipy import linalg

def make_test_matrix(d, m, rand_seed = 1):
    np.random.seed(rand_seed)

    #Uはランダムに設定
    U = np.random.randn(d, m)
    #Vの2つのベクトルはランダム
    v1 = np.random.randn(d, 1)
    v2 = np.random.randn(d, 1)
    #他の3つは既存のUとVから構成(これで共通部分は3次元になるはず)
    v3 = ((7 * U[:, 0]) + (5 * U[:, 1]))[:, np.newaxis]
    v4 = ((1 * U[:, 1]) + (2 * U[:, 3]))[:, np.newaxis]
    v5 = ((3 * U[:, 3]) + (-2 * v1[:, 0]))[:, np.newaxis]
    V = np.concatenate([v1, v2, v3, v4, v5] , axis = 1)
    return U, V

if __name__ == '__main__':
    d = 15
    m = 5

    #テストデータの生成
    U, V = make_test_matrix(d, m, 1)

    #UとVを連結した行列の零空間の基底を計算
    W = np.concatenate([U, V], axis=1)
    C = linalg.null_space(W)
    #Uの係数のみ取り出して基底を計算
    A_hat = C[:m, :]
    basis = U @ A_hat

    #結果の表示
    print("A_hat")
    print(A_hat)
    print("basis of intersection")
    print(basis)

結果

A_hat
[[-7.16152592e-01  4.38912368e-01  1.34809890e-01]
 [-6.53812424e-01 -2.10948597e-01 -7.03781271e-02]
 [-1.63931368e-16 -2.84223600e-16 -1.73499453e-16]
 [-1.05903826e-01 -6.16391057e-01  6.45398134e-01]
 [-3.74700271e-16  5.09575021e-17 -3.66568755e-16]]
basis of intersection
[[-0.64967372  1.50336269 -0.43045986]
 [ 0.47368583 -1.57489224 -0.22715953]
 [ 0.34052603  1.31304872  0.09422749]
 [ 0.89595509 -0.47240257 -0.10889637]
 [-0.01343938 -1.03428625  0.09537105]
 [ 0.59837101 -0.10904928 -0.25641896]
 [ 0.84424754  0.30109334 -0.61081408]
 [ 0.56380177 -0.79295082  1.14816012]
 [ 0.53848895 -0.94016879  1.12891538]
 [ 0.31863787 -0.39392353 -0.02175922]
 [ 0.0523343   0.42138681 -0.16020853]
 [-0.99889429 -0.09553967  0.20435405]
 [-0.24730755 -0.41166361 -0.38226338]
 [-0.91721327 -1.61905413  1.32073652]
 [ 1.27124225 -1.06748573  0.40630033]]

期待通り3つの基底ベクトルが得られた。

方法2：正準角を使う方法

たぶん上で書いた方法がスタンダードなやり方だと思うが、次のような方法でも求められる。簡単のため各列ベクトルはそれぞれが張る部分空間の正規直交基底とする。つまり、 ${\bf U}^T {\bf U} = {\bf I}, \ \ {\bf V}^T {\bf V} = {\bf I}$ 。もし与えられているベクトルが正規直交じゃない場合でもQR分解を使って直交化すればよい*2。

今、両部分空間の元 ${\bf U} {\bf a}$ と ${\bf V} {\bf b}$ のなす角 $\theta$ を考えると

$\displaystyle \cos \theta = \frac{{\bf a}^T {\bf U}^T {\bf V} {\bf b}}{\parallel {\bf U} {\bf a} \parallel_2^{2} \parallel {\bf V} {\bf b} \parallel_2^{2}} = \frac{{\bf a}^T {\bf U}^T {\bf V} {\bf b}}{ ({\bf a}^T {\bf a}) ( {\bf b}^T {\bf b} )}$

ここで、 ${\bf U}$ と ${\bf V}$ に共通部分があるなら、上式の $\cos \theta$ を $1$ （つまりなす角を $0$ ）にできる ${\bf a}, {\bf b}$ が存在するはずであり、それを ${\hat {\bf a}}_1, {\hat {\bf b}}_1$ とすると ${\bf U}{\hat {\bf a}}_1$ と ${\bf V}{\hat {\bf b}}_1$ は共通部分の元である。さらに、 ${\hat {\bf a}}_1, {\hat {\bf b}}_1$ と直交して、かつ、上式 $\cos \theta$ を $1$ にできる ${\hat {\bf a}}_2, {\hat {\bf b}}_2$ が存在するなら、 ${\bf U}{\hat {\bf a}}_2, {\bf V}{\hat {\bf b}}_2$ も共通部分の元である*3。同様にして、互いに直交し上記 $\cos \theta$ を $1$ にできる ${\bf a}, {\bf b}$ が $k$ 組み見つかったとすると、 ${\bf U}{\hat {\bf a}}_1, \cdots, {\bf U}{\hat {\bf a}}_k$ が共通部分の基底となる。このなす角 $\theta$ は正準角（principal angle）と呼ばれ、このブログでも過去に少しだけ触れている（9年越しでタイトルの疑問に回答・・・）。

yamagensakam.hatenablog.com

この記事の内容より、結局のところ ${\bf U}^T {\bf V}$ を特異値分解して、その特異値が $1$ に対応する左特異ベクトルが上記の ${\hat {\bf a}}_1, \cdots, {\hat {\bf a}}_k$ に当たるため、後はそれぞれ左から ${\bf U}$ をかければ良い。

こちらも以下のコードで実験してみた（テストデータの生成関数make_test_matrixは方法1のコードと同じ）。

if __name__ == '__main__':
    d = 15
    m = 5
    #テストデータの生成
    U, V = make_test_matrix(d, m)

    #直交化
    UQ , _ = np.linalg.qr(U)
    VQ , _ = np.linalg.qr(V)

    #特異値分解でcos(theta)と係数a, bを計算
    A_hat, D, B_hat = np.linalg.svd(UQ.T @ VQ)

    #特異値が1の左特異ベクトルを使って共通部分の基底ベクトルを計算
    A_hat = A_hat[:, np.abs(D - 1.0) < 1e-15]
    basis = UQ @ A_hat

    #結果の表示
    print("singular values")
    print(D)
    print("A_hat")
    print(A_hat)
    print("basis of intersection")
    print(basis)

結果

singular values
[1.         1.         1.         0.90457795 0.29262578]
A_hat
[[ 6.44393662e-01 -1.84065423e-01 -7.42210703e-01]
 [ 6.08863932e-01  7.10720557e-01  3.52364871e-01]
 [-2.07087630e-01  3.03916180e-01 -2.55165300e-01]
 [ 4.13710326e-01 -6.07150036e-01  5.09757726e-01]
 [ 2.49800181e-16 -9.54097912e-16  8.53483950e-16]]
basis of intersection
[[-0.37542124  0.14849663  0.16837461]
 [ 0.22018893 -0.18075024 -0.51658548]
 [ 0.00558178  0.36511633  0.39239376]
 [ 0.26960343  0.1627379  -0.16594012]
 [ 0.0955871  -0.22341439 -0.22149631]
 [ 0.13222487  0.16072238 -0.14863213]
 [ 0.10764232  0.333674   -0.21546023]
 [ 0.41641461 -0.06119148  0.37256877]
 [ 0.41790505 -0.09802718  0.32441074]
 [ 0.11692351  0.01010546 -0.1092835 ]
 [-0.04649923  0.11043067  0.03039536]
 [-0.23646815 -0.31480919  0.06458789]
 [-0.10098034 -0.13695922 -0.30038741]
 [ 0.09864014 -0.66361076  0.22550194]
 [ 0.50983159  0.12044854 -0.0594418 ]]

こちらも3つの特異値が1になり、期待通り3次元分の基底が得られた。

*1:ランダムでたまたま線形従属なベクトルが生成されなければ

*2:直交化しなくても、そのまま共通部分の基底を求めることも可能だが、少し計算が煩雑になるのでここでは書かない。このあたりについては、例えば以下の書籍の6章なんかに書いてあるので興味あればそちらを参照。

入門機械学習による異常検知―Rによる実践ガイド

作者:剛, 井手
発売日: 2015/02/19
メディア: 単行本

*3: ${\bf U}$ が直交行列なので ${\hat {\bf a}}_1\perp {\hat {\bf a}}_2$ ならば ${\bf U}{\hat {\bf a}}_1 \perp {\bf U}{\hat {\bf a}}_2$ 。同様に ${\bf V}$ が直交行列なので ${\hat {\bf b}}_1\perp {\hat {\bf b}}_2$ ならば ${\bf V}{\hat {\bf b}}_1 \perp {\bf V}{\hat {\bf b}}_2$ 。

2020-01-28

MCMCでガウス過程分類（ロジスティック回帰編）

ガウス過程ベイズ推論機械学習

2020/2/2追記

ソースコードに大きな間違いを見つけたので修正した。2か所間違いがあってそれらが奇跡的に作用することでそれっぽい結果に見えていた。詳細は実装の節を参照。

はじめに

以下の記事の中で、ガウス過程と機械学習サポートページにあるサンプルコードをもとにコーシー分布を観測モデルに基づくガウス過程回帰を実装（というより微修正）した。

yamagensakam.hatenablog.com

今回はこれを少しだけ変えて、ガウス過程ロジスティック回帰をMCMCで行い２クラス分類問題を解いてみる。

やったこと

といっても、本質的には観測モデル $p(y|f({\bf x}))$ を

$\displaystyle p(y|f({\bf x})) = {\rm Bern}(y| \sigma(f({\bf x}) ) = \left\{ \sigma(f({\bf x}))\right\}^y \left\{1 - \sigma(f({\bf x}))\right\}^{(1 - y)}$

というベルヌーイ分布にするだけ。あとはテスト入力に対する $y$ の予測値は、テスト入力 ${\bar {\bf x}}$ に対する関数出力値の事後分布からのサンプリング値 $f_1({\bar {\bf x}}), \cdots,f_S({\bar {\bf x}})$ を用いて、

$\displaystyle{\bar y} = \frac{1}{S} \sum_{i = 1}^{S} \sigma(f_i({\bar {\bf x}}))$

として近似する。

実装

以下、ソースコード。

2020/2/20追記

冒頭でも書いたが大ポカを犯していた。間違いは

ベルヌーイ分布の対数尤度にマイナスをつけていた。
学習データをプロットするときに、クラスの色を逆にしていた。

という2点。対数尤度にマイナスをつけてしまっていたので予測値としては逆の結果が得られるが、プロットする際クラスラベルを逆にするというミスも犯していたから結果的にそれっぽい出力が得られ気が付かなかった。ということで、修正したソースコードは以下の通り。

import sys
import numpy as np
from pylab import *
from numpy import exp,log
from elliptical import elliptical
from numpy.linalg import cholesky as chol

def gpr_bell(f, param):
    X,y,Kinv = param
    N = X.shape[1]
    return gpr_bell_lik (y[0:N], f[0:N], Kinv)

def gpr_bell_lik(y, f, Kinv):
    log_sigf = log(1.0/(1.0+exp(-f)))

    # 2020/2/2 ベルヌーイ分布の対数尤度にマイナスをつけていたがない方が正しいので修正
    # return -(np.sum(y @ log_sigf) + np.sum((1 - y) @ (1 - log_sigf)))\
    #         - np.dot (f, np.dot(Kinv, f)) / 2            
    return np.sum((y * log_sigf) + ((1 - y) * (1 - log_sigf)))\
            - (f @ Kinv @ f / 2)

def GetDistanceMatrix(X):
    H = np.tile(np.diag((X.T @ X)) , (np.size(X , 1) , 1))
    G = X.T @ X
    DistMat = H - 2 * G + H.T
    return DistMat

def kgauss(tau, sigma):
    return lambda X: tau * exp (-GetDistanceMatrix(X) / sigma)

def kernel_matrix(X, kernel):
    N = X.shape[1]
    eta = 1e-6
    return kernel(X) + eta * np.eye(N)


def gpr_mcmc(X, y, cov_func, iters, remove_num):
    
    N = X.shape[1]
    K_NN = kernel_matrix(X, cov_func)
    K_NNinv = inv(K_NN)

    #こうすると多次元正規分布からサンプリングすることになる
    P = chol (K_NN)
    f = P @ randn(N)
    
    #MCMCサンプルのうち最初のremove_num個は使わない
    S = iters - remove_num
    f_posterior = np.zeros((f.shape[0], S))

    for iter in range(iters):
        #楕円スライスサンプリング
        f,lik = elliptical (f, P, gpr_bell, (X, y, K_NNinv))
        print ('\r[iter %2d]' % (iter + 1),)
        if iter >= remove_num:
            f_posterior[:, iter - remove_num] = f
    return K_NNinv, f_posterior

def predict(X_train, X_test, cov_func, K_NNinv, f_posterior, sampling_num):
    N = X_train.shape[1]
    M = X_test.shape[1]
    XX = np.concatenate([X_train, X_test], axis = 1)

    K = kernel_matrix(XX, cov_func)
    K_NM = K[:N, N:N + M]
    K_MM = K[N:N + M, N:N + M]
    y_test = np.zeros(M)
    S = f_posterior.shape[1]

    for t in range(sampling_num):
        s = np.random.randint(S)
        f_n = f_posterior[:, s]
        mu = K_NM.T@ K_NNinv @ f_n
        sig = K_MM - K_NM.T @ K_NNinv @ K_NM
        P = chol(sig)
        f_m = P @ np.random.randn(M) + mu
        y_test += (1.0/(1.0+np.exp(-f_m)))

    y_test /= sampling_num
    return y_test

def make_data(xmin, xmax, n):
    np.random.seed(2222)
    cnt = 0
    t = 0
    N = n * 9
    X_train = np.zeros((2, N))
    y_train = np.zeros(N)
    center1 =  np.linspace(-(-xmin[0] + xmax[0])/3.5, (-xmin[0] + xmax[0])/3.5, 3)
    center2 =  np.linspace(-(-xmin[1] + xmax[1])/3.5, (-xmin[1] + xmax[1])/3.5, 3)
    for i in range(3):
        for j in range(3):
            y_train[t:t + n] = cnt % 2
            X_train[:, t:t + n] = 1.5 * np.random.randn(2, n) \
                                    + np.array([center1[i], center2[j]])[:, None]
            t += n
            cnt += 1


    #テスト入力点生成
    T1 = np.linspace(xmin[0], xmax[0], 60)
    T2 = np.linspace(xmin[1], xmax[1], 60)
    X_test = np.zeros((2, T1.shape[0] * T2.shape[0]))
    i = 0
    for t1 in T1:
        for t2 in T2:
            X_test[:, i] = np.array([t1, t2])
            i += 1
    return X_train, y_train, X_test
    
def main ():    
    xmin = np.array([-9, -9])
    xmax = np.array([9, 9])
    
    #学習データとテストデータの生成
    n = 20
    X_train, y_train, X_test = make_data(xmin, xmax, n)

    #共分散関数の設定
    cov_func = kgauss(3, 2.5)
    
    #MCMCで事後分布からのサンプリング
    iters = 1000
    remove_num = 200
    K_NNinv, f_posterior = gpr_mcmc(X_train, y_train, cov_func, iters, remove_num)

    #予測値の計算
    sampling_num = 50
    y_test = predict(X_train, X_test, cov_func, K_NNinv, f_posterior, sampling_num)

    #結果の表示
    size = 100
    sc = plt.scatter(X_test[0, :], X_test[1, :], size, y_test, marker = ',', cmap='bwr')
    plt.colorbar(sc)

    # 2020/2/2 プロットの色が逆だったので修正
    # plt.scatter(X_train[0, y_train == 1], X_train[1, y_train == 1] ,
    #             edgecolors="black", marker="o" , color='b')
    # plt.scatter(X_train[0, y_train == 0], X_train[1, y_train == 0] ,
    #             edgecolors="black", marker="o" , color='r')
    plt.scatter(X_train[0, y_train == 1], X_train[1, y_train == 1] ,
                edgecolors="black", marker="o" , color='r')
    plt.scatter(X_train[0, y_train == 0], X_train[1, y_train == 0] ,
                edgecolors="black", marker="o" , color='b')

    plt.axis([xmin[0], xmax[0], xmin[1], xmax[1]])

    plt.show()

if __name__ == "__main__":
    main ()

上で書いたこと以上に前回のコードから変更しているが、本質的な変更は尤度計算がgpr_bell_likになっておりベルヌーイ分布を使った計算になっていること。他は見れば大体わかると思うが、1つだけ補足するとカーネル行列を作るところで以前はfor文を回していたが遅いのでfor文使わないように修正している。詳細はこの記事の『各データ間の距離の2乗を計算』ってとこ参照。

yamagensakam.hatenablog.com

結果

f:id:YamagenSakam:20200202021523p:plain

丸マーカーでプロットしているのが学習データ、これをもとに各座標位置の所属クラスを予測している。ちゃんと非線形なクラス分類ができている（2020/2/2 修正に伴い当然結果も変わる）。

最後に

ロジスティック回帰編と書いているってことは、プロビット回帰編もいつかやる・・・はず。

2020-01-24

ガウス過程回帰で学習データにない任意の入力点に対する関数出力値をMCMCサンプリングする方法

ガウス過程ベイズ推論機械学習

はじめに

表題についてちょっと調べた＆やってみたのでメモ。やりたいことは表題の通りなんだけど、少し数式で表現して整理。

今、学習データとして $\left\{{\bf X}, {\bf y} \right\} = \left\{({\bf x}_1, y_1), \cdots, ({\bf x}_N, y_N)\right\}$ の $N$ 個得られているとする。
入力 ${\bf x}$ と関数 $f$ が与えられたとき、出力 $y$ は $p(y | f({\bf x}))$ に従うとする（ハイパーパラメータを持つ場合もあるが今回は固定値として扱うので表記から省略）。
そして、この関数はガウス過程に従う確率変数とする（は共分散関数）。
- $p(y | f({\bf x}))$ がガウス分布でないなら、関数の事後分布 $p(f|{\bf X}, {\bf y}) \propto \prod_{i = 1}^{N}p( y_i | f({\bf x}_i))p(f)$ は解析的に得られないのでそういう場合は近似するかMCMC等でサンプリングを行う。
やりたいのは、個の任意の未知入力が与えられて、事後分布からサンプリングされた関数の出力を得たい。
- 予測分布からのサンプリングも得たいという話もあるけど今回は対象外。
  - とりあえず、事後分布からの出力値 $f({\bar {\bf x}})$ が得られれば、これから $p(y |f({\bar {\bf x}}))$ が計算できるので、それっぽい $y$ は得られそう。

ようするに、学習入力点 $N$ 個に対する関数 $f$ の出力を並べたベクトルを ${\bf f}_N = \left[f({\bf x})_1 \cdots, f({\bf x}_N)\ \right]^T$ とし、未知入力点 $M$ 個に対する関数 $f$ の出力を並べたベクトルを ${\bar {\bf f}}_M = \left[f({\bar {\bf x}})_1 \cdots, f({\bar {\bf x}}_M)\ \right]^T$ とした場合、

$\displaystyle p(f|{\bf X}, {\bf y}) = p({\bf f}_N | {\bf y}) \propto p({\bf y}| {\bf f}_N) p({\bf f}_N)$

からのMCMCサンプリングはできるのは分かるけど、 $p({\bar {\bf f}}_M | {\bf y})$ からのサンプリングはどうやったらいいんだ？って話。

方法1 ~MCMCサンプリングの時、同時分布 $p({\bf f}_N, {\bar {\bf f}}_M )$ からサンプリング~

非常にお世話になっている書籍

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

作者:持橋大地,大羽成征
出版社/メーカー: 講談社
発売日: 2019/03/09
メディア: 単行本（ソフトカバー）

のサポートページに、コーシー分布を観測モデルとしたガウス過程回帰のpythonコードgpr-cauchy.py, elliptical.pyがある。これは楕円スライスサンプリングというMCMCサンプリングの一手法を用いて関数出力値の推定を行っているが、これを見るとサンプルを得るところで

$p({\bf f}_N, {\bar {\bf f}}_M ) = {\mathcal N} \left( {\bf f}_N, {\bar {\bf f}}_M| 0, \begin{pmatrix} {\bf K}_{NN} & {\bf K}_{NM} \\ {\bf K}_{NM}^T & {\bf K}_{MM} \\ \end{pmatrix} \right)$

（ただし、 ${\bf K}_{NN}$ は $K_{i,j} = , k({\bf x}_i, {\bf x}_j) \; (i,j \in \left\{1, \cdots, N \right\})$ という要素を持つ $N \times N$ の行列、 ${\bf K}_{NM}$ は $K_{i,j} = , k({\bf x}_i, {\bar {\bf x}}_j) \; (i \in \left\{1, \cdots, N \right\}, j \in \left\{1, \cdots, M \right\})$ という要素を持つ $N \times M$ の行列、 ${\bf K}_{MM}$ は $K_{i,j} = , k({\bar{\bf x}}_i, {\bar{\bf x}}_j) \; (i,j \in \left\{1, \cdots, M \right\})$ という要素を持つ $M \times M$ の行列）

からのサンプルを使って更新している。そして、そのサンプルを採択するか否かの尤度計算は $p({\bf y}| {\bf f}_N) p({\bf f}_N)$ で計算、つまり学習入力点だけで行っている*1。

なるほど、確かにこうすれば $p({\bf f}_N, {\bar {\bf f}}_M | {\bf y})$ からサンプリングしたことになり、そのサンプルから ${\bf f}_N$ を捨てれば周辺化することになるから $p({\bar {\bf f}}_M | {\bf y})$ からサンプリングしたことになる。

方法2 ~ $p({\bf f}_N | {\bf y})$ からのサンプルを用いて $p({\bar {\bf f}}_M | {\bf y})$ からサンプリング~

方法1だと、学習時に評価したい未知入力 $\left\{{\bar {\bf x}}_1, \cdots, {\bar{\bf x}}_M \right\}$ がわかっている必要があるが、場合によってはそれが困難だったりする。そういう場合はどうするか？

$p({\bf f}_N | {\bf y})$ からのサンプルが $S$ 個 ${\bf f}_N^{(1)}, \cdots, {\bf f}_N^{(S)}$ 得られているとすると、

$\displaystyle p({\bar {\bf f}}_M | {\bf y}) = \int p({\bar {\bf f}}_M | {\bf f}_N) p({\bf f}_N | {\bf y})d {\bf f} _N \sim \frac{1}{S} \sum_{s = 1}^{S}p({\bar {\bf f}}_M | {\bf f}_N^{(s)})$

と近似できる。ここでガウス過程の性質より、

$p({\bar {\bf f}}_M | {\bf f}_N^{(s)}) = {\mathcal N} ({\bar {\bf f}}_M | {\bf K}_{NM}^T {\bf K}_{NN}^{-1} {\bf f}_N^{(s)} , {\bf K}_{MM} - {\bf K}_{NM}^T {\bf K}_{NN}^{-1} {\bf K}_{NM} )$

なので、 $p({\bar {\bf f}}_M | {\bf y})$ は全て同じ混合比の混合正規分布で近似できる。従って $s \in \left\{1 \cdots S \right\}$ の乱数を発生させた後、それに紐づく多次元正規分布 ${\mathcal N} ({\bar {\bf f}}_M | {\bf K}_{NM}^T {\bf K}_{NN}^{-1} {\bf f}_N^{(s)} , {\bf K}_{MM} - {\bf K}_{NM}^T {\bf K}_{NN}^{-1} {\bf K}_{NM} )$ からサンプリングすればよい。

実装

方法2によるコーシー分布ガウス過程回帰を上記サポートページにあるgpr-cauchy.pyを改変する形でやってみた。

import sys
import putil
import numpy as np
from pylab import *
from numpy import exp,log
from elliptical import elliptical
from numpy.linalg import cholesky as chol

def gpr_cauchy (f, param):
    x,y,gamma,Kinv = param
    M = len(x)
    return gpr_cauchy_lik (y[0:M], f[0:M], gamma, Kinv)

def gpr_cauchy_lik (y, f, gamma, Kinv):
    return - np.sum (log (gamma + (y - f)**2 / gamma)) \
           - np.dot (f, np.dot(Kinv, f)) / 2

def kgauss (tau,sigma):
    return lambda x,y: exp(tau) * exp (-(x - y)**2 / exp(sigma))

def kernel_matrix (xx, kernel):
    N = len(xx)
    eta = 1e-6
    return np.array (
        [kernel (xi, xj) for xi in xx for xj in xx]
    ).reshape(N,N) + eta * np.eye(N)

def gpr_mcmc(x, y, iters, xmin, xmax, gamma):
    
    N = len(x)
    K_NN = kernel_matrix(x, kgauss(1,1))
    K_NNinv = inv(K_NN)

    #こうすると多次元正規分布からサンプリングすることになる
    P = chol (K_NN)
    f = P @ randn(N)
    
    #MCMCサンプルのうち最初の200個は使わない
    remove_num = 200
    S = iters - remove_num
    f_posterior = np.zeros((f.shape[0], S))
    
    for iter in range(iters):
        #楕円スライスサンプリング
        f,lik = elliptical (f, P, gpr_cauchy, (x,y,gamma,K_NNinv))
        print ('\r[iter %2d]' % (iter + 1),)
        if iter >= remove_num:
            f_posterior[:, iter - remove_num] = f
    
    #未知入力点
    x_new = np.linspace(xmin, xmax, 100)
    M = len(x_new)
    xx = np.hstack((x, x_new))

    K = kernel_matrix(xx, kgauss(1, 1))
    K_NM = K[:N, N:N + M]
    K_MM = K[N:N + M, N:N + M]

    g = np.zeros (len(x_new))
    #50個サンプリング
    sampling_num = 50
    for t in range(sampling_num):
        s = np.random.randint(S)
        f_n = f_posterior[:, s]
        mu = K_NM.T@ K_NNinv @ f_n
        sig = K_MM - K_NM.T @ K_NNinv @ K_NM
        P = chol(sig)
        f_m = P @ np.random.randn(M) + mu
        g = g + f_m
        plot (x_new, f_m)
    plot (x, y, 'bx', markersize=14)
    plot (x_new, g/sampling_num, 'k', linewidth=3)

def usage ():
    print ('usage: gpr-cauchy.py data.xyf iters [output]')
    sys.exit (0)

def main ():
    xmin = -5
    xmax =  5
    ymin = -7.5
    ymax = 12.5
    gamma = 0.2
    
    if len(sys.argv) < 3:
        usage ()
    else:
        [x, y, f] = np.loadtxt (sys.argv[1]).T
        iters = int (sys.argv[2])

    gpr_mcmc (x, y, iters, xmin, xmax, gamma)
    axis ([xmin, xmax, ymin, ymax])
    show ()

if __name__ == "__main__":
    main ()

楕円スライスサンプリングを行う関数ellipticalについては、サポートページのelliptical.pyから変えてないので省略。また多次元正規分布からはコレスキー分解を用いてサンプリングしている（この辺が参考になる）。これを、

python .\gpr-cauchy_new.py  .\gpr-cauchy.dat 1000

として、実行すると以下のような結果が得られた。

f:id:YamagenSakam:20200123235146p:plain

うん、いい感じ。

今後やりたいこと

ハイパーパラメータも事前分布与えてMCMCサンプリング
予測分布からのサンプリング
ガウス過程分類
変分近似

*1:事前分布の方は $p({\bf f}_N, {\bar {\bf f}}_M)$ の同時分布を評価するべきかと思ったがどうなんだろ？実際両方やってみたけど、結果はあまり変わらなかった

2020-01-06

EMアルゴリズムと変分ベイズEMアルゴリズムのメモ

機械学習ベイズ推論雑メモ

表題の通り。この辺りややこしくていつもわからなくなるので、身につけるためにも自分なりの理解をまとめたメモを書く。

対数尤度と下限とKLダイバージェンス

まず、既知の変数を $X$ 、未知の変数（ここでは潜在変数か未知パラメータかは区別していない）を $Z$ とすると以下が成り立つ。

$\displaystyle \ln p(X) = KL\left[q(Z)\parallel p(Z| X) \right] + {\mathcal L}\left[q(Z)\right] \tag{1}$

ただし、 $KL\left[q(Z) \parallel p(Z | X) \right]$ はKLダイバージェンスで、

$\displaystyle KL\left[q(Z) \parallel p(Z | X) \right] =\int q(Z) \ln \frac{q(Z)}{p(Z | X)} dZ = - \int q(Z) \ln \frac{p(Z | X)}{q(Z)} dZ \tag{2}$

${\mathcal L}\left[q(Z)\right]$ は、対数尤度の下限で

$\displaystyle {\mathcal L}\left[q(Z)\right] = \int q(Z) \ln \frac{p(X, Z)}{q(Z)} dZ \tag{3}$

であり、常に

$\displaystyle \ln p(X) \geq {\mathcal L}\left[q(Z)\right] \tag{4}$

が成り立つ。これは関数 $\ln$ が上に凸でありイェンセン不等式より $\int p(x)\ln f(x)dx \leq \ln \int p(x)f(x)dx$ であるため、 $\int q(Z) \ln \frac{p(X, Z)}{q(Z)} dZ \leq \ln \int q(Z)\frac{p(X, Z)}{q(Z)} dZ =\ln \int p(X, Z) dZ = \ln p(X)$ として得られる。

また、式（1）は、式（2）と式（3）を実際に足してみると得られる。

EMアルゴリズム

考え方

得られているものは既知の変数の $X$ 、求めたいのは未知パラメータ $\Theta$ 。
目的関数は対数尤度であり、これを最大化するを求めたい。
- ところが、これを直接的に最適化するのは一般に困難。
ここで、潜在変数を導入し、という量を考える。
- 式（4）より、 $\ln p(X|\Theta) \geq {\mathcal L}\left[q(Z), \Theta \right]$ が成り立つ。
そのため、を最大化する関数とパラメータを求めることを考える。
- ただし、この両者を一気に最適化する解析解は得られないので、座標降下法の考えで一方を固定し他方を最適化するという操作を繰り返して求める。

EMアルゴリズムの更新式

ということで、更新アルゴリズムは以下の通り。

まずは、パラメータを既知の値と固定した場合に、を最大にする関数を求める。
- 式（1）より、 $\displaystyle {\mathcal L}\left[q(Z), {\hat \Theta} \right] = \ln p(X| {\hat \Theta}) - KL\left[q(Z) | p(Z| X, {\hat \Theta}) \right]$ 。
- $KL\left[q(Z) \parallel p(Z| X, {\hat \Theta})\right] \geq 0$ なので、結局 $q(Z) = p(Z| X, {\hat \Theta})$ とするのが最適。
- 色々分け方があるようだが、これをEステップと呼ぶことが多いらしい。
次に、関数を上で求めたに固定し、を最大化するを求める。
- ${\mathcal L} \left[p(Z| X, {\hat \Theta}), \Theta \right] = \int p(Z| X, {\hat \Theta})\ln \frac{p(X, Z|\Theta)}{p(Z| X, {\hat \Theta})} dZ = \int p(Z| X, {\hat \Theta})\ln p(X, Z|\Theta) dZ + {\rm const}$ なので、 $\int p(Z| X, {\hat \Theta})\ln p(X, Z|\Theta) dZ = \langle \ln p(X, Z|\Theta)\rangle _{p(Z| X, {\hat \Theta})}$ を最大化する $\Theta$ を求めればよい。
- 期待値の最適化なので $\ln p(X, Z|\Theta)$ さえ最適化できるなら適用可能。
- こちらをMステップと呼ぶ。

変分ベイズ EMアルゴリズム

変分近似の考え方

得られているものは既知の変数の $X$ 、未知変数は潜在変数もパラメータもひっくるめて $Z$ とする。
モデルが与えられていて、未知変数の事後分布を求めたい。
- これは、 $p(Z | X) = \frac{p(X, Z)}{p(X)}$ であり、厳密に求めるのは困難（ $p(X)$ が難しい）。
そこで、を別の分布で近似することを考える。
- 当然 $q(Z)$ は $p(Z | X)$ に近い方がいいので、式（2）のKLダイバージェンスを最小化する関数 $q$ を求めることになる。
その際、に何の制約も課さないと、という無意味な解が得られるので、に何らの制約を課して表現能力を限定する。
- 典型的なのは、 $q(Z) = \prod_{i=1}^{D} q(z_i | X)$ のように未知変数間の独立性を制約にする近似（平均場近似）。
式（1）より、KLダイバージェンスを最小化はを最大化することと等価。
- $\ln p(X)$ は定数なので。
ということで、定めた制約のもと ${\mathcal L}\left[q(Z)\right]$ を最大化する $q$ を求めればよい。

ここまでは変分近似の考え方。おそらく、変分ベイズ EMアルゴリズムと言うと、以下のようなもう少し限定されたアルゴリズムのことを言うと思う。

変分EMアルゴリズムの考え方

今、未知変数がパラメータ $\Theta$ と潜在変数 $Z$ だとする。
これらの事後分布を求めたいが厳密に求めるのは困難なので、という分布で近似することを考える。
- $q$ の制約としては、 $q(\Theta, Z) = q(\Theta)q(Z)$ という独立性を制約とする。
- モデルに応じてどういう独立性を仮定するのが適切なのかが変わってくる（ $\Theta$ をいくつかの要素が独立だと仮定したり）。
式（3）より、。
- この式はベイズの定理よりと展開している。
  - これはベイズの定理から導かれる展開で、モデル次第ではもっと色々展開できる可能性がある（ $p(Z|\Theta)=p(Z)$ とできたりとか）。
後はこれを座標降下法の要領で、 $q(\Theta)$ を固定して $q(Z)$ に関する最適化と $q(Z)$ を固定して $q(\Theta)$ に関する最適化を交互に繰り返せばよい。

変分EMアルゴリズムの更新式の導出

ここからはもう少し、この式を展開する。

$\begin{align} {\mathcal L} \left[q(\Theta)q(Z) \right] &= \int \int q(\Theta)q(Z) \ln \frac{p(X|\Theta, Z)p(Z | \Theta)p(\Theta)}{q(\Theta)q(Z)} dZ d\Theta \\ &= \int \int q(\Theta)q(Z) \ln p(X|\Theta, Z) dZ d\Theta + \int \int q(\Theta)q(Z)\ln p(Z | \Theta) dZ d\Theta \\ &\quad +\int \int q(\Theta) \ln \frac{p(\Theta)}{q(\Theta)}d\Theta + \int q(Z)\ln\frac{1}{ q(Z)} dZ \end{align}$

この式において、まず $q(\Theta) = {\hat q(\Theta)}$ という既知の関数として、 $q(Z)$ について最大化することを考える。

$\displaystyle \begin{align} {\mathcal L} \left[{\hat q}(\Theta) q(Z) \right] &=\int \int {\hat q}(\Theta)q(Z) \ln p(X|\Theta, Z) dZ d\Theta + \int \int {\hat q}(\Theta)q(Z)\ln p(Z | \Theta) dZ d\Theta \\ &\quad + \int q(Z)\ln\frac{1}{ q(Z)} dZ + {\rm const} \\ &=\int q(Z) \langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)} dZ + \int q(Z)\langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)} dZ \\ &\quad + \int q(Z)\ln\frac{1}{ q(Z)} dZ + {\rm const} \\ &=\int q(Z) \ln \exp(\langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)}) dZ + \int q(Z)\ln \exp(\langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)}) dZ \\ &\quad + \int q(Z)\ln\frac{1}{ q(Z)} dZ + {\rm const} \\ &= \int q(Z)\ln\frac{\exp \left(\langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)} + \langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)} \right)}{ q(Z)} dZ + {\rm const} \end{align}$

この最後の式をよく見ると、分子を $Z$ で積分した時 $1$ になれば（つまり確率分布であれば）、これはKLダイバージェンス（に $-1$ をかけたもの）である。そのため、

$\displaystyle r(Z) \propto \exp \left(\langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)} + \langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)} \right)$

$\displaystyle \int r(Z) dZ = 1$

と定義すれば、

$\displaystyle {\mathcal L} \left[{\hat q}(\Theta) q(Z) \right] = -KL\left[ q(Z) \parallel r(Z) \right] + {\rm const}$

になり、下限を最大化するためにはKLダイバージェンスを0にする $q(Z)$ を求めればよい。すなわち、

$\displaystyle q(Z) = r(Z) \propto exp \left(\langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)} + \langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)} \right)$

とすればよい。この潜在変数側の更新を変分Eステップという。ちなみに両辺logを取ると、

$\displaystyle \ln q(Z) = \langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)} + \langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)} + {\rm const}$

と書ける。 $q(Z)$ がどんな分布かを求める際、このlogを取った形式で計算をすることが多い。

同様の式展開を $q(Z) = {\hat q}(Z)$ という既知の関数として $q(\Theta)$ について最大化する場合に適用すると、

$\displaystyle \begin{align} {\mathcal L} \left[q(\Theta) {\hat q}(Z) \right] = \int q(\Theta)\ln\frac{\exp \left(\langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(Z)} + \langle \ln p(Z | \Theta) \rangle_{ {\hat q}(Z)} \right) p(\Theta)}{ q(\Theta)} d\Theta + {\rm const} \end{align}$

なので、 $q(\Theta)$ に関する更新式は、

$\displaystyle q(\Theta) \propto \exp \left( \langle \ln p(X | \Theta, Z) \rangle_{{\hat q}(Z)} + \langle \ln p(Z | \Theta) \rangle_{{\hat q}(Z)} \right) p(\Theta)$

となる。こちらは変分Mステップと呼ばれる。こちらも両辺logを取った式を書いておくと、

$\displaystyle \ln q(\Theta) = \langle \ln p(X | \Theta, Z) \rangle_{{\hat q}(Z)} + \langle \ln p(Z | \Theta) \rangle_{{\hat q}(Z)} + \ln p(\Theta) + {\rm const}$

2019-02-10

ベイズ推論による教師あり分類～ガウス混合分布編～

機械学習ベイズ推論

はじめに

以前、ベイズ推論による機械学習入門という書籍を読み、以下のような記事を公開しました。

yamagensakam.hatenablog.com

この書籍には、多次元ガウス混合モデルによるクラスタリング手法が詳しく述べられていますが、多次元ガウス混合モデルによる教師あり分類の手法は少し触れる程度にとどまっています（クラスタリングで述べられている導出過程を流用すれば、比較的簡単に導出できるためだと思います）。
ということで、今回は多次元ガウス混合モデルによる教師あり分類の導出、実装をやってみました。

ガウス混合モデル

モデル自体は教師ありだろうが、教師なしだろうが同じモデルを考えることができます。まず、全て独立でサンプリングされる $D$ 次元のデータが ${\bf X} = \left\{ {\bf x}_1 , \cdots, {\bf x}_N \right\}$ と $N$ 個あるとし、そのそれぞれに対応するラベルデータを ${\bf S} = \left\{ {\bf s}_1 , \cdots, {\bf s}_N \right\}$ とします。ここで、 ${\bf s}$ は $1$ of $K$ 表現と呼ばれ、クラス数が $K$ の場合、所属するクラスに対応する要素のみ $1$ で他は全て $0$ の $K$ 次元ベクトルです。この ${\bf s}$ をサンプルするための分布として、パラメータ ${\bf \pi}$ を持つカテゴリ分布を考えます。
また、 $K$ 個のガウス分布が混合したモデルを考えるので、それぞれの分布の平均を ${\bf \mu}_1, \cdots, {\bf \mu}_K$ 、精度行列（共分散行列の逆行列）を ${\bf \Lambda}_1,\cdots, {\bf \Lambda}_K$ とします。なお、各パラメータを個別に扱わない場面では表記を簡単にするため、これらのガウス混合分布のパラメータを全てまとめて ${\bf \Theta}$ と表記します。
更に、カテゴリ分布のパラメータ ${\bf \pi}$ や各ガウス分布の ${\bf \mu}_k$ 、 ${\bf \Lambda}_k$ といったパラメータもベイズ推論の枠組みでは、確率分布に従うと考えます。ここでは、それぞれ共役事前分布に従うことにします。具体的には、

$\displaystyle p({\bf \pi}) = {\rm dir}({\bf \pi} | {\bf \alpha}) \tag{1}$
$\displaystyle p({\bf \mu}_k, {\bf \Lambda}_k) = {\rm NW} ({\bf \mu}_k, {\bf \Lambda}_k | {\bf m}, \beta, \nu, {\bf W} ) \tag{2}$

という分布に従うとします。ここで、 ${\rm dir}(\cdot)$ はディリクレ分布、 ${\rm NW}(\cdot)$ はガウス・ウィシャート分布です（具体的な分布の式は書籍を参照ください）。ここで事前分布として共役事前分布を採用しておくと事後分布や予測分布が解析的に求めることができます。 ${\bf \alpha}$ や ${\bf m}, \beta, \nu, {\bf W}$ は超パラメータで今回は事前に与えるものとします。また、 ${\bf m}, \beta, \nu, {\bf W}$ については、それぞれのガウス分布で異なる値を与えることもできますが、簡単のため今回は全て共通の値とします。
さて、ここまで出てきた確率変数 ${\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta}$ の同時分布を考えると、以下のようになります。
${ \displaystyle \begin{eqnarray} p({\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta}) &=& p({\bf X} | {\bf S}, {\bf \Theta}) p ({\bf S}|{\bf \pi})p({\bf \Theta})p({\bf \pi} ) \\ &=& \left\{\prod_{i = 1}^{N} \prod_{k = 1}^{K} {\rm N}({\bf x}_i| {\bf \mu}_k, {\bf \Lambda}_k^{-1})^{{\bf s}_{i, k}} {\rm Cat}({\bf s}_k | {\bf \pi})\right\} \left\{ \prod_{k = 1}^{K} p({\bf \mu}_k, {\bf \Lambda}_k) \right\} p({\bf \pi}) \end{eqnarray} } \tag{3}$

ここで、 ${\rm N}(\cdot)$ はガウス分布、 ${\rm Cat}(\cdot)$ はカテゴリ分布を表し、 $p({\bf \pi})$ と $p({\bf \mu}_k, {\bf \Lambda}_k)$ は式（1）、式（2）で表される分布です。
このモデルのグラフィカルモデルを下図に示します。この図で ${\bf X}$ と ${\bf S}$ の四角で囲まれているのが学習データとして与えられるデータ、 ${\bf x}_*$ と ${\bf s}_*$ がクラスを予測する対象のデータ（テストデータ）です。

このように図示しておくと有効分離の考え方を用いて、変数間の条件付き独立性がわかりやすくなり、後々式展開する際に便利です。有効分離については、以下の記事がわかりやすいと思います。

条件付き独立と有向分離を用いた統計モデルの妥当性チェック - StatModeling Memorandum

推論

上述式（3）のモデルから、パラメータの事後分布とそれを元に、新たなデータ ${\bf x}_*$ が与えられた時の ${\bf s}_*$ の確率分布である予測分布を導出します。このモデルの場合、上述のように共役事前分布を事前分布として採用しているため、事後分布および予測分布は解析的に求めることができます。ここでは、逐次的にデータが与えられるとして、オンラインで予測分布を更新していけるような式を導出します。

事後分布の導出

今回扱う問題では ${\bf X}$ と ${\bf S}$ がデータとして与えられるので、これらが与えられている下でのパラメータの事後分布を求めます。つまり、

$\displaystyle p({\bf \Theta}, {\bf \pi}|{\bf X}, {\bf S}) = \frac{p({\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta})}{p({\bf X}, {\bf S})} \varpropto p({\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta}) \tag{4}$

を求めることになります。なお、クラスタリングの場合、 ${\bf S}$ がデータとして与えられているわけではないので、求めるべき事後分布としては $p({\bf \Theta}, {\bf \pi}, {\bf S}|{\bf X})$ となり、こちらは解析的に求めることが難しくなります（なので、ギブスサンプリングや変分推論などの近似手法を使う）。
ここで式（3）を見ると、 $\Theta$ と $\pi$ は完全に分離していることがわかります。つまり、 ${\bf X}$ と ${\bf S}$ が与えられた下では、 $\Theta$ と $\pi$ は独立であり、

$\displaystyle p({\bf \Theta} |{\bf X}, {\bf S}) \varpropto p({\bf X} | {\bf S}, {\bf \Theta})p({\bf \Theta}) \tag{5}$
$\displaystyle p({\bf \pi} |{\bf X}, {\bf S}) =p({\bf \pi} |{\bf S}) \varpropto p ({\bf S}|{\bf \pi})p({\bf \pi} ) \tag{6}$

を求めればよいことになります。
それでは、式（5）を展開してみます。式（5）の右辺について、対数をとってみると、

${ \displaystyle \begin{eqnarray} \ln p({\bf X} | {\bf S}, {\bf \Theta}) p({\bf \Theta}) &=& \ln \prod_{i = 1}^{N} \prod_{k = 1}^{K} {\rm N}({\bf x}_i| {\bf \mu}_k, {\bf \Lambda}_k^{-1})^{{\bf s}_{i, k}} p({\bf \mu}_k, {\bf \Lambda}_k)\\ &=& \sum_{k = 1}^{K} \left\{ \sum_{i = 1}^{N} s_{i,k} \ln {\rm N}({\bf x}_i| {\bf \mu}_k, {\bf \Lambda}_k^{-1}) + {\rm NW} ({\bf \mu}_k, {\bf \Lambda}_k | {\bf m}, \beta, \nu, {\bf W} ) \right\} \end{eqnarray} } \tag{7}$

となります。これを見ると各 $k$ に対して分離しているのがわかるので、各 $k$ に対し ${\bf \mu}_k, {\bf \Lambda}_k$ それぞれ独立に事後分布を求めればよいことになります。また、式（7）の中括弧の中はガウス・ウィシャート分布を事前分布として用いたガウス分布パラメータの事後分布推定の計算となります。この計算方法は件の書籍に詳しく載っているので、ここでは割愛しますが、最終的な形は以下のようになります。
$\displaystyle p({\bf \mu}_k, {\bf \Lambda}_k | {\bf X}, {\bf S}) = {\rm WA}({\bf \mu}_k, {\bf \Lambda}_k | \hat{{\bf m}}_k, {\hat \beta}_k, {\hat \nu}_k, {\hat {\bf W}}_k) \tag{8}$
ただし、

${ \displaystyle \begin{eqnarray} {\hat \beta}_{k} &=& \sum_{i = 1}^{N} s_{i, k} + \beta \\ {\hat \nu}_{k} &=& \sum_{i = 1}^{N} s_{i, k} + \nu \\ {\hat {\bf m}}_{k} &=& \frac{\sum_{i = 1}^{N} s_{i, k} {\bf x}_{i} + \beta {\bf m} } {{\hat \beta}_{k}} \\ {\hat {\bf W}}_{k} ^{-1} &=& \sum_{i = 1}^{N} s_{i, k} {\bf x}_{i} {\bf x}_{i}^T + \beta {\bf m}{\bf m}^T - {\hat \beta}_{k}{\hat {\bf m}}_{k} {\hat {\bf m}}_{k}^{T} + {\bf W}^{-1}\\ \end{eqnarray} } \tag{9}$

です。同様に式（6）はディリクレ分布を事前分布としてカテゴリ分布の事後分布そのものであるため、

$\displaystyle p({\bf \pi} |{\bf S}) = {\rm dir}({\bf \pi}| {\hat {\bf \alpha}}) \tag{10}$

というディリクレ分布になります。ただし、

$\displaystyle \alpha_{k} = \sum_{i = 1}^{N} s_{i, k} + \alpha \tag{11}$

です。

ここで、式（9）および式（11）の各超パラメータをデータが与えられるたびに逐次更新できる形に変更すると、

${ \displaystyle \begin{eqnarray} {\hat \beta}_{n, k} &=& {\hat \beta}_{n-1, k} + s_{n, k} \\ {\hat \nu}_{n, k} &=& {\hat \nu}_{n-1, k} + s_{i, k} \\ {\hat {\bf m}}_{n, k} &=& \frac{s_{n, k} {\bf x}_{n} + {\hat \beta}_{n-1, k} {\hat {\bf m}}_{n-1, k} } {{\hat \beta}_{n, k}}\\ {\hat {\bf W}}_{n, k}^{-1} &=& { \hat {\bf W}}_{n-1, k}^{-1} + s_{i, k}{\bf x}_{i} {\bf x}_{i}^T +{\hat \beta}_{n-1, k}{\hat {\bf m}}_{n-1, k} {\hat {\bf m}}_{n-1, k}^{T} - {\hat \beta}_{n, k}{\hat {\bf m}}_{n, k} {\hat {\bf m}}_{n, k}^{T} \\ {\hat \alpha}_{n, k} &=& {\hat \alpha}_{n-1, k} + s_{i, k} \end{eqnarray} } \tag{12}$

となり、これにより $n-1$ の時点で、計算した超パラメータを用いて、 $n$ 番目のデータが来たときの超パラメータに更新することが可能です（初期値は事前分布の超パラメータ ${\bf m}, \beta, \nu, {\bf W}, \alpha$ ）。

予測分布の導出

最終的に欲しいのは、新たなデータ ${\bf x}_*$ が与えられた時のクラス分類結果 ${\bf s}_*$ です。そのため、 ${\bf x}_*$ が与えられた下での ${\bf s}_*$ の予測分布を求めます。予測分布は、

$\displaystyle p({\bf s}_* | {\bf x}_*, {\bf X}, {\bf S} ) = \frac{p({\bf s}_* , {\bf x}_*, {\bf X}, {\bf S} )}{p( {\bf x}_*, {\bf X}, {\bf S} )} = \frac{p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S} )p( {\bf X}|{\bf s}_*, {\bf S})p({\bf s}_*|{\bf S})p({\bf S})}{{p( {\bf x}_*, {\bf X}, {\bf S} )}} \tag{13}$

となります。ここで、 ${\bf s}_*$ に関係のない項は無視していいので、 $p({\bf S})$ と $p( {\bf x}_*, {\bf X}, {\bf S} )$ は考慮から外します。更に、 $p( {\bf X}|{\bf s}_*, {\bf S})$ について、グラフィカルモデルから有効分離性を読み解くと ${\bf S}$ が与えられた下では、 ${\bf X}$ と ${\bf s}_*$ は独立であり（ ${\bf X}$ に繋がる ${\bf S}$ が観測され、かつ、head-to-head型を形成している ${\bf s}_*, {\bf x}_*, {\bf \Theta}$ において ${\bf x}_*$ が観測されていないため）、 $p( {\bf X}|{\bf s}_*, {\bf S})= p( {\bf X}| {\bf S})$ となるため、こちらも考慮しなくてよくなります。まとめると、式（13）は

$\displaystyle p({\bf s}_* | {\bf x}_*, {\bf X}, {\bf S} ) \varpropto p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S} ) p({\bf s}_*|{\bf S}) \tag{14}$

と、2つの項のみ考えればよくなります。1つ目の項は、
$\displaystyle p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S} ) = \int p({\bf x}_* | {\bf s}_* ,{\bf X},{\bf S}, {\bf \Theta} ) p({\bf \Theta}| {\bf s}_*,{\bf X},{\bf S}) d {\bf \Theta} \tag{15}$

であり、ここでもグラフィカルモデルから有効分離性を考えると $p({\bf x}_* | {\bf s}_* ,{\bf X},{\bf S}, {\bf \Theta} ) =p({\bf x}_* | {\bf s}_* , {\bf \Theta} )$ 、 $p({\bf \Theta}| {\bf s}_*,{\bf X},{\bf S}) = p({\bf \Theta}| {\bf X},{\bf S})$ とすることができます。従って、

$\displaystyle p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S} ) = \int p({\bf x}_* | {\bf s}_* , {\bf \Theta} ) p({\bf \Theta}| {\bf X},{\bf S}) d {\bf \Theta} \tag{16}$

となります。
$p({\bf \Theta}| {\bf X},{\bf S})$ は、式（5）、式（8）、式（9）に記している事後分布そのものです。今、各 $k$ に対して ${s}_{*,k} = 1$ であるそれぞれの予測分布を求めたいので、式（16）より

$\displaystyle p({\bf x}_*| {\bf s}_{*,k} = 1 , {\bf X}, {\bf S} ) = \int \int p({\bf x}_* | {\bf \mu}_k, {\bf \Lambda}_k ) p({\bf \mu}_k, {\bf \Lambda}_k | {\bf X},{\bf S}) d {\bf \mu}_k d {\bf \Lambda}_k \tag{17}$

を求めればよいことになります。これは、ガウス・ウィッシャート分布を事前分布として用いた場合のガウス分布の予測分布計算になります。こちらも予測分布導出は書籍に載っているので、導出過程は省略し最終形を示します。

$\displaystyle p({\bf x}_*| {\bf s}_{*,k} = 1 , {\bf X}, {\bf S} ) = {\rm St}({\bf x}_* |{\hat {\bf m}}_{n, k}, \frac{(1-D+ {\hat \nu}_{n, k}){\hat \beta}_{n, k} }{1 + {\hat \beta}_{n, k} } {\hat {\bf W}}_{n,k}, 1-D+ {\hat \nu}_{n, k}) \tag{18}$

ただし、 ${\rm St}(\cdot)$ は多次元のスチューデントのt分布で、 $D$ は ${\bf x}$ の次元です。
次に、式（14）の $p({\bf s}_*|{\bf S})$ ですが、こちらも同様に、

$\displaystyle p({\bf s}_*|{\bf S}) = \int p({\bf s}_*|{\bf S},{\bf \pi})p({\bf \pi}| {\bf S}) d {\bf \pi} = \int p({\bf s}_*|{\bf \pi})p({\bf \pi}| {\bf S}) d {\bf \pi} \tag{19}$

となります。ただし、第2式から第3式への展開はグラフィカルモデルより ${\bf \pi}$ を観測した下での ${\bf S}$ と ${\bf s}_*$ の条件付き独立性を利用しました。そしてこの式（19）は、ディリクレ分布を事前分布として使ったカテゴリ分布の予測分布です。ということで、この予測分布についても導出過程は割愛し最終形を示します。

$\displaystyle p({\bf s}_*|{\bf S}) = {\rm Cat}({\bf s}_*|{\hat {\bf \eta}}_{n}) \tag{20}$

ただし、 ${\hat {\bf \eta}_{n}}$ の各要素 ${\hat \eta}_{n,k}$ は
$\displaystyle {\hat \eta}_{n,k}= \frac{{\hat \alpha}_{n, k}}{\sum_{i=1}^K {\hat \alpha}_{n, i}} \tag{21}$
です。
ここまでの話をまとめると、 $n$ 番目のデータが与えら手た時の予測分布の更新式は、式（14）、式（17）、式（20）、式（21）より

$\displaystyle p(s_{*, k} = 1 | {\bf x}_*, {\bf X}, {\bf S} ) \varpropto {\hat \eta}_{n,k} {\rm St}({\bf x}_* |{\hat {\bf m}}_{n, k}, \frac{(1-D+ {\hat \nu}_{n, k}){\hat \beta}_{n, k} }{1 + {\hat \beta}_{n, k} } {\hat {\bf W}}_{n,k}, 1-D+ {\hat \nu}_{n, k}) \tag{22}$

となります（各超パラメータは式（12）、式（21））。そして、各クラスの所属確率の平均を得たいときは、 $k = 1 \cdots K$ に対し、式（22）の右辺に ${\bf x}_*$ を代入し、得られる値を合計が1になるように正規化すればよいです。

実装

ということで実装です。

import numpy as np
import scipy as sci
import math
import random
import matplotlib.pyplot as plt
import matplotlib.animation as anm

class multivariate_student_t():
    def __init__(self, mu, lam, nu):
        self.D   = mu.shape[0]
        self.mu  = mu
        self.lam = lam
        self.nu  = nu

    def pdf(self, x):
        temp1 = np.exp( math.lgamma((self.nu+self.D)/2) - math.lgamma(self.nu/2) )
        temp2 = np.sqrt(np.linalg.det(self.lam)) / (np.pi*self.nu)**(self.D/2)
        temp3 = 1 + np.dot(np.dot((x-self.mu), self.lam),  (x-self.mu))/self.nu
        temp4 = -(self.nu+self.D)/2
        return temp1*temp2*((temp3)**temp4)

class mixture_gaussian_classifier:
    def __init__(self, alpha0, m0, invW0, beta0, nu0):
        self.K = alpha0.shape[0]

        #予測分布（多次元t分布）の生成
        self.pred_distr = (lambda m, W, beta, nu:
                                      multivariate_student_t
                                      (m,
                                      (1 - m.shape[0] + nu) * beta * W/(1 + beta),
                                      (1 - m.shape[0] + nu)))

        #超パラメータの初期化
        self.alpha = alpha0
        self.eta = alpha0/np.sum(alpha0)
        self.m = [None] * self.K
        self.invW = [None] * self.K
        self.beta = np.zeros(self.K)
        self.nu = np.zeros(self.K)
        for k in range(self.K):
            self.m[k] = m0
            self.invW[k] = invW0
            self.beta[k] = beta0
            self.nu[k] = nu0



    def train(self, x, s):
        #各超パラメータの更新
        k = np.where(s == 1)[0][0]
        self.alpha[k] += 1.0
        self.eta = self.alpha/np.sum(self.alpha)
        pre_beta = self.beta[k]
        self.beta[k] += 1.0
        self.nu[k] += 1.0
        pre_m = self.m[k]
        self.m[k] = (x + pre_m * pre_beta)/self.beta[k]
        self.invW[k] = (self.invW[k] + np.dot(x.reshape(-1,1), x.reshape(-1,1).T)
                        - self.beta[k] * np.dot(self.m[k].reshape(-1,1), self.m[k].reshape(-1,1).T)
                        + pre_beta * np.dot(pre_m.reshape(-1,1), pre_m.reshape(-1,1).T))

    def test(self, x):
        s = np.zeros(self.K)
        #各クラスの予測平均を算出
        for k in range(self.K):
            s[k] = (self.eta[k] *
                    self.pred_distr(self.m[k], np.linalg.inv(self.invW[k]),
                                    self.beta[k], self.nu[k]).pdf(x))
        s /= np.sum(s)
        return s

クラスmixture_gaussian_classifierは、ここまで述べてきた学習（train）と予測平均の算出（test）を行うクラスです。また、multivariate_student_t多次元のスチューデントのt分布のクラスで、この実装については、以下の記事の内容を参考にさせてもらいました。

<ベイズ推論> 多次元ガウス分布の学習 - Helve’s Python memo

実験・結果

ということで、以下のようなコードを組んで実験してみました。クラス数は $K=3$ で、学習データ数はそれぞれ $N_1=50, N_2 = 60, N_3 = 70$ と少し偏りを持たせています。
また、せっかく事後分布を逐次学習できるようにしたので、学習データが与えられるたびに予測平均がどのように変化するかをアニメーションにしてみました。

global used_idx

def update(i, idx, mgc, X, S):
    plt.cla
    test_range = range(0,60)
    X_new = np.zeros((2, len(test_range)**2))
    S_new = np.zeros((3, len(test_range)**2))
    i = idx[i]

    #データの学習
    mgc.train(X[:, i],S[:, i])

    used_idx.append(i)
    c = S[:, used_idx]

    plt.scatter(X[0, used_idx], X[1, used_idx],
                edgecolors="black", marker="o" , color=S[:, used_idx].T)

    count = 0
    S_res = np.zeros((len(test_range), len(test_range), 3))
    for x1 in test_range:
        for x2 in test_range:
            #各座標位置に対する予測平均を算出
            X_new[:, count] = np.array([float(x1), float(x2)])
            S_new[:, count] = mgc.test(X_new[:, count])
            S_res[int(x2 - test_range[0]), int(x1 - test_range[0]), :] = S_new[:, count]
            count += 1

    plt.imshow(S_res)


def main():
    N1 = 50
    N2 = 60
    N3 = 70
    X = np.zeros((2, N1+N2+N3))
    S = np.zeros((3, N1+N2+N3), 'int8')

    #クラス1の学習データ
    mu1 = np.array([45, 25])
    A = np.array([[3,0],[0,8]])
    sigma1 = np.dot(A.T, A)
    X[:, :N1] = np.random.multivariate_normal(mu1, sigma1, N1).T
    S[0,:N1] =1

    #クラス2の学習データ
    mu2 = np.array([35, 40])
    A = np.array([[5, 3],[3,7]])
    sigma2 = np.dot(A.T, A)
    X[:, N1:N1 + N2] = np.random.multivariate_normal(mu2, sigma2, N2).T
    S[1,N1:N1 + N2] = 1

    #クラス3の学習データ
    mu3 = np.array([30, 25])
    A = np.array([[7,-5],[-5,7]])
    sigma3 = np.dot(A.T, A)
    X[:, N1 + N2:]= np.random.multivariate_normal(mu3, sigma3, N3).T
    S[2,N1 + N2:] = 1

    #事前分布の超パラメータ
    alpha0 = np.random.randn(3)
    m0 = (mu1 + mu2 + mu3)/2
    invW0 = np.linalg.inv((sigma1 + sigma2+ sigma3)/3)
    beta0 = 1.0
    nu0 = m0.shape[0]

    #推論オブジェクトの生成
    mgc = mixture_gaussian_classifier(alpha0, m0,invW0, beta0, nu0)

    fig = plt.figure()
    idx = [i for i in range(X.shape[1])]
    random.shuffle(idx)

    global used_idx
    used_idx = []
    #結果のアニメーション表示（300msごとに関数updateを呼び出す）
    ani = anm.FuncAnimation(fig, update, fargs = (idx, mgc, X, S),
                            interval = 300, frames = N1 + N2+ N3)
    ani.save("result.gif", writer="imagemagick")

if __name__ == '__main__':
    main()

そして、こちらが結果です。

f:id:YamagenSakam:20190210212509g:plain

この図で、各色がそれぞれのクラスを表し、追加されていく各点が学習データを表しています。また、背景色が予測分布が算出した各座標点の各クラスへの所属確率を表しています。
特に左上に辺りの領域に着目すると、最初の方は青色のクラスと予測されていましたが、青と緑のクラスのデータを観測していくにつれ、緑色のクラスだと予測されるようになっていくのがわかります。これは青色のクラスが左下から右上にかけての共分散が大きいのに対し、緑色のクラスの共分散が左上から右下への共分散が大きく、その特徴をサンプルが増えるにつれ、正しく推論できるようになっていくためだと考えられます。

まとめ

今回、ガウス混合分布を用いた教師あり分類をベイズ推論の枠組みを適用して実装してみました。教師ありの分類の場合は、ガウス混合分布モデルによる予測は解析的に求めることができます。
ベイズ推論の枠組みで分類を行う場合は他にもロジスティック回帰をベイズの枠組みで定式化する方法などがあります。今後は、このベイズロジスティック回帰を含めもう少し複雑なモデルを学び実装していきたいと思います。

はじめに

解きたい問題

最適解

単体法の考え方

交換する変数の選択と停止、非有界の判定

初期化（2段階単体法）

処理の流れまとめ

2段階単体法の流れ

単体法の流れ

実装

最後に

はじめに

テストデータが主部分空間の元である場合

テストデータが主部分空間の元でない場合

再構成誤差の等高線

まとめ

はじめに

方法1：零空間を使う方法

方法2：正準角を使う方法

2020/2/2追記

はじめに

やったこと

実装

2020/2/20追記

結果

最後に

はじめに

方法1 ~MCMCサンプリングの時、同時分布からサンプリング~

方法2 ~からのサンプルを用いて からサンプリング~

実装

今後やりたいこと

対数尤度と下限とKLダイバージェンス

考え方

EMアルゴリズムの更新式

変分ベイズEMアルゴリズム

変分近似の考え方

変分EMアルゴリズムの考え方

変分EMアルゴリズムの更新式の導出

はじめに

ガウス混合モデル

推論

事後分布の導出

予測分布の導出

実装

実験・結果

まとめ

テストデータ ${\bf x}_{*}$ が主部分空間の元である場合

テストデータ ${\bf x}_{*}$ が主部分空間の元でない場合

方法1 ~MCMCサンプリングの時、同時分布 $p({\bf f}_N, {\bar {\bf f}}_M )$ からサンプリング~

方法2 ~ $p({\bf f}_N | {\bf y})$ からのサンプルを用いて $p({\bar {\bf f}}_M | {\bf y})$ からサンプリング~

変分ベイズ EMアルゴリズム