ml-gaussian-process

最終更新: 2018-05-12 22:50

ガウス過程とは。関数y(x)y(\mathbf{x})上の確率分布で任意の点集合x1,,xN\mathbf{x}_1,\,\cdots,\,\mathbf{x}_Nに対するy(x)y(\mathbf{x})の値の同時分布がガウス分布に従うものである。
ガウス過程を用いて観測データから関数の事後分布を求めることで、回帰や識別を行うことができる。

例: 線形結合モデルを用いた回帰

問題設定とモデル

入力値{x1,x2,,xN}\{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}に対して出力値{t1,t2,,tN}\{t_1, t_2, \cdots, t_N\}が与えられたときに、入出力の関係を与える未知の関数ffに対する回帰問題を考える。ここで、出力値tit_iは真の関数値yi=f(xi)y_i=f(\mathbf{x}_i)に正規分布に従うノイズが加わったものとする。

ここで、関数のモデルとして

y(x)=wTϕ(x)y(\mathbf{x}) = \mathbf{w}^T \mathbf{\phi}(\mathbf{x})

という形の基底関数の線形結合を考えることにする。ここで、ϕ(x)\mathbf{\phi}(\mathbf{x})はあらかじめ固定されたMM次元の基底関数のベクトルである。w\mathbf{w}MM次元のパラメータベクトルで、w\mathbf{w}の事前分布は

p(w)=N(w0,α1I)p(\mathbf{w}) = \mathcal{N}(\mathbf{w} \mid \mathbf{0}, \alpha^{-1}\mathbf{I})

という形のガウス分布とする。

ガウス分布の線形結合はガウス分布になることからy(x)y(\mathbf{x})はガウス分布に従うことが分かる。よって、このモデルはガウス過程の例になっている。

事前分布の計算

各観測点における関数値をy=(y(x1)y(xN))\mathbf{y} =\left( \begin{array}{c} y(\mathbf{x}_1) \\ \vdots \\ y(\mathbf{x}_N) \end{array} \right) というベクトルで表すことにする。
また、Φ\PhiΦnk=ϕk(xn)\Phi_{nk} = \phi_k(\mathbf{x}_n)を要素に持つ行列とする。y=Φw\mathbf{y} = \Phi \mathbf{w}と書けることから、y\mathbf{y}の分布について

E[y]=ΦE[w]=0cov[y]=E[yyT]=ΦE[wwT]ΦT=α1ΦΦT=K\begin{aligned} \mathbb{E}[\mathbf{y}] &= \Phi \mathbb{E}[\mathbf{w}] = 0 \\ \\ \mathrm{cov}[\mathbb{y}] &= \mathbb{E}[\mathbf{y}\mathbf{y}^T] = \Phi \mathbb{E}[\mathbf{w}\mathbf{w}^T]\Phi^T = \alpha^{-1}\Phi\Phi^T = \mathbf{K} \end{aligned}

が成り立つことが分かる。
y\mathbf{y}はガウス分布に従うことから、y\mathbf{y}p(y)=N(y0,K)p(\mathbf{y}) = \mathcal{N}(\mathbf{y} \mid \mathbf{0}, \mathbf{K})というガウス分布に従うことが分かる。
なお、K\mathbf{K}Kn,m=α1ϕ(xn)Tϕ(xm)K_{n, m} = \alpha^{-1} \phi(\mathbf{x}_n)^T\phi(\mathbf{x}_m)を要素に持つ行列(グラム行列)である。
またこれ以降k(xn,xm)k(\mathbf{x}_n, \mathbf{x}_m)k(xn,xm)=α1ϕ(xn)Tϕ(xm)k(\mathbf{x}_n, \mathbf{x}_m) = \alpha^{-1}\phi(\mathbf{x}_n)^T\phi(\mathbf{x}_m)を満たすカーネル関数とする。

観測値tit_iは真の関数値にyiy_iにノイズが加わったものだから、

p(ty)=N(ty,β1IN)p(\mathbf{t} \mid \mathbf{y}) = \mathcal{N}(\mathbf{t} \mid \mathbf{y}, \beta^{-1}\mathbf{I}_N)

のように書ける。

