読者です 読者をやめる 読者になる 読者になる

めもめも

このブログに記載の内容は個人の見解であり、必ずしも所属組織の立場、戦略、意見を代表するものではありません。

PRML6.4.2 Gaussian processes for regressionのメモ

基底を固定した場合の計算量

一般には、N✕N行例の逆行列計算が必要になるところが、有限個の基底を固定するとM✕M行列の逆行列計算に計算量が減少する理由を具体的な計算で確認。


Figure 6.10を再現するコード

わかりやすさ優先で冗長なコードにしてあります。

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from numpy.random import normal

np.random.seed(20160317)
num_data = 100
iter_num = 50
rate = 0.005
beta = 36.0
t0 = 1.0
eta1, eta2, eta3 = 1.0, 1.0, 1.0

def prep_data():
    x1 = normal(0.0, 0.5, num_data)
    x2 = x1 + normal(0.0, 0.05, num_data)
    x3 = normal(0.0, 0.5, num_data)
    t = np.sin(2.0*np.pi*x1) + normal(0, 1.0/np.sqrt(beta), num_data)
    return t, x1, x2, x3

def calc_mean(t, xs):
    m = np.zeros(len(xs))
    for i, x in enumerate(xs):
        k = np.zeros(num_data)
        for n in range(num_data):
            k[n] = t0 * np.exp(-0.5*eta1*(x1[n]-x)*(x1[n]-x))
        c = matrix_c()
        c_inv = np.linalg.inv(c)
        m[i] = np.dot(k, np.dot(c_inv, t))
    return m

def matrix_c():
    c = np.zeros(num_data*num_data).reshape(num_data,num_data)
    for n in range(num_data):
        for m in range(num_data):
            c[n][m] = t0 * np.exp(-0.5 * (
                              eta1*(x1[n]-x1[m])*(x1[n]-x1[m]) +
                              eta2*(x2[n]-x2[m])*(x2[n]-x2[m]) +
                              eta3*(x3[n]-x3[m])*(x3[n]-x3[m]) ))
    for n in range(num_data):
        c[n][n] += 1.0/beta
    return c

def deru_c1():
    c = matrix_c()
    for n in range(num_data):
        c[n][n] -= 1.0/beta
    for n in range(num_data):
        for m in range(num_data):
            c[n][m] *= -0.5 * (x1[n]-x1[m])*(x1[n]-x1[m])
    return c

def deru_c2():
    c = matrix_c()
    for n in range(num_data):
        c[n][n] -= 1.0/beta
    for n in range(num_data):
        for m in range(num_data):
            c[n][m] *= -0.5 * (x2[n]-x2[m])*(x2[n]-x2[m])
    return c

def deru_c3():
    c = matrix_c()
    for n in range(num_data):
        c[n][n] -= 1.0/beta
    for n in range(num_data):
        for m in range(num_data):
            c[n][m] *= -0.5 * (x3[n]-x3[m])*(x3[n]-x3[m])
    return c

# main
if __name__ == '__main__':
    t, x1, x2, x3 = prep_data()
    eta1s, eta2s, eta3s = [eta1], [eta2], [eta3]

    for i in range(iter_num):
        c = matrix_c()
        c_inv = np.linalg.inv(c)
        d1, d2, d3 = deru_c1(), deru_c2(), deru_c3()

        tmp0 = np.dot(t, np.dot(c_inv, np.dot(d1, np.dot(c_inv, t))))
        delta = -0.5 * np.trace(np.dot(c_inv, d1)) + 0.5 * tmp0
        eta1 += rate * eta1 * delta

        tmp0 = np.dot(t, np.dot(c_inv, np.dot(d2, np.dot(c_inv, t))))
        delta = -0.5 * np.trace(np.dot(c_inv, d2)) + 0.5 * tmp0
        eta2 += rate * eta2 * delta

        tmp0 = np.dot(t, np.dot(c_inv, np.dot(d3, np.dot(c_inv, t))))
        delta = -0.5 * np.trace(np.dot(c_inv, d3)) + 0.5 * tmp0
        eta3 += rate * eta3 * delta

        eta1s.append(eta1)
        eta2s.append(eta2)
        eta3s.append(eta3)
        
        print eta1, eta2, eta3

    fig = plt.figure()
    subplot = fig.add_subplot(2,1,1)
    subplot.set_yscale('log')
    subplot.plot(np.arange(0,iter_num+1), eta3s)
    subplot.plot(np.arange(0,iter_num+1), eta2s)
    subplot.plot(np.arange(0,iter_num+1), eta1s)

    subplot = fig.add_subplot(2,1,2)
    xs = np.linspace(-0.5,0.5,50)
    subplot.set_xlim(-0.5,0.5)

    subplot.scatter(x1,t)
    eta2, eta3 = 0, 0
    ms = calc_mean(t, xs)
    subplot.plot(xs, ms)
    eta1 = 1.0
    ms = calc_mean(t, xs)
    subplot.plot(xs, ms)
    plt.show()

実行結果はこちら。

詳細は、テキストを読んでいただくことにして・・・、ざっくり説明すると、次の通りです。

3つの特徴変数 (x_1,x_2,x_3) を用いて、t(x_1,x_2,x_3) の値を推定する問題です。実際には、x_1 だけが t と相関を持っており、x_2 は、x_1 に乱数を加えて擬似的な相関をもたせています。x_3 はまったく無関係な乱数です。

このデータに対して、下記のカーネル関数を用いて、ガウス過程による回帰を実施します。

 k({\mathbf x},{\mathbf x'}) = \theta_0\exp\left\{-\frac{1}{2}\sum_{i=1}^3\eta_i(x_i-x_i')^2\right\}

この時、ハイパーパラメーター (\eta_1,\eta_2,\eta_3) の最適化を勾配降下法で実施したのが上のグラフで、\eta_1(赤)と \eta_2(緑)に対して、\eta_3(青)は値が0に近くなります。

下のグラフは、\eta_2 = \eta_3 = 0 と置いて、x_1 に対する t の推定結果を図示しています。緑は、ハイパーパラメーターを最適化する前の結果で、青は、最適化した後の結果です。ハイパーパラメーターの最適化で大きく結果が変わることが分かります。