Hatena::ブログ(Diary)

きちめも

2010-05-19

RVMで回帰

| 22:21

5月も残り僅か...


新歓や講義も落ち着き始めたのでPRMLを再開。今回は7章後半の関連ベクトルマシン(Relevance Vector Machine)。


感想とか

  • p57に「モデル(7.78)を用いた場合は, 計画行列Φは,...」のところで、自分はN×(N+1)の行列(1列目は全て1、残りのN×NはK_nm=k(x_n, x_{m-1})な感じの行列を想像したのだが、本文を見ると(N+1)×(N+1)の対象カーネル行列と書いてある。何か勘違いしてるのかな...?
  • SVNとあんまし関係無い気g(ry
  • 本題とは外れるが(C.7)のWoodBuryの公式の便利さをとようやく理解...
  • p64の「各繰り返しにかかる計算時間はO(M^3)」ってのは、αiが無限大だとΣのi行i列はほぼ0になるので、モデルに含まれる基底ベクトルだけでΦを構成して、M×MのΣにしても計算結果あんま変わらなくて大丈夫→やった計算減るよ!(特にQとかSの) …みたいな感じで理解は合ってるのかな
  • パラメータは勝手に決めてくれるけど、カーネル関数パラメータ(ガウスカーネルとか)とかはどうすんだろ。…交差検定?

数式はちゃんと紙とペンで導出したし満足…、と言いつつほんとに動くんかいな。てなわけで以下今回書いたコード。

行列の計算が多かったので、scipyが気持ち良い。分類問題はまた明日、というか確率と画像処理論の課題ががががが。


まずは(7.87)でαの更新を行なうパターン。

αの上限値を定めたのは、値がInfになると逆行列が求められないから。

カーネル関数ガウスカーネルを用いた。

#!/usr/bin/python
# -*- coding: utf-8 -*-

import scipy as sp
import scipy.linalg as spla
import itertools as it
import functools as fn

class RVM(object):
    u"""
    relevance vector machine
    PRML 7.2.1を参考にした
    """

    def __init__(self,
                 kernel=lambda x, y: sp.exp(-sp.square(spla.norm(x-y))/0.05)):
        self._kernel = kernel

    def learn(self, X, t, tol=0.01, amax=1e10):
        u"""学習"""

        N = X.shape[0]
        a = sp.ones(N+1) # hyperparameter
        b = 1.0
        phi = sp.ones((N, N+1)) # design matrix
        phi[:,1:] = [[self._kernel(xi, xj) for xj in X] for xi in X]

        diff = 1
        while diff >= tol:
            sigma = spla.inv(sp.diag(a) + b * sp.dot(phi.T, phi))
            m = b * sp.dot(sigma, sp.dot(phi.T, t))
            gamma = sp.ones(N+1) - a * sigma.diagonal()
            anew = gamma / (m * m)
            bnew = (N -  gamma.sum()) / sp.square(spla.norm(t - sp.dot(phi, m)))
            anew[anew >= amax] = amax
            adiff, bdiff = anew - a, bnew - b
            diff = (adiff * adiff).sum() + bdiff * bdiff
            a, b = anew, bnew

        self._a = a
        self._b = b
        self._X = X
        self._m = m
        self._sigma = sigma
        self._amax = amax

    def mean(self, x):
        u"""予測値の平均"""

        ret = 0
        phi = sp.append([1.0], [self._kernel(x, xi) for xi in X])
        for i in range(len(self._m)):
            if self._a[i] < self._amax:
                ret += self._m[i] * phi[i]
        return ret

    def variance(self, x):
        u"""予測値の分散"""

        phi = sp.append([1.0], [self._kernel(x, xi) for xi in X])
        return 1.0/self._b + sp.dot(phi.T, sp.dot(self._sigma, phi))

    def _get_a(self):
        return self._a
    # 超パラメータ
    a = property(_get_a)

    def _get_rv_index(self):
        return self._a[1:] < self._amax
    # 関連ベクトルのインデックス
    rv_index = property(_get_rv_index)


if __name__ == '__main__':
    # めんどうなのでデータはここに;-p
    data = sp.array([
        [0.000000, 0.349486],
        [0.111111, 0.830839],
        [0.222222, 1.007332],
        [0.333333, 0.971507],
        [0.444444, 0.133066],
        [0.555556, 0.166823],
        [0.666667, -0.848307],
        [0.777778, -0.445686],
        [0.888889, -0.563567],
        [1.000000, 0.261502],
        ])

    X = data[:, 0]
    t = data[:, 1]
    rvm = RVM()
    rvm.learn(X, t)

    print "a=", rvm.a

    # 描画
    import matplotlib.pyplot as plt
    x = sp.linspace(0, 1, 50)

    # +-標準偏差1つ分の幅の表示
    meshx, meshy = sp.meshgrid(sp.linspace(0, 1, 200), sp.linspace(-1.5, 1.5, 200))
    meshz = [[abs(rvm.mean(meshx[j][i]) - meshy[j][i]) <= sp.sqrt(rvm.variance(meshx[j][i]))
              for i in range(len(meshx[0]))] for j in range(len(meshx))]
    plt.contour(meshx, meshy, meshz, 1)
    plt.spring()

    # 入力
    plt.scatter(X, t, label="input")

    # 関連ベクトルの描画
    plt.scatter(X[rvm.rv_index], t[rvm.rv_index], marker='d', color='r', label="relevance vector")

    # 元の曲線を表示
    y = sp.sin(2*sp.pi*x)
    plt.plot(x, y, label="sin(2pix)")

    # RVMの予測の平均
    y_  = [rvm.mean(xi) for xi in x]
    plt.plot(x, y_, label="RVM regression")

    # 表示範囲の調整
    plt.xlim(-0.05, 1.05)
    plt.ylim(-1.5, 1.5)

    # label表示
    plt.legend()
    plt.show()

出力は

a= [  1.00000000e+10   1.00000000e+10   1.00000000e+10   8.99348648e-01
   1.00000000e+10   1.00000000e+10   1.00000000e+10   1.00000000e+10
   1.95794077e+00   1.00000000e+10   1.09871526e+01]

と確かに疎。

グラフは以下。ピンクの線は平均から標準偏差1つ分に相当する幅で等高線で引いたもの。赤い点は関連ベクトル

f:id:se-kichi:20100519220138p:image

…上手く動いてるようなので、逐次的疎ベイジアン学習アルゴリズムとやらを試してみる。

#!/usr/bin/python
# -*- coding: utf-8 -*-

import scipy as sp
import scipy.linalg as spla
import itertools as it
import functools as fn


class RVM(object):
    u"""
    relevance vector machine
    PRML 7.2.2を参考にした
    """
    
    def __init__(self,
                 kernel=lambda x, y: sp.exp(-sp.square(spla.norm(x-y))/0.05)):
        self._kernel = kernel

    def learn(self, X, t, tol=0.01, amax=1e10):
        u"""学習"""

        N = X.shape[0]
        a = sp.ones(N+1)*amax # hyperparameter

        phiT = sp.ones((N+1, N)) # design matrix(*transposed*)
        phiT[1:] = [[self._kernel(xi, xj) for xj in X] for xi in X]

        bases = set([0]) # basis func indexes
        a[0] = 1.0
        b = 1

        diff = 1
        while diff >= tol:
            indexes = list(bases)
            phiT_ = phiT[indexes]
            sigma = spla.inv(sp.diag(a[indexes]) + b * sp.dot(phiT_, phiT_.T))
            m = b * sp.dot(sigma, sp.dot(phiT_, t))
            S = Q = b * phiT - b*b*sp.dot(phiT, sp.dot(phiT_.T, sp.dot(sigma, phiT_)))
            Q = sp.dot(Q, t)
            S = sp.dot(S, phiT.T).diagonal()
            q = Q / (1 - S / a)
            s = S / (1 - S / a)
            anew = a.copy()

            for i in range(N+1):
                q2 = q[i]*q[i]
                if q2 > s[i] and a[i] < amax:
                    anew[i] = s[i]*s[i]/(q2-s[i])
                elif q2 > s[i] and a[i] >= amax:
                    bases.add(i)
                    anew[i] = s[i]*s[i]/(q2-s[i])
                elif q2 <= s[i] and a[i] < amax:
                    bases.remove(i)
                    anew[i] = amax

            bnew = (sp.dot(a[indexes], sigma.diagonal())+N-len(indexes)) / sp.square(spla.norm(t - sp.dot(phiT_.T, m)))
            adiff, bdiff = anew - a, bnew - b
            diff = (adiff * adiff).sum() + bdiff * bdiff
            a, b = anew, bnew

        self._a = a
        self._b = b
        self._m = m
        self._sigma = sigma
        self._indexes = list(bases)

    def mean(self, x):
        u"""予測値の平均"""
        phi = sp.append([1.0], [self._kernel(x, xi) for xi in X])[self._indexes]
        return sp.dot(phi, self._m)

    def variance(self, x):
        u"""予測値の分散"""
        phi = sp.append([1.0], [self._kernel(x, xi) for xi in X])[self._indexes]
        return 1.0/self._b + sp.dot(phi.T, sp.dot(self._sigma, phi))

    def _get_a(self):
        return self._a
    a = property(_get_a)

    def _get_b(self):
        return self._b_
    b = property(_get_b)

    def _get_m(self):
        return self._a
    m = property(_get_m)

    def _get_rv_index(self):
        indexes = [x-1 for x in self._indexes]
        if len(indexes) > 0 and indexes[0] == -1:
            indexes = indexes[1:]
        return indexes
    # 関連ベクトルのインデックス
    rv_index = property(_get_rv_index)


if __name__ == '__main__':
    # めんどうなのでデータはここに;-p
    data = sp.array([
        [0.000000, 0.349486],
        [0.111111, 0.830839],
        [0.222222, 1.007332],
        [0.333333, 0.971507],
        [0.444444, 0.133066],
        [0.555556, 0.166823],
        [0.666667, -0.848307],
        [0.777778, -0.445686],
        [0.888889, -0.563567],
        [1.000000, 0.261502],
        ])

    X = data[:, 0]
    t = data[:, 1]
    rvm = RVM()
    rvm.learn(X, t)

    print "a=", rvm.a

    # 描画
    import matplotlib.pyplot as plt
    x = sp.linspace(0, 1, 50)

    # +-標準偏差1つ分の幅の表示
    meshx, meshy = sp.meshgrid(sp.linspace(0, 1, 200), sp.linspace(-1.5, 1.5, 200))
    meshz = [[abs(rvm.mean(meshx[j][i]) - meshy[j][i]) <= sp.sqrt(rvm.variance(meshx[j][i]))
              for i in range(len(meshx[0]))] for j in range(len(meshx))]
    plt.contour(meshx, meshy, meshz, 1)
    plt.spring()

    # 入力
    plt.scatter(X, t, label="input")

    # 関連ベクトルの描画
    plt.scatter(X[rvm.rv_index], t[rvm.rv_index], marker='d', color='r', label="relevance vector")

    # 元の曲線を表示
    y = sp.sin(2*sp.pi*x)
    plt.plot(x, y, label="sin(2pix)")

    # RVMの予測の平均
    y_  = [rvm.mean(xi) for xi in x]
    plt.plot(x, y_, label="RVM regression")

    # 表示範囲の調整
    plt.xlim(-0.05, 1.05)
    plt.ylim(-1.5, 1.5)

    # label表示
    plt.legend()
    plt.show()

出力は以下。先程と殆どかわらず。ほぇー。

a= [  1.00000000e+10   1.00000000e+10   1.00000000e+10   8.99430644e-01
   1.00000000e+10   1.00000000e+10   1.00000000e+10   1.00000000e+10
   1.96108669e+00   1.00000000e+10   1.10242579e+01]

グラフは

f:id:se-kichi:20100519220139p:image


ちなみにP59に「データが存在しない領域で予測分散が小さく...」と書いてあったので、データを一部削除して試してみた。

f:id:se-kichi:20100523165813p:image

…確かに図6.8のガウス過程さんに比べて残念っぽい。

しましましましま 2010/05/22 11:51 PRML翻訳チームのメンバーです.読んでいただいてどうもです.

> p57に「モデル(7.78)を用いた場合は, 計画行列Φは,...」のところで、自分はN×(N+1)の行列(1列目は全て1、残りのN×NはK_nm=k(x_n, x_{m-1})な感じの行列を想像したのだが、本文を見ると(N+1)×(N+1)の対象カーネル行列と書いてある。何か勘違いしてるのかな...?

確かにちょっとおかしいですね Φ^T Φ がカーネル行列とかいう意図か,他の意図か?Bishop先生に確認してみます.

> SVNとあんまし関係無い気g(ry

関係ないです.だから,章のタイトルも「疎な解を持つカーネルマシン」ですよね.
Sparse Bayesian Learning という呼び方の方がメジャーかも…
SVM とか RVM とか machine で終わらせるのは特許をとる都合という大人の事情です.

> p64の「各繰り返しにかかる計算時間はO(M^3)」ってのは、αiが無限大だとΣのi行i列はほぼ0になるので、モデルに含まれる基底ベクトルだけでΦを構成して、M×MのΣにしても計算結果あんま変わらなくて大丈夫→やった計算減るよ!(特にQとかSの) …みたいな感じで理解は合ってるのかな

O(M^3) の計算量はΣの計算で逆行列を求めるためですね.

SVM,RVM,L1正則化などで特徴を疎にしたいという要求は以下のようなところからでてきます.
- 一般に,機械学習は,役に立つ特徴だけが厳選されているよいモデル空間で学習すると汎化誤差は減る
- 画像処理とか特徴の値を計算するのが面倒だったり,医療データで特徴量をとるとき検査しないといけないとか,特徴の値を得るのに手間やコストがかかるときには,特徴数が少ないときの方がよい
- どんな特徴が分類に効いてくるかとかを見たいという目的にも使えたりする

> 超パラメータは勝手に決めてくれるけど、カーネル関数のパラメータ(ガウスカーネルとか)とかはどうすんだろ。…交差検定?

普通はそうですね.
エビデンス近似とかはできるかな?面倒そうですね…やめましょう (;^_^A


ところで numpy/scipy とかpythonの数値計算系はお詳しいですか?
私は numpy 修行中の身なのですが,演算の broadcasting とかがよく分からず,for文に手をだしてしまいます.

for x in xrange(m.nx):
for y in xrange(m.ny):
q[:, x, y] = m.pz * m.px[x, :] * m.py[y, :]

とかを1行で書こうとして挫折しました.scipy のサイトや,Python Scripting for Computational Science 本以外に良い資料をご存じでしたら教えていただけませんか?

se-kichise-kichi 2010/05/23 04:57 コメントありがとうございます、色々とスッキリしました。

> ところで numpy/scipy とかpythonの数値計算系はお詳しいですか?
全然詳しくない、というか適当に使ってます... -_-);
あまり知らないで使っているとどこかでハマりそう+冗長な書き方を気付かづにしていそう...
ということで私も何か勉強しないとと思ってる現状です;

しましましましま 2010/07/01 12:14 >> p57に「モデル(7.78)を用いた場合は, 計画行列Φは,...」のところで、自分はN×(N+1)の行列(1列目は全て1、残りのN×NはK_nm=k(x_n, x_{m-1})な感じの行列を想像したのだが、本文を見ると(N+1)×(N+1)の対象カーネル行列と書いてある。何か勘違いしてるのかな...?

> 確かにちょっとおかしいですね Φ^T Φ がカーネル行列とかいう意図か,他の意図か?Bishop先生に確認してみます.

Bishop先生よりお返事があり,次のページのP.57のよう変更することとなりました.該当文は削除です.
http://ibisforest.org/index.php?PRML%2Ferrata2

トラックバック - http://d.hatena.ne.jp/se-kichi/20100519/1274275285