このとき、周辺分布p(t)p(\mathbf{t})C=K+β1IN\mathbf{C}=\mathbf{K}+\beta^{-1}\mathbf{I}_Nを用いて

p(t)=N(t0,C)p(\mathbf{t}) = \mathcal{N}(\mathbf{t} \mid \mathbf{0}, \mathbf{C})

と書ける。ここでガウス分布について

p(x)=N(xμ,Λ1)p(yx)=N(yAx+b,L1)\begin{aligned} p(\mathbf{x}) &= \mathcal{N}(\mathbf{x} \mid \mu, \Lambda^{-1}) \\ p(\mathbf{y} \mid \mathbf{x}) &= \mathcal{N}(\mathbf{y} \mid A\mathbf{x}+\mathbf{b}, L^{-1}) \end{aligned}

であるとき、

p(y)=N(yAμ+b,L1+AΛ1AT)p(\mathbf{y}) = \mathcal{N}(\mathbf{y} \mid A\mu+\mathbf{b}, L^{-1}+A\Lambda^{-1}A^T)

が成立することを用いた。

事後分布の計算

計算した事前分布を用いて、NN個の入力x1,,xN\mathbf{x}_1, \cdots, \mathbf{x}_Nと観測値tN=(t1,,tN)T\mathbf{t}_N = (t_1, \cdots, t_N)^Tからなる訓練データから入力xN+1\mathbf{x}_{N+1}に対する観測値tN+1t_{N+1}の分布を求めることにする。

まず今までの議論から、tN+1\mathbf{t}_{N+1}の周辺分布はp(tN+1)=N(tN+10,CN+1)p(\mathbf{t}_{N+1}) = \mathcal{N}(\mathbf{t}_{N+1} \mid \mathbf{0}, \mathbf{C}_{N+1})と書ける。
ここでCN+1\mathbf{C}_{N+1}

CN+1=(CNkkTc)\mathbf{C}_{N+1} = \left( \begin{array}{cc} \mathbf{C}_N & \mathbf{k} \\ \mathbf{k}^T & c \end{array} \right)

という形に分解することができる。
k\mathbf{k}kn=k(xn,xN+1)k_n = k(\mathbf{x}_n, \mathbf{x}_{N+1})という要素を持つベクトルで、c=k(xN+1,xN+1)+β1c = k(\mathbf{x}_{N+1}, \mathbf{x}_{N+1})+\beta^{-1}である。

これを用いると、tN+1t_{N+1}の事後分布p(tN+1tN)p(t_{N+1} \mid \mathbf{t}_N)の平均と分散は

μ(xN+1)=kTCN1tNσ2(xN+1)=ckTCN1k\begin{aligned} \mu(\mathbf{x}_{N+1}) &= \mathbf{k}^T\mathbf{C}_N^{-1}\mathbf{t}_N \\ \\ \sigma^2(\mathbf{x}_{N+1}) &= c - \mathbf{k}^T\mathbf{C}_N^{-1}\mathbf{k} \\ \end{aligned}

と書ける。この計算には後述する分割されたガウス分布の性質を用いた。

以上より、p(tN+1tN)=N(tN+1μ(xN+1),σ(xN+1))p(t_{N+1} \mid \mathbf{t}_N) = \mathcal{N}(t_{N+1} \mid \mu(\mathbf{x}_{N+1}), \sigma(\mathbf{x}_{N+1}))というガウス分布で表すことが出来ることが分かった。

分割されたガウス分布の性質 [PRML (2.113)-(2.115)]

x\mathbf{x}p(x)=N(μ,Σ)p(\mathbf{x}) = \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma})に従うベクトルであるとして、x=(xaxb)\mathbf{x} = \left(\begin{array}{c}\mathbf{x}_a \\ \mathbf{x}_b\end{array}\right)のように分割することを考える。

この分割に対応して、平均ベクトルμ\mathbf{\mu}の分割を

μ=(μaμb)\mathbf{\mu} = \left(\begin{array}{c}\mathbf{\mu}_a \\ \mathbf{\mu}_b\end{array}\right)

共分散行列の分割を

Σ=(ΣaaΣabΣbaΣbb)\mathbf{\Sigma} = \left(\begin{array}{cc}\mathbf{\Sigma}_{aa} & \mathbf{\Sigma}_{ab} \\ \mathbf{\Sigma}_{ba} & \mathbf{\Sigma}_{bb} \end{array}\right)

のように与えることにする。

このとき、条件付き分布p(xaxb)p(\mathbf{x}_a \mid \mathbf{x}_b)はガウス分布となり、平均と共分散は

μab=μa+ΣabΣbb1(xbμb)Σab=ΣaaΣabΣbb1Σba\begin{aligned} \mathbf{\mu}_{a \mid b} &= \mathbf{\mu}_a + \Sigma_{ab}\Sigma_{bb}^{-1}(\mathbf{x}_b - \mathbf{\mu}_b) \\ \\ \mathbf{\Sigma}_{a \mid b} &= \mathbf{\Sigma}_{aa} - \Sigma_{ab}\Sigma_{bb}^{-1}\mathbf{\Sigma}_{ba} \end{aligned}

となる。[PRML (2.94)-(2.97)]

対数尤度の計算

この回帰モデルでの対数尤度関数を考える。tN\mathbf{t}_NN(t0,CN)\mathcal{N}(\mathbf{t} \mid \mathbf{0}, \mathbf{C}_N)に従うことから、パラメータをθ\mathbf{\theta}とすれば対数尤度関数は

lnp(tθ)=12lnCN12tTCN1tN2ln(2π)\ln p(\mathbf{t}\mid\mathbf{\theta})=-\frac{1}{2}\ln|\mathbf{C}_N|-\frac{1}{2}\mathbf{t}^T\mathbf{C}_N^{-1}\mathbf{t}-\frac{N}{2}\ln(2\pi)

と書ける。

これをθi\theta_iで微分すると

θilnp(tθ)=12Tr(CN1CNθi)+12tTCN1CNθiCN1t=12Tr((ααTCN1)CNθi)\begin{aligned} \frac{\partial}{\partial \theta_i} \ln p(\mathbf{t}\mid\mathbf{\theta}) &= -\frac{1}{2}\mathrm{Tr}\left(\mathbf{C}_N^{-1}\frac{\partial\mathbf{C}_N}{\partial\theta_i}\right)+\frac{1}{2}\mathbf{t}^T\mathbf{C}_N^{-1}\frac{\partial\mathbf{C}_N}{\partial\theta_i}\mathbf{C}_N^{-1}\mathbf{t} \\ \\ &=\frac{1}{2} \mathrm{Tr}\left((\mathbf{\alpha}\mathbf{\alpha}^T-\mathbf{C}_N^{-1})\frac{\partial \mathbf{C}_N}{\partial \theta_i}\right) \end{aligned}

となる。[GPML (5.9)]ここでα\alphaα=CN1t\mathbf{\alpha}=\mathbf{C}_N^{-1}\mathbf{t}を満たすベクトルである。

この勾配を用いて勾配法でパラメータの最適化を行うことができる。

アルゴリズム

以上の理論から、tN+1t_{N+1}の事後分布p(tN+1tN)p(t_{N+1} \mid \mathbf{t}_N)の平均と分散は

μ(xN+1)=kTCN1tNσ2(xN+1)=ckTCN1k\begin{aligned} \mu(\mathbf{x}_{N+1}) &= \mathbf{k}^T\mathbf{C}_N^{-1}\mathbf{t}_N \\ \\ \sigma^2(\mathbf{x}_{N+1}) &= c - \mathbf{k}^T\mathbf{C}_N^{-1}\mathbf{k} \\ \end{aligned}

と書けることが分かったので

を計算すれば回帰を行える。

数値計算を行うときに逆行列を陽に計算しないことは常識なので、コレスキー分解を用いて以下のように計算を行う。[GPML ALgorithm 2.1](グラム行列は対称正定値なのでCN\mathbf{C}_Nにはコレスキー分解を適用することができる)

  1. L:=cholesky(CN)L := \mathrm{cholesky}(\mathbf{C}_N)
  2. α:=LT\(L\tN)\mathbf{\alpha} := L^T\backslash(L\backslash\mathbf{t}_N)
  3. μ(xN+1):=kTα\mu(\mathbf{x}_{N+1}) : = \mathbf{k}^T\mathbf{\alpha}
  4. v:=L\k\mathbf{v} := L \backslash \mathbf{k}
  5. σ2(xN+1):=cvTv\sigma^2(\mathbf{x}_{N+1}) := c - \mathbf{v}^T\mathbf{v}
  6. logp(tθ):=12tTαilogLiin2log2π\log p(\mathbf{t} \mid \mathbf{\theta}) := -\frac{1}{2}\mathbf{t}^T\mathbf{\alpha}-\sum_{i}\log L_{ii} - \frac{n}{2}\log 2\pi

なお、A\bA \backslash \mathbf{b}Ax=bA\mathbf{x} = \mathbf{b}を満たすx\mathbf{x}を表す。

対数尤度の勾配を求めるときには逆行列を回避する方法は思い浮かばなかった……。

カーネル関数について

上のアルゴリズムからも分かるように、ガウス過程においては基底関数ϕ(x)\phi(\mathbf{x})を具体的に求める必要はなく、カーネル関数だけを考えることになる。(カーネルトリック)

ガウス過程の表現力はカーネル関数の選択に依存する。よく使われるカーネルを紹介する。

ガウスカーネル

k(x,x)=exp(xx22σ2)k(\mathbf{x}, \mathbf{x}')= \exp\left(-\frac{||\mathbf{x}-\mathbf{x}'||^2}{2\sigma^2}\right)

Matern カーネル

Maternカーネルはr=xxr=||\mathbf{x}-\mathbf{x}'||として

k(r)=21νΓ(ν)(2νrl)νKν(2νrl)k(r) = \frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}r}{l}\right)^\nu K_{\nu}\left(\frac{\sqrt{2 \nu}r}{l}\right)

という形で書けるカーネルの総称である。Γ\Gammaはガンマ関数、KνK_\nuは変形ベッセル関数である。

ν=32,52\nu=\frac{3}{2}, \frac{5}{2}のときがよく使われており

kν=3/2(r)=(1+3rl)exp(3rl)kν=5/2(r)=(1+5rl+5r23l2)exp(5rl)\begin{aligned} k_{\nu=3/2}(r) &= \left(1+\frac{\sqrt{3}r}{l}\right)\exp\left(-\frac{\sqrt{3}r}{l}\right) \\ \\ k_{\nu=5/2}(r) &= \left(1+\frac{\sqrt{5}r}{l}+\frac{5r^2}{3l^2}\right)\exp\left(-\frac{\sqrt{5}r}{l}\right) \end{aligned}

という形になる。[GPML (4.17)]

実装

ガウス過程回帰 (PRML 第 6 章の図版再現)を参考にPRML6章の図と似たものを出力するC++のコードを書いた。行列演算ライブラリとしてEigenを用いた。また出力にはgnuplotを用いた。

#include <iostream>
#include <vector>
#include <cstdio>
#include <cmath>
#include <chrono>

#include "Eigen/Dense"
#include "Eigen/Cholesky"

// ガウスカーネル
double kernel(double x1, double x2, double t=1.0) {
    double r = pow(x1 - x2, 2);
    return exp(-t * r);
}

// ガウスカーネルのtに関する微分
double d_kernel(double x1, double x2, double t=1.0) {
    double r = pow(x1 - x2, 2);
    return -r * exp(-t * r);
}

int main(void) {
    // 入力データ
    const int N = 7;
    Eigen::VectorXd x_train(N), y_train(N);
    x_train << 
        0.000000,
        0.111111,
        0.222222,
        0.333333,
        0.444444,
        0.555556,
        0.666667;
    y_train <<
        0.349486,
        0.830839,
        1.007332,
        0.971507,
        0.133066,
        0.166823,
        -0.848307;

    // 対数尤度最大化によるパラメータ最適化
    double t = 1.0;
    std::vector<double> x_plot, y_pred_plot, y_true_plot, s_plot;
    const int max_trial = 100;
    auto start = std::chrono::system_clock::now();
    for(int trial=0; trial<max_trial; ++trial) {
        Eigen::MatrixXd Cn(N, N);
        const double beta = 30.0;
        for(int i=0; i<N; ++i) {
            for(int j=0; j<N; ++j) {
                Cn(i, j) = kernel(x_train(i), x_train(j), t) + (i==j ? 1.0/beta : 0);
            }
        }

        Eigen::LLT<Eigen::MatrixXd> chol(Cn);
        auto alpha = chol.solve(y_train);
        auto L = chol.matrixL();

        // 最後だけ結果を記録する
        if(trial == max_trial-1) {
            for(double p=0; p<1000; p++) {
                double x_pred = p / 1000.0;
                Eigen::VectorXd k(N);
                for(int i=0; i<N; ++i) {
                    k(i) = kernel(x_train(i), x_pred, t);
                }

                double y_pred = k.transpose() * alpha;
                auto v = L.solve(k);
                double s_pred = sqrt(kernel(x_pred, x_pred, t) + 1.0/beta - v.transpose() * v);
                double y_true = sin(2 * M_PI * x_pred);
                x_plot.push_back(x_pred);
                y_pred_plot.push_back(y_pred);
                y_true_plot.push_back(y_true);
                s_plot.push_back(s_pred);
            }
        }
        // 対数尤度の表示
        double tmp = 0;
        Eigen::MatrixXd MatL(L);
        for(int i=0; i<N; ++i) tmp += std::log(MatL(i, i));
        std::cout << "log likelihood: " 
                  << -0.5 * y_train.transpose()*alpha - tmp - 0.5 * N * std::log(2.0+M_PI) 
                  << "\t(" << "t = " << t << ")" << std::endl;
        // 対数尤度の微分の計算
        Eigen::MatrixXd dCn(N, N);
        for(int i=0; i<N; ++i) {
            for(int j=0; j<N; ++j) {
                dCn(i, j) = d_kernel(x_train(i), x_train(j), t);
            }
        }
        // 対角成分のみ計算すれば良いので高速化の余地がある
        // Cn.inverse()をコレスキー分解経由にするのは高速化にそれほど寄与しなかった
        double d = 0.1 * ((alpha * alpha.transpose() - Cn.inverse())* dCn).trace();
        // パラメータの更新
        t += d;
    }
    auto end = std::chrono::system_clock::now();

    // 経過時間をミリ秒単位で表示
    auto diff = end - start;
    std::cout << "elapsed time = "
        << std::chrono::duration_cast<std::chrono::milliseconds>(diff).count()
        << " msec."
        << std::endl;

    // gnuplotによる表示
    FILE *gp = popen("gnuplot -persist", "w");
    fprintf(gp, "set multiplot\n");
    fprintf(gp, "set nokey\n");
    fprintf(gp, "set xrange [0:1]\n");
    fprintf(gp, "set yrange [-1.4:1.4]\n");
    fprintf(gp, "plot '-' with filledcurves lc rgb \"pink\"\n");
    for(size_t i=0; i<x_plot.size(); ++i) {
        fprintf(gp, "%f\t%f\t%f\n", x_plot[i], y_pred_plot[i]-s_plot[i]*2,
                                               y_pred_plot[i]+s_plot[i]*2);
    }
    fprintf(gp, "e\n");

    fprintf(gp, "plot '-' with lines lt 1 lc rgb \"red\"\n");
    for(size_t i=0; i<x_plot.size(); ++i) {
        fprintf(gp, "%f\t%f\n", x_plot[i], y_pred_plot[i]);
    }
    fprintf(gp, "e\n");

    fprintf(gp, "plot '-' with lines lt 1 lc rgb \"dark-green\"\n");
    for(size_t i=0; i<x_plot.size(); ++i) {
        fprintf(gp, "%f\t%f\n", x_plot[i], y_true_plot[i]);
    }
    fprintf(gp, "e\n");

    fprintf(gp, "plot '-' with points pt 7 lc rgb \"blue\"\n");
    for (size_t i=0; i<N; ++i) {
        fprintf(gp, "%f\t%f\n", x_train[i], y_train[i]);
    }
    fprintf(gp, "e\n");
    fprintf(gp, "set nomultiplot\n");
    fflush(gp);
    getchar();
    fprintf(gp, "exit\n");
    pclose(gp);
}

最適化の結果t=4.94229が最適なパラメータとなり、以下のような図が出力された。

参考文献